00001
00002
00003
00004
00005
00006
00007
00008
00009
00010 #include <stdio.h>
00011 #include <stdlib.h>
00012 #include <math.h>
00013 #include "constrel.h"
00014 #include "C30baroghel.h"
00015 #include "globalt.h"
00016
00017 C30barmat::C30barmat()
00018 {
00019 mw = 18.01528;
00020 ma = 28.9645;
00021 gasr = 8314.41;
00022
00023 t0 = 273.15;
00024 t00 = 293.15;
00025 p0 = 101325.0;
00026 tcr = 647.3;
00027
00028 w1 = 100.0;
00029 c1 = 200.0;
00030
00031
00032 phi0 = 0.1368;
00033 aphi = 7.8e-5;
00034
00035
00036 k0 = 3.2e-18;
00037 ak = 5.0e-3;
00038
00039
00040 scr = 1.0;
00041 ag = 1.0;
00042
00043
00044 sir = 0.2;
00045 aw = 2.0;
00046 bw = 6.0;
00047
00048
00049 rhos = 2625.8;
00050
00051 fs = 1.0;
00052
00053
00054 lambdas0 = 2.0;
00055 alam = -1.017e-3;
00056
00057
00058 ac = 0.35;
00059 cps0 = 940.0;
00060
00061
00062 hydren = 0.5e+6;
00063
00064 finv = 0.4;
00065
00066 fste = 0.24;
00067
00068
00069 dld = 4.0E-2;
00070
00071
00072 emod0 = 3.0e+10;
00073
00074
00075 vcoeff = 0.2;
00076
00077
00078 betas = 9.0e-6;
00079
00080
00081 alpha = 0.5;
00082
00083
00084
00085
00086 at = 0.8;
00087 bt = 20000.0;
00088
00089 acc = 1.5;
00090
00091 ddbw0 = 1.0e-20;
00092 }
00093
00094 C30barmat::~C30barmat()
00095 {}
00096
00097
00098
00099
00100
00101
00102
00103
00104 double C30barmat::sat(double pc,double t)
00105 {
00106 double sw,n,a,b,tt,q0,q2,q3,z,e,e0,g;
00107
00108 b = 2.27;
00109 q3 = 18.62e6;
00110 q2 = 7.0e6;
00111 n = 1.2;
00112 z = 0.5;
00113
00114 if (t <= 373.15)
00115 a = q3;
00116 else{
00117 tt=(t-373.15)/(647.15-373.15);
00118 q0=(q3-q2)*((2.*(tt*tt*tt))-(3.0*(tt*tt))+1.0);
00119 q2 = 25.0e6;
00120 a = q0 + q2;
00121 }
00122
00123 if (t <= tcr)
00124 e = pow(((tcr - t0)/(tcr - t)),n);
00125 else{
00126 e0=pow(((tcr-293.15)/z),n);
00127 e = n/z*e0*t + e0 - n/z*e0*(tcr - z);
00128 }
00129
00130 g = pow((e/a*pc),(b/(b-1.0)));
00131 sw = pow((g + 1.0),(-1.0/b));
00132
00133 return(sw);
00134
00135 }
00136
00137
00138
00139
00140
00141
00142
00143
00144 double C30barmat::dsat_dpc(double pc,double t)
00145 {
00146 double dsw_dpc;
00147 double n,a,b,dg_dpc,tt,q0,q2,q3,z,e,e0,g;
00148
00149 b = 2.27;
00150 q3 = 18.62e6 ;
00151 q2 = 7.0e6;
00152 n = 1.2;
00153 z = 0.5;
00154
00155 if (t <= 373.15)
00156 a = q3;
00157 else{
00158 tt=(t-373.15)/(647.15-373.15);
00159 q0=(q3-q2)*((2.*(tt*tt*tt))-(3.0*(tt*tt))+1.0);
00160 q2 = 25.0e6;
00161 a = q0 + q2;
00162 }
00163
00164 if (t <= tcr){
00165 e = pow(((tcr - t0)/(tcr - t)),n);
00166 }
00167 else{
00168 e0=pow(((tcr-293.15)/z),n);
00169 e = n/z*e0*t + e0 - n/z*e0*(tcr - z);
00170 }
00171 g = pow((e/a*pc),(b/(b-1.0)));
00172 dg_dpc = (b/(b-1.0))*e/a*pow((e/a*pc),(b/(b-1.0)-1.0));
00173 dsw_dpc = (-1.0/b)*pow((g + 1.0),(-1.0-1.0/b))*dg_dpc;
00174
00175 return(dsw_dpc);
00176 }
00177
00178
00179
00180
00181
00182
00183
00184
00185 double C30barmat::dsat_dt(double pc,double t)
00186 {
00187 double dsw_dt,e,b,n,de_dt,g,a,da_dt,dea_dt,dg_dt,tt,dtt_dt,q0,dq0_dt,q2,q3,e0,z;
00188
00189 b = 2.27;
00190 q3 = 18.62e6 ;
00191 q2 = 7.0e6;
00192 n = 1.2;
00193 z = 0.5;
00194
00195 if (t <= 373.15){
00196 a = q3;
00197 da_dt = 0.0;
00198 }
00199 else{
00200 tt=(t-373.15)/(647.15-373.15);
00201 q0=(q3-q2)*((2.*(tt*tt*tt))-(3.0*(tt*tt))+1.0);
00202 dtt_dt = 1.0/(647.15-373.15);
00203 dq0_dt = (q3-q2)*(6.0*tt*tt*dtt_dt - 6.0*tt*dtt_dt);
00204 q2 = 25.0e6;
00205 a = q0 + q2;
00206 da_dt = dq0_dt;
00207 }
00208
00209 if (t <= tcr){
00210 e = pow(((tcr - t0)/(tcr - t)),n);
00211 de_dt = n*pow(((tcr - t0)/(tcr - t)),(n-1.0))*(tcr - t0)/(tcr - t)/(tcr - t);
00212 }
00213 else{
00214 e0=pow(((tcr-293.15)/z),n);
00215 e = n/z*e0*t + e0 - n/z*e0*(tcr - z);
00216 de_dt = n/z*e0;
00217 }
00218
00219 g = pow((e/a*pc),(b/(b-1.0)));
00220 dea_dt = (de_dt*a - da_dt*e)/a/a;
00221 dg_dt = (b/(b-1.0))*pc*pow((e/a*pc),(b/(b-1.0)-1.0))*dea_dt;
00222 dsw_dt = (-1.0/b)*pow((g + 1.0),(-1.0-1.0/b))*dg_dt;
00223
00224 return(dsw_dt);
00225 }
00226
00227
00228
00229
00230
00231
00232 double C30barmat::ssp()
00233 {
00234 return(0.55);
00235 }
00236
00237
00238
00239
00240
00241
00242
00243
00244 double C30barmat::C30bar_phi(double t)
00245 {
00246 double phi;
00247
00248 phi = phi0 + aphi*(t - t00);
00249
00250 return(phi);
00251 }
00252
00253
00254
00255
00256
00257
00258
00259
00260
00261
00262 double C30barmat::C30bar_kintr(double pc,double pg,double t,double dam)
00263 {
00264 double kintr,bk,ad;
00265 state_eq tt;
00266
00267 bk = log(5.0/3.0)/log(4.0);
00268
00269 kintr = k0*exp(log(10.0)*ak*(t-298.15))*pow((pg/p0),bk);
00270
00271
00272 ad = 4.0;
00273 kintr = kintr + k0*pow(10.0,(ad*dam));
00274
00275
00276
00277
00278 return(kintr);
00279 }
00280
00281
00282
00283
00284
00285
00286
00287
00288 double C30barmat::C30bar_krg(double pc,double t)
00289 {
00290 double krg,s;
00291
00292 s = sat(pc,t);
00293
00294 krg = 1.0 - pow((s/scr),ag);
00295
00296 return(krg);
00297 }
00298
00299
00300
00301
00302
00303
00304
00305
00306
00307 double C30barmat::C30bar_krw(double pc,double t,double rh)
00308 {
00309 double krw,help,s;
00310
00311 s = sat(pc,t);
00312
00313 if (rh < 0.75)
00314 krw = pow(((s-sir)/(1.0-sir)),aw);
00315 else{
00316 help = 1.0 + pow(((1.0-rh)/0.25),bw);
00317 krw = pow(help,-1.0)*pow(s,aw);
00318 }
00319
00320 return(krw);
00321 }
00322
00323
00324
00325
00326
00327
00328
00329
00330
00331 double C30barmat::C30bar_dd(double pc,double t)
00332 {
00333 double tau,dd;
00334
00335 tau = C30bar_tau(pc,t);
00336
00337 dd = tau*mw*ma/gasr;
00338
00339
00340 return(dd);
00341 }
00342
00343
00344
00345
00346
00347
00348
00349
00350
00351 double C30barmat::C30bar_deff(double pc,double pg,double t)
00352 {
00353
00354 double deff,dd,phi,s,cdiff;
00355 state_eq tt;
00356
00357 dd = C30bar_dd(pc,t);
00358 phi = C30bar_phi(t);
00359 s = sat(pc,t);
00360 cdiff = tt.get_cdiff(pc,pg,t);
00361
00362 deff = dd*phi*(1.0 - s)*fs*cdiff;
00363
00364 return(deff);
00365 }
00366
00367
00368
00369
00370
00371
00372
00373 double C30barmat::C30bar_cps(double t)
00374 {
00375 double cps;
00376
00377 if (t < tcr)
00378 cps = cps0 + ac*(t - 298.15);
00379 else
00380 cps = cps0 + ac*(tcr - 298.15);
00381
00382 return(cps);
00383 }
00384
00385
00386
00387
00388
00389
00390
00391
00392
00393 double C30barmat::C30bar_rhocp(double pc,double pg,double t)
00394 {
00395 double s,phi,rhos,rhocp,rhow,rhog,rhogw,cps,cpw,cpga,cpgw;
00396 state_eq tt;
00397
00398 s = sat(pc,t);
00399 phi = C30bar_phi(t);
00400 rhow = tt.get_rhow(t);
00401 rhog = tt.get_rhog(pc,pg,t);
00402 rhogw = tt.get_rhogw(pc,t);
00403 rhos = C30bar_rhos();
00404 cps = C30bar_cps(t);
00405 cpw = tt.get_cpw();
00406 cpga = tt.get_cpga();
00407 cpgw = tt.get_cpgw();
00408
00409 rhocp = (1.0-phi)*rhos*cps + phi*(s*rhow*cpw + (1.0-s)*(rhog*cpga + rhogw*(cpgw-cpga)));
00410
00411 return(rhocp);
00412 }
00413
00414
00415
00416
00417
00418
00419
00420
00421 double C30barmat::C30bar_tau(double pc,double t)
00422 {
00423 double phi,s,tau;
00424
00425 s = sat(pc,t);
00426 phi = C30bar_phi(t);
00427
00428 tau = pow(phi,(1.0/3.0))*pow((1.0-s),(7.0/3.0));
00429
00430 return(tau);
00431 }
00432
00433
00434
00435
00436
00437
00438
00439
00440
00441
00442 double C30barmat::C30bar_lambdas(double t)
00443 {
00444 double lambdas;
00445
00446 if (t < tcr)
00447 lambdas = lambdas0 + alam*(t-298.15);
00448 else
00449 lambdas = lambdas0 + alam*(tcr-298.15);
00450
00451 return(lambdas);
00452 }
00453
00454
00455
00456
00457
00458
00459
00460
00461
00462 double C30barmat::C30bar_lambdaeff(double pc,double pg,double t)
00463 {
00464 double lambdaeff,lambdas,s,phi,rhow,rhos;
00465 state_eq tt;
00466
00467 s = sat(pc,t);
00468 phi = C30bar_phi(t);
00469 rhow = tt.get_rhow(t);
00470 lambdas = C30bar_lambdas(t);
00471 rhos = C30bar_rhos();
00472
00473 lambdaeff = lambdas*(1.0 + 4.0*phi*rhow*s/(1.0-phi)/rhos);
00474
00475 return(lambdaeff);
00476 }
00477
00478
00479
00480
00481
00482
00483 double C30barmat::C30bar_rhos()
00484 {
00485 return(rhos);
00486 }
00487
00488
00489
00490
00491
00492
00493 double C30barmat::C30bar_betas()
00494 {
00495 return(betas);
00496 }
00497
00498
00499
00500
00501
00502
00503
00504
00505
00506
00507
00508
00509 double C30barmat::C30bar_hydw(double pc,double pg,double t)
00510 {
00511 double fhy,hydw;
00512
00513 if((t-t0)< 105.0){
00514 fhy = 0.0;
00515 }
00516 else{
00517 fhy = (1.0+sin(3.1416/2.0*(1.0-2.0*exp(-0.004*((t-t0)-105.0)))))/2.0;
00518 }
00519
00520 hydw = fste*finv*c1*fhy;
00521
00522 return(hydw);
00523 }
00524
00525
00526
00527
00528
00529
00530
00531
00532
00533 double C30barmat::C30bar_dehydw_dt(double pc,double pg,double t)
00534 {
00535 double dfhyt,dehydw_dt;
00536
00537 if((t-t0)< 105.0){
00538 dfhyt = 0.0;
00539 }
00540 else{
00541 dfhyt = (3.1416*0.004/2.0)*cos(3.1416/2.0*(1.0-2.0*exp(-0.004*((t-t0)-105.0))))*exp(-0.004*((t-t0)-105.0));
00542 }
00543
00544 dehydw_dt = fste*finv*c1*dfhyt;
00545
00546 return(dehydw_dt);
00547 }
00548
00549
00550
00551
00552
00553
00554
00555
00556
00557
00558 double C30barmat::C30bar_hydren(double pc,double pg,double t)
00559 {
00560 return(hydren);
00561 }
00562
00563
00564
00565
00566
00567
00568
00569
00570
00571 double C30barmat::C30bar_fste(double pc,double pg,double t)
00572 {
00573 return(fste);
00574 }
00575
00576
00577
00578
00579
00580
00581
00582
00583
00584 double C30barmat::C30bar_ddbw(double pc,double pg,double t)
00585 {
00586 double ddbw;
00587
00588 ddbw = ddbw0*exp(-t/(273.15+23));
00589 if (t > tcr)
00590 ddbw = ddbw0*exp(-tcr/(273.15+23));
00591
00592 return(ddbw);
00593 }
00594
00595
00596
00597
00598
00599
00600
00601
00602
00603
00604
00605 double C30barmat::C30bar_emod(double pc,double pg,double t)
00606 {
00607 double emod,ttc;
00608
00609 ttc = t - 273.15;
00610
00611 if (ttc <= 50.0)
00612 emod = emod0;
00613 else
00614 emod = emod0*(3.14e-6*pow(ttc,2.0)-3.77e-3*(ttc)+1.1806);
00615 if (ttc >= 600.0)
00616 emod = emod0*0.05;
00617
00618 return(emod);
00619 }
00620
00621
00622
00623
00624
00625
00626
00627
00628
00629 double C30barmat::C30bar_fct(double pc,double pg,double t)
00630 {
00631 double fct,ttc;
00632
00633 ttc = t - 273.15;
00634
00635 if (ttc <= 600.0)
00636 fct = (6.0 - 8.56e-3*ttc)*1.0e6;
00637 else
00638 fct = 0.864*1.0e6;
00639
00640 return(fct);
00641 }
00642
00643
00644
00645
00646
00647
00648 double C30barmat::C30bar_xk0()
00649 {
00650 double xk0;
00651
00652 xk0 = 1.0e-4;
00653
00654 return(xk0);
00655 }
00656
00657
00658
00659
00660
00661
00662 double C30barmat::C30bar_bcc()
00663 {
00664 double bcc;
00665
00666 bcc = 2000.0;
00667
00668 return(bcc);
00669 }
00670
00671
00672
00673
00674
00675
00676 double C30barmat::C30bar_alpha()
00677 {
00678 double alpha;
00679
00680 alpha = 0.5;
00681
00682 return(alpha);
00683 }
00684
00685
00686
00687
00688
00689
00690 double C30barmat::C30bar_nu()
00691 {
00692 return(vcoeff);
00693 }
00694
00695
00696
00697
00698
00699
00700 void C30barmat::read(XFILE *in)
00701 {}
00702
00703
00704
00705
00706
00707
00708 void C30barmat::print(FILE *out)
00709 {}
00710
00711
00712
00713
00714
00715
00716
00717
00718
00719
00720
00721
00722
00723
00724 void C30barmat::give_reqntq(long *antq)
00725 {
00726
00727 antq[scal_iso_damage-1] = 1;
00728 }