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 "C30baroghelc.h"
00015 #include "globalt.h"
00016
00017 C30barmatc::C30barmatc()
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 C30barmatc::~C30barmatc()
00095 {}
00096
00097
00098
00099
00100
00101
00102
00103
00104 double C30barmatc::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 C30barmatc::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 C30barmatc::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 C30barmatc::ssp()
00233 {
00234 return(0.55);
00235 }
00236
00237
00238
00239
00240
00241
00242
00243
00244 double C30barmatc::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 C30barmatc::C30bar_kintr(double pc,double pg,double t,double dam)
00263 {
00264 double kintr,bk;
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
00273
00274
00275
00276
00277
00278 return(kintr);
00279 }
00280
00281
00282
00283
00284
00285
00286
00287
00288 double C30barmatc::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 C30barmatc::C30bar_krw(double pc,double t,double rh)
00308 {
00309 double krw,help,s;
00310
00311 s = sat(pc,t);
00312
00313
00314 if (rh < 0.75)
00315 krw = pow(((s-sir)/(1.0-sir)),aw);
00316 else{
00317 help = 1.0 + pow(((1.0-rh)/0.25),bw);
00318 krw = pow(help,-1.0)*pow(s,aw);
00319 }
00320
00321 return(krw);
00322 }
00323
00324
00325
00326
00327
00328
00329
00330
00331
00332 double C30barmatc::C30bar_dd(double pc,double t)
00333 {
00334 double tau,dd;
00335
00336 tau = C30bar_tau(pc,t);
00337
00338 dd = tau*mw*ma/gasr;
00339
00340
00341 return(dd);
00342 }
00343
00344
00345
00346
00347
00348
00349
00350
00351
00352 double C30barmatc::C30bar_deff(double pc,double pg,double t)
00353 {
00354
00355 double deff,dd,phi,s,cdiff;
00356 state_eq tt;
00357
00358 dd = C30bar_dd(pc,t);
00359 phi = C30bar_phi(t);
00360 s = sat(pc,t);
00361 cdiff = tt.get_cdiff(pc,pg,t);
00362
00363 deff = dd*phi*(1.0 - s)*fs*cdiff;
00364
00365 return(deff);
00366 }
00367
00368
00369
00370
00371
00372
00373
00374 double C30barmatc::C30bar_cps(double t)
00375 {
00376 double cps;
00377
00378 if (t < tcr)
00379 cps = cps0 + ac*(t - 298.15);
00380 else
00381 cps = cps0 + ac*(tcr - 298.15);
00382
00383 return(cps);
00384 }
00385
00386
00387
00388
00389
00390
00391
00392
00393
00394 double C30barmatc::C30bar_rhocp(double pc,double pg,double t)
00395 {
00396 double s,phi,rhos,rhocp,rhow,rhog,rhogw,cps,cpw,cpga,cpgw;
00397 state_eq tt;
00398
00399 s = sat(pc,t);
00400 phi = C30bar_phi(t);
00401 rhow = tt.get_rhow(t);
00402 rhog = tt.get_rhog(pc,pg,t);
00403 rhogw = tt.get_rhogw(pc,t);
00404 rhos = C30bar_rhos();
00405 cps = C30bar_cps(t);
00406 cpw = tt.get_cpw();
00407 cpga = tt.get_cpga();
00408 cpgw = tt.get_cpgw();
00409
00410 rhocp = (1.0-phi)*rhos*cps + phi*(s*rhow*cpw + (1.0-s)*(rhog*cpga + rhogw*(cpgw-cpga)));
00411
00412 return(rhocp);
00413 }
00414
00415
00416
00417
00418
00419
00420
00421
00422 double C30barmatc::C30bar_tau(double pc,double t)
00423 {
00424 double phi,s,tau;
00425
00426 s = sat(pc,t);
00427 phi = C30bar_phi(t);
00428
00429 tau = pow(phi,(1.0/3.0))*pow((1.0-s),(7.0/3.0));
00430
00431 return(tau);
00432 }
00433
00434
00435
00436
00437
00438
00439
00440
00441
00442
00443 double C30barmatc::C30bar_lambdas(double t)
00444 {
00445 double lambdas;
00446
00447 if (t < tcr)
00448 lambdas = lambdas0 + alam*(t-298.15);
00449 else
00450 lambdas = lambdas0 + alam*(tcr-298.15);
00451
00452 return(lambdas);
00453 }
00454
00455
00456
00457
00458
00459
00460
00461
00462
00463 double C30barmatc::C30bar_lambdaeff(double pc,double pg,double t)
00464 {
00465 double lambdaeff,lambdas,s,phi,rhow,rhos;
00466 state_eq tt;
00467
00468 s = sat(pc,t);
00469 phi = C30bar_phi(t);
00470 rhow = tt.get_rhow(t);
00471 lambdas = C30bar_lambdas(t);
00472 rhos = C30bar_rhos();
00473
00474 lambdaeff = lambdas*(1.0 + 4.0*phi*rhow*s/(1.0-phi)/rhos);
00475
00476 return(lambdaeff);
00477 }
00478
00479
00480
00481
00482
00483
00484 double C30barmatc::C30bar_rhos()
00485 {
00486 return(rhos);
00487 }
00488
00489
00490
00491
00492
00493
00494 double C30barmatc::C30bar_betas()
00495 {
00496 return(betas);
00497 }
00498
00499
00500
00501
00502
00503
00504
00505
00506
00507
00508
00509
00510 double C30barmatc::C30bar_hydw(double pc,double pg,double t)
00511 {
00512 double fhy,hydw;
00513
00514 if((t-t0)< 105.0){
00515 fhy = 0.0;
00516 }
00517 else{
00518 fhy = (1.0+sin(3.1416/2.0*(1.0-2.0*exp(-0.004*((t-t0)-105.0)))))/2.0;
00519 }
00520
00521 hydw = fste*finv*c1*fhy;
00522
00523 return(hydw);
00524 }
00525
00526
00527
00528
00529
00530
00531
00532
00533
00534 double C30barmatc::C30bar_dehydw_dt(double pc,double pg,double t)
00535 {
00536 double dfhyt,dehydw_dt;
00537
00538 if((t-t0)< 105.0){
00539 dfhyt = 0.0;
00540 }
00541 else{
00542 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));
00543 }
00544
00545 dehydw_dt = fste*finv*c1*dfhyt;
00546
00547 return(dehydw_dt);
00548 }
00549
00550
00551
00552
00553
00554
00555
00556
00557
00558
00559 double C30barmatc::C30bar_hydren(double pc,double pg,double t)
00560 {
00561 return(hydren);
00562 }
00563
00564
00565
00566
00567
00568
00569
00570
00571
00572 double C30barmatc::C30bar_fste(double pc,double pg,double t)
00573 {
00574 return(fste);
00575 }
00576
00577
00578
00579
00580
00581
00582
00583
00584
00585 double C30barmatc::C30bar_ddbw(double pc,double pg,double t)
00586 {
00587 double ddbw;
00588
00589 ddbw = ddbw0*exp(-t/(273.15+23));
00590 if (t > tcr)
00591 ddbw = ddbw0*exp(-tcr/(273.15+23));
00592
00593 return(ddbw);
00594 }
00595
00596
00597
00598
00599
00600
00601
00602
00603
00604
00605
00606 double C30barmatc::C30bar_emod(double pc,double pg,double t)
00607 {
00608 double emod,ttc;
00609
00610 ttc = t - 273.15;
00611
00612 if (ttc <= 50.0)
00613 emod = emod0;
00614 else
00615 emod = emod0*(3.14e-6*pow(ttc,2.0)-3.77e-3*(ttc)+1.1806);
00616 if (ttc >= 600.0)
00617 emod = emod0*0.05;
00618
00619 return(emod);
00620 }
00621
00622
00623
00624
00625
00626
00627
00628
00629
00630 double C30barmatc::C30bar_fct(double pc,double pg,double t)
00631 {
00632 double fct,ttc;
00633
00634 ttc = t - 273.15;
00635
00636 if (ttc <= 600.0)
00637 fct = (6.0 - 8.56e-3*ttc)*1.0e6;
00638 else
00639 fct = 0.864*1.0e6;
00640
00641 return(fct);
00642 }
00643
00644
00645
00646
00647
00648
00649 double C30barmatc::C30bar_xk0()
00650 {
00651 double xk0;
00652
00653 xk0 = 1.0e-4;
00654
00655 return(xk0);
00656 }
00657
00658
00659
00660
00661
00662
00663 double C30barmatc::C30bar_bcc()
00664 {
00665 double bcc;
00666
00667 bcc = 2000.0;
00668
00669 return(bcc);
00670 }
00671
00672
00673
00674
00675
00676
00677 double C30barmatc::C30bar_alpha()
00678 {
00679 double alpha;
00680
00681 alpha = 0.5;
00682
00683 return(alpha);
00684 }
00685
00686
00687
00688
00689
00690
00691 double C30barmatc::C30bar_nu()
00692 {
00693 return(vcoeff);
00694 }
00695
00696
00697
00698
00699
00700
00701 void C30barmatc::read(XFILE *in)
00702 {}