![(please configure the [header_logo] section in trac.ini)](/trac/issm/chrome/common/trac_banner.png)
Changeset 26744 for issm/trunk/src/m/classes/materials.m
- Timestamp:
- 12/22/21 10:39:44 (3 years ago)
- Location:
- issm/trunk
- Files:
-
- 3 edited
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk
- Property svn:mergeinfo changed
/issm/trunk-jpl merged: 25837-25866,25868-25993,25995-26330,26332-26733,26736-26739,26741
- Property svn:mergeinfo changed
-
issm/trunk/src
- Property svn:mergeinfo changed
-
issm/trunk/src/m/classes/materials.m
r25836 r26744 37 37 self.addprop('thermalconductivity'); 38 38 self.addprop('temperateiceconductivity'); 39 self.addprop('effectiveconductivity_averaging'); 39 40 self.addprop('meltingpoint'); 40 41 self.addprop('beta'); … … 52 53 self.addprop('burgers_viscosity'); 53 54 self.addprop('burgers_mu'); 54 self.addprop('isburgers'); 55 self.addprop('ebm_alpha'); 56 self.addprop('ebm_delta'); 57 self.addprop('ebm_taul'); 58 self.addprop('ebm_tauh'); 59 self.addprop('rheologymodel'); 55 60 self.addprop('density'); 56 61 self.addprop('issolid'); … … 59 64 self.addprop('rho_water'); 60 65 self.addprop('rho_freshwater'); 61 self.addprop('earth_density');62 66 otherwise 63 67 error('materials constructor error message: nature of the material not supported yet! (''ice'' or ''litho'' or ''hydro'')'); 64 68 end 65 69 end 70 self.addprop('earth_density'); 71 66 72 %set default parameters: 67 73 self.setdefaultparameters(); … … 90 96 91 97 %ice latent heat of fusion L (J/kg) 92 self.latentheat=3.34*1 0^5;98 self.latentheat=3.34*1e5; 93 99 94 100 %ice thermal conductivity (W/m/K) … … 98 104 self.temperateiceconductivity=.24; 99 105 106 %computation of effective conductivity 107 self.effectiveconductivity_averaging=1; 108 100 109 %the melting point of ice at 1 atmosphere of pressure in K 101 110 self.meltingpoint=273.15; 102 111 103 112 %rate of change of melting point with pressure (K/Pa) 104 self.beta=9.8*1 0^-8;113 self.beta=9.8*1e-8; 105 114 106 115 %mixed layer (ice-water interface) heat capacity (J/kg/K) … … 108 117 109 118 %thermal exchange velocity (ice-water interface) (m/s) 110 self.thermal_exchange_velocity=1.00*1 0^-4;119 self.thermal_exchange_velocity=1.00*1e-4; 111 120 112 121 %Rheology law: what is the temperature dependence of B with T … … 114 123 self.rheology_law='Paterson'; 115 124 125 %Rheology fields default: 126 self.rheology_B = 1 * 1e8; 127 self.rheology_n = 3; 128 116 129 case 'litho' 117 130 %we default to a configuration that enables running GIA solutions using giacaron and/or giaivins. … … 123 136 124 137 self.viscosity=[1e21;1e40]; %mantle and lithosphere viscosity (respectively) [Pa.s] 125 self.lame_mu=[1.45*1e11;6.7*1 0^10]; % (Pa) %lithosphere and mantle shear modulus (respectively) [Pa]138 self.lame_mu=[1.45*1e11;6.7*1e10]; % (Pa) %lithosphere and mantle shear modulus (respectively) [Pa] 126 139 self.lame_lambda=self.lame_mu; % (Pa) %mantle and lithosphere lamba parameter (respectively) [Pa] 127 140 self.burgers_viscosity=[NaN;NaN]; 128 141 self.burgers_mu=[NaN;NaN]; 129 self.isburgers=[0;0]; 142 143 self.ebm_alpha=[NaN;NaN]; 144 self.ebm_delta=[NaN;NaN]; 145 self.ebm_taul=[NaN;NaN]; 146 self.ebm_tauh=[NaN;NaN]; 147 self.rheologymodel=[0;0]; 130 148 self.density=[5.51*1e3;5.50*1e3]; % (Pa) %mantle and lithosphere density [kg/m^3] 131 149 self.issolid=[1;1]; % is layer solid or liquid. … … 137 155 %ocean water density (kg/m^3) 138 156 self.rho_water=1023.; 139 % average density of the Earth (kg/m^3) 140 self.earth_density=5512; 141 157 142 158 %fresh water density (kg/m^3) 143 159 self.rho_freshwater=1000.; … … 146 162 error('materials constructor error message: nature of the material not supported yet! (''ice'' or ''litho'' or ''hydro'')'); 147 163 end 164 165 % average density of the Earth (kg/m^3) 166 self.earth_density=5512; 167 148 168 end 149 169 end % }}} … … 155 175 switch nat 156 176 case 'ice' 157 disp(sprintf(' \nIce:'));177 disp(sprintf('\n Ice:')); 158 178 fielddisplay(self,'rho_ice','ice density [kg/m^3]'); 159 179 fielddisplay(self,'rho_water','ocean water density [kg/m^3]'); … … 172 192 fielddisplay(self,'rheology_law',['law for the temperature dependance of the rheology: ''None'', ''BuddJacka'', Cuffey'', ''CuffeyTemperate'', ''Paterson'', ''Arrhenius'', ''LliboutryDuval'', ''NyeCO2'', or ''NyeH2O''']); 173 193 case 'litho' 174 disp(sprintf(' \nLitho:'));175 fielddisplay(self,'numlayers','number of layers (default 2)');194 disp(sprintf('\n Litho:')); 195 fielddisplay(self,'numlayers','number of layers (default: 2)'); 176 196 fielddisplay(self,'radius','array describing the radius for each interface (numlayers+1) [m]'); 177 197 fielddisplay(self,'viscosity','array describing each layer''s viscosity (numlayers) [Pa.s]'); … … 180 200 fielddisplay(self,'burgers_viscosity','array describing each layer''s transient viscosity, only for Burgers rheologies (numlayers) [Pa.s]'); 181 201 fielddisplay(self,'burgers_mu','array describing each layer''s transient shear modulus, only for Burgers rheologies (numlayers) [Pa]'); 182 fielddisplay(self,'isburgers','array describing whether we adopt a MaxWell (0) or Burgers (1) rheology (default 0)'); 202 203 fielddisplay(self,'ebm_alpha','array describing each layer''s exponent parameter controlling the shape of shear modulus curve between taul and tauh, only for EBM rheology (numlayers)'); 204 fielddisplay(self,'ebm_delta','array describing each layer''s amplitude of the transient relaxation (ratio between elastic rigity to pre-maxwell relaxation rigity), only for EBM rheology (numlayers)'); 205 fielddisplay(self,'ebm_taul','array describing each layer''s starting period for transient relaxation, only for EBM rheology (numlayers) [s]'); 206 fielddisplay(self,'ebm_tauh','array describing each layer''s array describing each layer''s end period for transient relaxation, only for Burgers rheology (numlayers) [s]'); 207 208 209 fielddisplay(self,'rheologymodel','array describing whether we adopt a Maxwell (0), Burgers (1) or EBM (2) rheology (default: 0)'); 183 210 fielddisplay(self,'density','array describing each layer''s density (numlayers) [kg/m^3]'); 184 211 fielddisplay(self,'issolid','array describing whether the layer is solid or liquid (default 1) (numlayers)'); 185 212 case 'hydro' 186 disp(sprintf(' \nHydro:'));213 disp(sprintf('\n Hydro:')); 187 214 fielddisplay(self,'rho_ice','ice density [kg/m^3]'); 188 215 fielddisplay(self,'rho_water','ocean water density [kg/m^3]'); … … 217 244 md = checkfield(md,'fieldname','materials.density','NaN',1,'Inf',1,'size',[md.materials.numlayers 1],'>',0); 218 245 md = checkfield(md,'fieldname','materials.viscosity','NaN',1,'Inf',1,'size',[md.materials.numlayers 1],'>=',0); 219 md = checkfield(md,'fieldname','materials. isburgers','NaN',1,'Inf',1,'size',[md.materials.numlayers 1],'>=',0,'<=',1);246 md = checkfield(md,'fieldname','materials.rheologymodel','NaN',1,'Inf',1,'size',[md.materials.numlayers 1],'>=',0,'<=',2); 220 247 md = checkfield(md,'fieldname','materials.burgers_viscosity','Inf',1,'size',[md.materials.numlayers 1],'>=',0); 221 248 md = checkfield(md,'fieldname','materials.burgers_mu','Inf',1,'size',[md.materials.numlayers 1],'>=',0); 249 md = checkfield(md,'fieldname','materials.ebm_alpha','Inf',1,'size',[md.materials.numlayers 1],'>=',0); 250 md = checkfield(md,'fieldname','materials.ebm_delta','Inf',1,'size',[md.materials.numlayers 1],'>=',0); 251 md = checkfield(md,'fieldname','materials.ebm_taul','Inf',1,'size',[md.materials.numlayers 1],'>=',0); 252 md = checkfield(md,'fieldname','materials.ebm_tauh','Inf',1,'size',[md.materials.numlayers 1],'>=',0); 222 253 223 254 for i=1:md.materials.numlayers, 224 if md.materials.isburgers(i) & (isnan(md.materials.burgers_viscosity(i) | isnan(md.materials.burgers_mu(i)))), 225 error('materials checkconsistency error message: Litho burgers_viscosity or burgers_mu has NaN values, inconsistent with isburgers choice'); 255 if md.materials.rheologymodel(i)==1 & (isnan(md.materials.burgers_viscosity(i) | isnan(md.materials.burgers_mu(i)))), 256 error('materials checkconsistency error message: Litho burgers_viscosity or burgers_mu has NaN values, inconsistent with rheologymodel choice'); 257 end 258 if md.materials.rheologymodel(i)==2 & (isnan(md.materials.ebm_alpha(i)) | isnan(md.materials.ebm_delta(i)) | isnan(md.materials.ebm_taul(i)) | isnan(md.materials.ebm_tauh(i))), 259 error('materials checkconsistency error message: Litho ebm_alpha, ebm_delta, ebm_taul or ebm_tauh has NaN values, inconsistent with rheologymodel choice'); 226 260 end 227 261 end 228 262 if md.materials.issolid(1)==0 | md.materials.lame_mu(1)==0 229 263 error('First layer must be solid (issolid(1) > 0 AND lame_mu(1) > 0). Add a weak inner core if necessary.'); 230 264 end 231 265 ind=find(md.materials.issolid==0); 232 266 if sum(ismember(diff(ind),1)>=1) %if there are at least two consecutive indices that contain issolid = 0 233 267 error(['Fluid layers detected at layers #', num2str(ind'), ', but having 2 or more adjacent fluid layers is not supported yet. Consider merging them.']) 234 268 end 235 269 … … 250 284 %1: MatdamageiceEnum 2: MatestarEnum 3: MaticeEnum 4: MatenhancediceEnum 5: MaterialsEnum 251 285 WriteData(fid,prefix,'name','md.materials.nature','data',naturetointeger(self.nature),'format','IntMat','mattype',3); 252 WriteData(fid,prefix,'name','md.materials.type','data',5,'format','Integer'); %DANGER, this can evolve if you have classes. 253 286 WriteData(fid,prefix,'name','md.materials.type','data',5,'format','Integer'); %DANGER: this can evolve if you have classes. 254 287 for i=1:length(self.nature), 255 288 nat=self.nature{i}; … … 264 297 WriteData(fid,prefix,'object',self,'class','materials','fieldname','thermalconductivity','format','Double'); 265 298 WriteData(fid,prefix,'object',self,'class','materials','fieldname','temperateiceconductivity','format','Double'); 299 WriteData(fid,prefix,'object',self,'class','materials','fieldname','effectiveconductivity_averaging','format','Integer'); 266 300 WriteData(fid,prefix,'object',self,'class','materials','fieldname','meltingpoint','format','Double'); 267 301 WriteData(fid,prefix,'object',self,'class','materials','fieldname','beta','format','Double'); … … 279 313 WriteData(fid,prefix,'object',self,'class','materials','fieldname','density','format','DoubleMat','mattype',3); 280 314 WriteData(fid,prefix,'object',self,'class','materials','fieldname','viscosity','format','DoubleMat','mattype',3); 281 WriteData(fid,prefix,'object',self,'class','materials','fieldname',' isburgers','format','DoubleMat','mattype',3);315 WriteData(fid,prefix,'object',self,'class','materials','fieldname','rheologymodel','format','DoubleMat','mattype',3); 282 316 WriteData(fid,prefix,'object',self,'class','materials','fieldname','burgers_viscosity','format','DoubleMat','mattype',3); 283 317 WriteData(fid,prefix,'object',self,'class','materials','fieldname','burgers_mu','format','DoubleMat','mattype',3); 318 WriteData(fid,prefix,'object',self,'class','materials','fieldname','ebm_alpha','format','DoubleMat','mattype',3); 319 WriteData(fid,prefix,'object',self,'class','materials','fieldname','ebm_delta','format','DoubleMat','mattype',3); 320 WriteData(fid,prefix,'object',self,'class','materials','fieldname','ebm_taul','format','DoubleMat','mattype',3); 321 WriteData(fid,prefix,'object',self,'class','materials','fieldname','ebm_tauh','format','DoubleMat','mattype',3); 322 %compute earth density compatible with our layer density distribution: 323 earth_density=0; 324 for i=1:self.numlayers, 325 earth_density=earth_density + (self.radius(i+1)^3-self.radius(i)^3)*self.density(i); 326 end 327 earth_density=earth_density/self.radius(self.numlayers+1)^3; 328 self.earth_density=earth_density; 284 329 case 'hydro' 285 330 WriteData(fid,prefix,'object',self,'class','materials','fieldname','rho_ice','format','Double'); 286 331 WriteData(fid,prefix,'object',self,'class','materials','fieldname','rho_water','format','Double'); 287 WriteData(fid,prefix,'object',self,'class','materials','fieldname','earth_density','format','Double');288 332 WriteData(fid,prefix,'object',self,'class','materials','fieldname','rho_freshwater','format','Double'); 289 290 333 otherwise 291 334 error('materials constructor error message: nature of the material not supported yet! (''ice'' or ''litho'' or ''hydro'')'); 292 335 end 293 336 end 337 WriteData(fid,prefix,'data',self.earth_density,'name','md.materials.earth_density','format','Double'); 294 338 end % }}} 295 339 function self = extrude(self,md) % {{{ … … 317 361 writejsdouble(fid,[modelname '.materials.thermalconductivity'],self.thermalconductivity); 318 362 writejsdouble(fid,[modelname '.materials.temperateiceconductivity'],self.temperateiceconductivity); 363 writejsdouble(fid,[modelname '.materials.effectiveconductivity_averaging'],self.effectiveconductivity_averaging); 319 364 writejsdouble(fid,[modelname '.materials.meltingpoint'],self.meltingpoint); 320 365 writejsdouble(fid,[modelname '.materials.beta'],self.beta); … … 333 378 writejsdouble(fid,[modelname '.materials.density'],self.density); 334 379 writejsdouble(fid,[modelname '.materials.viscosity'],self.viscosity); 335 writejsdouble(fid,[modelname '.materials. isburgers'],self.isburgers);380 writejsdouble(fid,[modelname '.materials.rheologymodel'],self.rheologymodel); 336 381 writejsdouble(fid,[modelname '.materials.burgers_viscosity'],self.burgers_viscosity); 337 382 writejsdouble(fid,[modelname '.materials.burgers_mu'],self.burgers_mu); 383 writejsdouble(fid,[modelname '.materials.ebm_alpha'],self.ebm_alpha); 384 writejsdouble(fid,[modelname '.materials.ebm_delta'],self.ebm_delta); 385 writejsdouble(fid,[modelname '.materials.ebm_taul'],self.ebm_taul); 386 writejsdouble(fid,[modelname '.materials.ebm_tauh'],self.ebm_tauh); 338 387 339 388 case 'hydro' … … 351 400 352 401 end % }}} 402 function setlitho2prem(self) % {{{ 403 %based on materials.radius and materials.issolid, replaces materials.density, materials.lame_mu, materials.lame_lambda 404 % with volumetric averages of a polynomial version of the PREM model 405 if ~strcmpi(self.nature, 'litho') 406 error('materials setlitho2prem error message: materials.nature must be ''litho''') 407 end 408 409 410 ra=self.radius(end); 411 r(1) = 0.0; r(2) = 1221.5; r(3) = 3480.0; r(4) = 3630.; 412 r(5) = 5600.0; r(6) = 5701.0; r(7) = 5771.0; r(8) = 5971.; 413 r(9) = 6151.0; r(10) = 6291.0; r(11) = 6346.6; 414 r(12) = 6356.0; r(13) = 6368.0; r(14) = 6371.; 415 416 if (r(14)*1e3 ~= ra) 417 error(['materials setlitho2prem error message: Earth surface radius in PREM r=' num2str(r(end)*1e3) ' does not match materials.radius(end)=' num2str(self.radius(end))]); 418 end 419 420 d = zeros(12,4); %polynomial coef for density 421 d(1,1) = 13.0885; d(1,3) = -8.8381 ; 422 d(2,1) = 12.5815; d(2,2) = -1.2638; d(2,3) = -3.6426; d(2,4) = -5.5281; 423 d(3,1) = 7.9565 ; d(3,2) = -6.4761; d(3,3) = 5.5283; d(3,4) = -3.0807; 424 d(4,1) = 7.9565 ; d(4,2) = -6.4761; d(4,3) = 5.5283; d(4,4) = -3.0807; 425 d(5,1) = 7.9565 ; d(5,2) = -6.4761; d(5,3) = 5.5283; d(5,4) = -3.0807; 426 d(6,1) = 5.3197 ; d(6,2) = -1.4836; 427 d(7,1) = 11.2494; d(7,2) = -8.0298; 428 d(8,1) = 7.1089 ; d(8,2) = -3.8045; 429 d(9,1) = 2.6910 ; d(9,2) = 0.6924; 430 d(10,1) = 2.6910; d(10,2) = 0.6924; 431 d(11,1) = 2.9 ; 432 d(12,1) = 2.6 ; 433 434 % ocean 435 if (~self.issolid(end)) 436 d(13,1) = 1.02 ; 437 438 % continental 439 else 440 d(13,1) = d(12,1); 441 end 442 443 p = zeros(13,4); %polynomial coef for Vp, pressure wave velocity 444 p(1,1) = 11.2622 ; p(1,3) = -6.3640; 445 p(2,1) = 11.0487 ; p(2,2) = -4.0362; p(2,3) = 4.8023; p(2,4) = -13.5732; 446 p(3,1) = 15.3891 ; p(3,2) = -5.3181; p(3,3) = 5.5242; p(3,4) = -2.5514; 447 p(4,1) = 24.952 ; p(4,2) = -40.4673; p(4,3) = 51.4832; p(4,4) = -26.6419; 448 p(5,1) = 29.2766 ; p(5,2) = -23.6027; p(5,3) = 5.5242; p(5,4) = -2.5514; 449 p(6,1) = 19.0957 ; p(6,2) = -9.8672; 450 p(7,1) = 39.7027 ; p(7,2) = -32.6166; 451 p(8,1) = 20.3926 ; p(8,2) = -12.2569; 452 p(9,1) = 4.1875 ; p(9,2) = 3.9382; 453 p(10,1) = 4.1875 ; p(10,2) = 3.9382; 454 p(11,1) = 6.8 ; 455 p(12,1) = 5.8; 456 % 457 % ocean 458 if (~self.issolid(end)) 459 p(13,1) = 1.45 ; 460 % 461 % continental 462 else 463 p(13,1) = p(12,1); 464 end 465 %---- 466 % 467 s = zeros(13,4); %polynomial coef for Vs, shear wave velocity 468 % 469 s(1,1) = 3.6678; s(1,3) = -4.4475; 470 471 s(3,1) = 6.9254; s(3,2) = 1.4672; s(3,3) = -2.0834; s(3,4) = 0.9783; 472 s(4,1) = 11.1671; s(4,2) = -13.7818; s(4,3) = 17.4575; s(4,4) = -9.2777; 473 s(5,1) = 22.3459; s(5,2) = -17.2473; s(5,3) = -2.0834; s(5,4) = 0.9783; 474 s(6,1) = 9.9839; s(6,2) = -4.9324; 475 s(7,1) = 22.3512; s(7,2) = -18.5856 ; 476 s(8,1) = 8.9496; s(8,2) = -4.4597; 477 s(9,1) = 2.1519; s(9,2) = 2.3481; 478 s(10,1) = 2.1519; s(10,2) = 2.3481; 479 s(11,1) = 3.9 ; 480 s(12,1) = 3.2 ; 481 % 482 % ocean (please don't modify) 483 if (~self.issolid(end)) 484 % 485 % continental 486 else 487 s(13,1) = s(12,1); 488 end 489 % 490 % 491 r = r*1e3; 492 493 %- handling the first layer : central sphere 494 rad = self.radius; 495 rad(1) = 0.; 496 497 for j = 1:self.numlayers 498 499 ro = 0.; 500 vp = 0.; 501 vs = 0.; 502 503 for i = 1:13 504 505 r1 = 0.; 506 r2 = 0.; 507 if ((rad(j) > r(i)) & (rad(j) <= r(i+1))) 508 if (rad(j+1) <= r(i+1)) 509 r2 = rad(j+1); 510 r1 = rad(j); 511 else 512 r2 = r(i+1); 513 r1 = rad(j); 514 end 515 elseif (rad(j) <= r(i)) 516 if ((rad(j+1) > r(i)) & (rad(j+1) <= r(i+1))) 517 r2 = rad(j+1); 518 r1 = r(i); 519 elseif (rad(j+1) > r(i+1)) 520 r2 = r(i+1); 521 r1 = r(i); 522 end 523 end 524 525 t1 = d(i,1)/3.; 526 t2 = d(i,2)/(ra*4.); 527 t3 = d(i,3)/((ra^2)*5.); 528 t4 = d(i,4)/((ra^3)*6.); 529 ro = ro + t1*(r2^3) + t2*(r2^4) + t3*(r2^5) + t4*(r2^6) - ( t1*(r1^3) + t2*(r1^4) + t3*(r1^5) + t4*(r1^6) ); 530 531 t1 = p(i,1)/3.; 532 t2 = p(i,2)/(ra*4.); 533 t3 = p(i,3)/((ra^2)*5.); 534 t4 = p(i,4)/((ra^3)*6.); 535 vp = vp + t1*(r2^3) + t2*(r2^4) + t3*(r2^5) + t4*(r2^6) - ( t1*(r1^3) + t2*(r1^4) + t3*(r1^5) + t4*(r1^6) ); 536 537 t1 = s(i,1)/3.; 538 t2 = s(i,2)/(ra*4.); 539 t3 = s(i,3)/((ra^2)*5.); 540 t4 = s(i,4)/((ra^3)*6.); 541 vs = vs + t1*(r2^3) + t2*(r2^4) + t3*(r2^5) + t4*(r2^6) - ( t1*(r1^3) + t2*(r1^4) + t3*(r1^5) + t4*(r1^6) ); 542 543 end 544 ro = ro*3 / (rad(j+1)^3-rad(j)^3); 545 vp = vp*3 /(rad(j+1)^3-rad(j)^3); 546 vs = vs*3 / (rad(j+1)^3-rad(j)^3); 547 mu = ro*vs.^2; 548 la = ro*vp.^2 - 2.*mu; 549 ro = ro*1e3; 550 la = la*1e9; 551 mu = mu*1e9; 552 553 self.density(j) = ro; 554 self.lame_lambda(j) = la; 555 self.lame_mu(j) = mu; 556 end 557 558 end % }}} 353 559 end 354 560 end … … 358 564 for i=1:length(strnat), 359 565 switch strnat{i}, 360 566 case 'damageice' 361 567 intnat(i)=1; 362 568 case 'estar'
Note:
See TracChangeset
for help on using the changeset viewer.
