00001 #include "globalt.h"
00002 #include "glasgowmat.h"
00003
00004
00005
00006
00007
00008
00009
00010 glasgowmat::glasgowmat ()
00011 {
00012 model=2;
00013 ra=287.0;
00014 rv=461.5;
00015 stef=5.67e-8;
00016
00017 rhoc = 0.0;
00018 rhos = 0.0;
00019 por0 = 0.0;
00020 k0 = 0.0;
00021 rhol0 = 0.0;
00022 sssp = 0.0;
00023 t0 = 0.0;
00024 rhov0 = 0.0;
00025 pginf = 0.0;
00026 emmi = 0.0;
00027 alph = 0.0;
00028 hq = 0.0;
00029 crhoair = 0.0;
00030 tfirestart = 0.0;
00031 }
00032
00033 glasgowmat::~glasgowmat ()
00034 {}
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048 void glasgowmat::matcond (matrix &d,long ri,long ci,long ipp)
00049 {
00050 long n;
00051 n = d.n;
00052
00053 switch (n){
00054 case 1:{
00055 matcond1d (d,ri,ci,ipp);
00056 break;
00057 }
00058 case 2:{
00059 matcond2d (d,ri,ci,ipp);
00060 break;
00061 }
00062 case 3:{
00063 matcond3d (d,ri,ci,ipp);
00064 break;
00065 }
00066 default:{
00067 fprintf (stderr,"\n unknown number of components of stress tensor is required");
00068 fprintf (stderr,"\n in function (%s, line %d).\n",__FILE__,__LINE__);
00069 }
00070 }
00071 }
00072
00073
00074
00075
00076
00077
00078
00079
00080
00081
00082 void glasgowmat::matcond1d (matrix &d,long ri,long ci,long ipp)
00083 {
00084 double k;
00085 double pg,t,rhov;
00086 k = 0.0;
00087
00088 rhov = Tm->ip[ipp].av[0];
00089 pg = Tm->ip[ipp].av[1];
00090 t = Tm->ip[ipp].av[2];
00091
00092 if((ri == 0) && (ci == 0))
00093 k = kmv (t,pg,rhov);
00094 if((ri == 0) && (ci == 1))
00095 k = kmp (t,pg,rhov);
00096 if((ri == 0) && (ci == 2))
00097 k = kmt (t,pg,rhov);
00098
00099 if((ri == 1) && (ci == 0))
00100 k = kav (t,pg,rhov);
00101 if((ri == 1) && (ci == 1))
00102 k = kap (t,pg,rhov);
00103 if((ri == 1) && (ci == 2))
00104 k = kat (t,pg,rhov);
00105
00106 if((ri == 2) && (ci == 0))
00107 k = ktv (t,pg,rhov);
00108 if((ri == 2) && (ci == 1))
00109 k = ktp (t,pg,rhov);
00110 if((ri == 2) && (ci == 2))
00111 k = ktt (t,pg,rhov);
00112
00113 d[0][0] = k;
00114 }
00115
00116
00117
00118
00119
00120
00121
00122
00123
00124 void glasgowmat::matcond2d (matrix &d,long ri,long ci,long ipp)
00125 {
00126 double k;
00127 double pg,t,rhov;
00128 k = 0.0;
00129
00130 rhov = Tm->ip[ipp].av[0];
00131 pg = Tm->ip[ipp].av[1];
00132 t = Tm->ip[ipp].av[2];
00133
00134 if((ri == 0) && (ci == 0))
00135 k = kmv (t,pg,rhov);
00136 if((ri == 0) && (ci == 1))
00137 k = kmp (t,pg,rhov);
00138 if((ri == 0) && (ci == 2))
00139 k = kmt (t,pg,rhov);
00140
00141 if((ri == 1) && (ci == 0))
00142 k = kav (t,pg,rhov);
00143 if((ri == 1) && (ci == 1))
00144 k = kap (t,pg,rhov);
00145 if((ri == 1) && (ci == 2))
00146 k = kat (t,pg,rhov);
00147
00148 if((ri == 2) && (ci == 0))
00149 k = ktv (t,pg,rhov);
00150 if((ri == 2) && (ci == 1))
00151 k = ktp (t,pg,rhov);
00152 if((ri == 2) && (ci == 2))
00153 k = ktt (t,pg,rhov);
00154
00155 fillm(0.0,d);
00156
00157 d[0][0] = k; d[0][1] = 0.0;
00158 d[1][0] = 0.0; d[1][1] = k;
00159 }
00160
00161
00162
00163
00164
00165
00166
00167
00168
00169
00170 void glasgowmat::matcond3d (matrix &d,long ri,long ci,long ipp)
00171 {
00172 double k;
00173 double pg,t,rhov;
00174 k = 0.0;
00175
00176 rhov = Tm->ip[ipp].av[0];
00177 pg = Tm->ip[ipp].av[1];
00178 t = Tm->ip[ipp].av[2];
00179
00180 if((ri == 0) && (ci == 0))
00181 k = kmv (t,pg,rhov);
00182 if((ri == 0) && (ci == 1))
00183 k = kmp (t,pg,rhov);
00184 if((ri == 0) && (ci == 2))
00185 k = kmt (t,pg,rhov);
00186
00187 if((ri == 1) && (ci == 0))
00188 k = kav (t,pg,rhov);
00189 if((ri == 1) && (ci == 1))
00190 k = kap (t,pg,rhov);
00191 if((ri == 1) && (ci == 2))
00192 k = kat (t,pg,rhov);
00193
00194 if((ri == 2) && (ci == 0))
00195 k = ktv (t,pg,rhov);
00196 if((ri == 2) && (ci == 1))
00197 k = ktp (t,pg,rhov);
00198 if((ri == 2) && (ci == 2))
00199 k = ktt (t,pg,rhov);
00200
00201 fillm(0.0,d);
00202
00203 d[0][0]=k; d[0][1]=0.0; d[0][2]=0.0;
00204 d[1][0]=0.0; d[1][1]=k; d[1][2]=0.0;
00205 d[2][0]=0.0; d[2][1]=0.0; d[2][2]=k;
00206 }
00207
00208
00209
00210
00211
00212
00213
00214
00215
00216
00217 void glasgowmat::matcap (double &c,long ri,long ci,long ipp)
00218 {
00219 double pg,t,rhov;
00220 c = 0.0;
00221
00222 rhov = Tm->ip[ipp].av[0];
00223 pg = Tm->ip[ipp].av[1];
00224 t = Tm->ip[ipp].av[2];
00225
00226 if((ri == 0) && (ci == 0))
00227 c = cmv (t,pg,rhov);
00228 if((ri == 0) && (ci == 1))
00229 c = cmp ();
00230 if((ri == 0) && (ci == 2))
00231 c = cmt (t,pg,rhov);
00232
00233 if((ri == 1) && (ci == 0))
00234 c = cav (t,pg,rhov);
00235 if((ri == 1) && (ci == 1))
00236 c = cap (t,pg,rhov);
00237 if((ri == 1) && (ci == 2))
00238 c = cat (t,pg,rhov);
00239
00240 if((ri == 2) && (ci == 0))
00241 c = ctv (t,pg,rhov);
00242 if((ri == 2) && (ci == 1))
00243 c = ctp ();
00244 if((ri == 2) && (ci == 2))
00245 c = ctt (t,pg,rhov);
00246 }
00247
00248
00249
00250
00251 void glasgowmat::read (XFILE *in)
00252 {
00253 xfscanf (in,"%lf %lf %lf %lf %lf %lf %lf %lf %lf",&rhoc,&rhos,&por0,&k0,&rhol0,&sssp,&t0,&rhov0,&pginf);
00254 xfscanf (in,"%lf %lf %lf %lf %lf",&emmi,&alph,&hq,&crhoair,&tfirestart);
00255 }
00256
00257
00258
00259 void glasgowmat::print (FILE *out)
00260 {
00261 fprintf (out," %e %e %e %e %e %e %e %e %e",rhoc,rhos,por0,k0,rhol0,sssp,t0,rhov0,pginf);
00262 fprintf (out,"%e %e %e %e %e",emmi,alph,hq,crhoair,tfirestart);
00263 }
00264
00265
00266
00267
00268
00269 double glasgowmat::f_rhovinf ()
00270 {
00271 double rhovinf;
00272
00273 rhovinf=rhov0*0.8;
00274
00275 return(rhovinf);
00276 }
00277
00278
00279
00280
00281
00282
00283 double glasgowmat::f_tinf (double time)
00284 {
00285 double tinf;
00286
00287
00288
00289
00290
00291 tinf=((t0-273.15)+345.0*log10(8.0*(time-tfirestart)/60.0+1.0))+273.15;
00292
00293 return(tinf);
00294 }
00295
00296
00297
00298
00299
00300 double glasgowmat::f_pv (double rhov, double t)
00301 {
00302 double pv;
00303
00304 pv=rhov*rv*t;
00305
00306 return(pv);
00307 }
00308
00309
00310
00311
00312
00313 double glasgowmat::f_psat(double t)
00314 {
00315 double psat;
00316
00317 psat = -0.000000001437422219446870*t*t*t*t*t*t +
00318 0.000004424390583021230*t*t*t*t*t -
00319 0.003928080821257910*t*t*t*t +
00320 1.591032529443030*t*t*t -
00321 325.8874385048470*t*t +
00322 32147.79527519750*t -
00323 1154663.603250870;
00324
00325 return(psat);
00326 }
00327
00328
00329
00330
00331
00332 double glasgowmat::f_dpsatdt (double t)
00333 {
00334 double dpsatdt;
00335
00336 dpsatdt=32147.7952751975 -
00337 651.774877009694*t +
00338 4.77309758832909*t*t -
00339 0.01571232328503164*t*t*t +
00340 0.00002212195291510615*t*t*t*t -
00341 8.62453331668122e-9*t*t*t*t*t;
00342
00343 return(dpsatdt);
00344 }
00345
00346
00347
00348
00349
00350 double glasgowmat::f_rhow (double t)
00351 {
00352 double rhow;
00353
00354
00355
00356 rhow=1000.0;
00357
00358
00359
00360
00361
00362
00363
00364
00365
00366
00367
00368
00369
00370
00371
00372
00373
00374
00375
00376
00377
00378
00379
00380
00381
00382
00383
00384
00385
00386
00387
00388
00389
00390
00391
00392
00393
00394
00395 return(rhow);
00396 }
00397
00398
00399
00400
00401
00402
00403 double glasgowmat::f_rhow0 (double t)
00404 {
00405 double rhow0;
00406
00407
00408
00409 rhow0=1000.0;
00410
00411
00412
00413
00414
00415
00416
00417
00418
00419
00420
00421
00422
00423
00424
00425
00426
00427
00428
00429
00430
00431
00432
00433
00434
00435
00436
00437
00438
00439
00440
00441
00442
00443
00444
00445
00446
00447
00448 return(rhow0);
00449 }
00450
00451
00452
00453
00454
00455 double glasgowmat::f_drhowdt (double t)
00456 {
00457 double drhowdt;
00458
00459
00460 drhowdt=0.0;
00461
00462
00463
00464
00465
00466
00467
00468
00469
00470
00471
00472
00473
00474
00475
00476
00477
00478
00479
00480
00481
00482
00483
00484
00485
00486
00487
00488
00489
00490
00491
00492
00493
00494 return(drhowdt);
00495 }
00496
00497
00498
00499
00500
00501 double glasgowmat::f_por (double t)
00502 {
00503 double tc,por,pora,porb,porc,pord;
00504
00505 tc=t-273.15;
00506
00507
00508
00509
00510
00511
00512
00513
00514
00515 pora= (-1.0/85750000.0);
00516 porb= (27.0/1715000.0);
00517 porc= (-24.0/8575.0);
00518 pord= (389.0/343.0);
00519
00520 if(tc < 100.0){
00521 por=por0;
00522 }
00523 if((tc >= 100.0) && (tc <= 800.0)){
00524 por=por0*((pora*tc*tc*tc)+(porb*tc*tc)+(porc*tc)+pord);
00525 }
00526 if(tc > 800.0){
00527 por=por0*3.0;
00528 }
00529
00530 return(por);
00531 }
00532
00533
00534
00535
00536
00537
00538
00539 double glasgowmat::f_dpordt (double t)
00540 {
00541 double dpordt,tc,pora,porb,porc,pord;
00542
00543 tc=t-273.15;
00544
00545
00546
00547
00548
00549
00550
00551
00552
00553 pora=-(1.0/85750000.0);
00554 porb=(27.0/1715000.0);
00555 porc=-(24.0/8575.0);
00556 pord=(389.0/343.0);
00557
00558 if(tc < 100.0){
00559 dpordt=0.0;
00560 }
00561 if((tc >= 100.0) && (tc <= 800.0)){
00562 dpordt=por0*((3*pora*tc*tc)+(2*porb*tc)+porc);
00563 }
00564 if(tc > 800.0){
00565 dpordt=0.0;
00566 }
00567
00568 return(dpordt);
00569 }
00570
00571
00572
00573
00574
00575
00576 double glasgowmat::f_fracl0 (double t)
00577 {
00578 double rhow0,fracl0;
00579
00580 rhow0 = f_rhow0(t);
00581
00582 fracl0=rhol0/rhow0;
00583 if(fracl0 > por0)
00584 fracl0=por0;
00585
00586 return(fracl0);
00587 }
00588
00589
00590
00591
00592
00593
00594 double glasgowmat::f_h (double rhov, double t)
00595 {
00596 double pv,psat,h;
00597
00598 pv = f_pv(rhov,t);
00599 psat = f_psat(t);
00600
00601 h=pv/psat;
00602
00603 return(h);
00604 }
00605
00606
00607
00608
00609
00610
00611 double glasgowmat::f_m (double t)
00612 {
00613 double m;
00614 double tc;
00615
00616 tc=t-273.15;
00617
00618 m=1.04-(((tc+10.0)*(tc+10.0))/((27317.5+((tc+10.0)*(tc+10.0)))));
00619
00620 if (m > 1.0){
00621 m = 1.0;
00622 }
00623
00624 return(m);
00625 }
00626
00627
00628
00629
00630
00631
00632
00633 double glasgowmat::f_mi (double t)
00634 {
00635 double m,mi;
00636 double tc;
00637
00638 tc=t-273.15;
00639
00640 m=1.04-(((tc+10.0)*(tc+10.0))/((27317.5+((tc+10.0)*(tc+10.0)))));
00641 mi=1.0/m;
00642
00643 if (m > 1.0){
00644 m = 1.0;
00645 mi = 1.0;
00646 }
00647
00648 return(mi);
00649 }
00650
00651
00652
00653
00654
00655
00656 double glasgowmat::f_dmdt (double t)
00657 {
00658 double m,dmdt;
00659 double tc;
00660
00661 tc=t-273.15;
00662
00663 m=1.04-(((tc+10.0)*(tc+10.0))/((27317.5+((tc+10.0)*(tc+10.0)))));
00664 dmdt= -1.0*((54635.0*(10.0+tc)))/(27417.5 + 20.0*tc + tc*tc)/(27417.5 + 20.0*tc + tc*tc);
00665
00666
00667 if (m > 1.0){
00668 dmdt = 0.0;
00669 }
00670
00671 return(dmdt);
00672 }
00673
00674
00675
00676
00677
00678
00679 double glasgowmat::f_dmidt (double t)
00680 {
00681 double m,dmdt,dmidt;
00682 double tc;
00683
00684 tc=t-273.15;
00685
00686 m=1.04-(((tc+10.0)*(tc+10.0))/((27317.5+((tc+10.0)*(tc+10.0)))));
00687 dmdt= -1.0*((54635.0*(10.0+tc)))/(27417.5 + 20.0*tc + tc*tc)/(27417.5 + 20.0*tc + tc*tc);
00688
00689 dmidt = -1.0/m/m*dmdt;
00690
00691
00692 if (m > 1.0){
00693 dmidt = 0.0;
00694 }
00695
00696 return(dmidt);
00697 }
00698
00699
00700
00701
00702
00703 double glasgowmat::f_fracl (double rhov, double t)
00704 {
00705 double fracl0,rhow0,fracl,h;
00706 double aaa,bbb,ccc,ddd;
00707 double mi,rhow,m,pv,psat,por;
00708
00709 h = f_h(rhov,t);
00710 fracl0 = f_fracl0(t);
00711 rhow0 = f_rhow0(t);
00712 mi = f_mi(t);
00713 rhow = f_rhow(t);
00714 m = f_m(t);
00715 pv = f_pv(rhov,t);
00716 psat = f_psat(t);
00717 por = f_por(t);
00718
00719
00720
00721 if(h <= 0.96){
00722
00723 fracl=exp(mi*log((fracl0*rhow0/rhoc)*h))*rhoc/rhow;
00724
00725 }
00726 if((h > 0.96) && (h < 1.04)){
00727
00728
00729
00730
00731
00732
00733
00734
00735
00736 */
00737
00738
00739
00740
00741
00742
00743
00744
00745
00746
00747
00748
00749
00750
00751
00752
00753 aaa=(-3887.499999985135*fracl0*rhow0*m + 3906.249999985064*exp(mi*log(0.96))*rhoc*
00754 exp(mi*log((fracl0*rhow0)/rhoc))*(0.04166666666666362 + m))/(m*rhow);
00755
00756 bbb=(11663.249999955411*fracl0*rhow0*m - 11718.749999955196*exp(mi*log(0.96))*rhoc*
00757 exp(mi*log((fracl0*rhow0)/rhoc))*(0.042222222222221294 + m))/(m*rhow);
00758
00759 ccc=(-11645.279999955485*fracl0*rhow0*m + 11699.999999955271*exp(mi*log(0.96))*rhoc*
00760 exp(mi*log((fracl0*rhow0)/rhoc))*(0.042824074074075444 + m))/(m*rhow);
00761
00762 ddd=(3870.02879998521*fracl0*rhow0*m - 3886.9999999851375*exp(mi*log(0.96))*rhoc*
00763 exp(mi*log((fracl0*rhow0)/rhoc))*(0.043478260869569095 + m))/(m*rhow);
00764
00765 fracl = aaa*(pv/psat)*(pv/psat)*(pv/psat) + bbb*(pv/psat)*(pv/psat) + ccc*(pv/psat) + ddd;
00766
00767
00768 }
00769 if(h >= 1.04){
00770 fracl=(fracl0*rhow0*(1.0+(0.12*(h-1.04))))/rhow;
00771 }
00772
00773
00774 if(fracl > por)
00775 fracl=por;
00776
00777 return(fracl);
00778 }
00779
00780
00781
00782
00783
00784 double glasgowmat::f_dfracldt (double rhov, double t)
00785 {
00786 double dfracldt;
00787 double fracl0,rhow0,ff,ffc,h;
00788 double aaa,bbb,ccc,daaadt,dbbbdt,dcccdt,dddddt;
00789 double a1,a2,a3,b1,b2,b3,c1,c2,c3,d1,d2,d3;
00790 double rhow,m,mi,pv,psat;
00791 double drhowdt,dmdt,dmidt,dpsatdt;
00792 double derf,derff,help;
00793
00794 h = f_h(rhov,t);
00795 fracl0 = f_fracl0(t);
00796 rhow0 = f_rhow0(t);
00797 rhow = f_rhow(t);
00798 m = f_m(t);
00799 mi = f_mi(t);
00800 pv = f_pv(rhov,t);
00801 psat = f_psat(t);
00802
00803 drhowdt = f_drhowdt(t);
00804 dmdt = f_dmdt(t);
00805 dmidt = f_dmidt(t);
00806 dpsatdt = f_dpsatdt(t);
00807
00808 a1 = -3887.499999985135;
00809 a2 = 3906.249999985064;
00810 a3 = 0.04166666666666362;
00811
00812 b1 = 11663.249999955411;
00813 b2 = -11718.749999955196;
00814 b3 = 0.042222222222221294;
00815
00816 c1 = -11645.279999955485;
00817 c2 = 11699.999999955271;
00818 c3 = 0.042824074074075444;
00819
00820 d1 = 3870.02879998521;
00821 d2 = -3886.9999999851375;
00822 d3 = 0.043478260869569095;
00823
00824 ff = fracl0*rhow0;
00825 ffc= fracl0*rhow0/rhoc;
00826
00827
00828
00829 if(h <= 0.96){
00830
00831
00832
00833
00834
00835
00836
00837
00838
00839
00840
00841
00842
00843
00844
00845 derf = ffc*rhov*rv/psat - ffc*pv/psat/psat*dpsatdt;
00846 derff = (exp(mi*log(ffc*pv/psat)))*(dmidt*log(ffc*pv/psat) +
00847 mi*derf/(ffc*pv/psat));
00848 dfracldt = derff*rhoc/rhow - exp(mi*log(ffc*pv/psat))*rhoc/rhow/rhow*drhowdt;
00849 }
00850 if((h > 0.96) && (h < 1.04)){
00851
00852
00853
00854
00855
00856
00857
00858
00859
00860
00861
00862
00863
00864 */
00865
00866
00867
00868
00869
00870
00871
00872
00873
00874
00875
00876
00877
00878 aaa=(-3887.499999985135*fracl0*rhow0*m + 3906.249999985064*exp(mi*log(0.96))*rhoc*
00879 exp(mi*log((fracl0*rhow0)/rhoc))*(0.04166666666666362 + m))/(m*rhow);
00880
00881 bbb=(11663.249999955411*fracl0*rhow0*m - 11718.749999955196*exp(mi*log(0.96))*rhoc*
00882 exp(mi*log((fracl0*rhow0)/rhoc))*(0.042222222222221294 + m))/(m*rhow);
00883
00884 ccc=(-11645.279999955485*fracl0*rhow0*m + 11699.999999955271*exp(mi*log(0.96))*rhoc*
00885 exp(mi*log((fracl0*rhow0)/rhoc))*(0.042824074074075444 + m))/(m*rhow);
00886
00887
00888
00889
00890
00891
00892
00893
00894
00895
00896
00897
00898
00899
00900
00901
00902
00903
00904
00905
00906
00907
00908
00909
00910
00911
00912
00913
00914
00915
00916
00917
00918
00919
00920
00921 daaadt = -1.0*a1*ff/rhow/rhow*drhowdt + a2*a3*rhoc*(exp(mi*log(0.96))*log(0.96)*dmidt*exp(mi*log(ffc))/m/rhow +
00922 exp(mi*log(0.96))*exp(mi*log(ffc))*log(ffc)*dmidt/m/rhow +
00923 -1.0*exp(mi*log(0.96))*exp(mi*log(ffc))/m/m/rhow*dmdt +
00924 -1.0*exp(mi*log(0.96))*exp(mi*log(ffc))/m/rhow/rhow*drhowdt) +
00925 a2*rhoc*(exp(mi*log(0.96))*log(0.96)*dmidt*exp(mi*log(ffc))/rhow + exp(mi*log(0.96))*exp(mi*log(ffc))*log(ffc)
00926 *dmidt/rhow +
00927 -1.0*exp(mi*log(0.96))*exp(mi*log(ffc))/rhow/rhow*drhowdt);
00928
00929
00930 dbbbdt = -1.0*b1*ff/rhow/rhow*drhowdt + b2*b3*rhoc*(exp(mi*log(0.96))*log(0.96)*dmidt*exp(mi*log(ffc))/m/rhow +
00931 exp(mi*log(0.96))*exp(mi*log(ffc))*log(ffc)*dmidt/m/rhow +
00932 -1.0*exp(mi*log(0.96))*exp(mi*log(ffc))/m/m/rhow*dmdt +
00933 -1.0*exp(mi*log(0.96))*exp(mi*log(ffc))/m/rhow/rhow*drhowdt) +
00934 b2*rhoc*(exp(mi*log(0.96))*log(0.96)*dmidt*exp(mi*log(ffc))/rhow + exp(mi*log(0.96))*exp(mi*log(ffc))*log(ffc)
00935 *dmidt/rhow +
00936 -1.0*exp(mi*log(0.96))*exp(mi*log(ffc))/rhow/rhow*drhowdt);
00937
00938
00939 dcccdt = -1.0*c1*ff/rhow/rhow*drhowdt + c2*c3*rhoc*(exp(mi*log(0.96))*log(0.96)*dmidt*exp(mi*log(ffc))/m/rhow +
00940 exp(mi*log(0.96))*exp(mi*log(ffc))*log(ffc)*dmidt/m/rhow +
00941 -1.0*exp(mi*log(0.96))*exp(mi*log(ffc))/m/m/rhow*dmdt +
00942 -1.0*exp(mi*log(0.96))*exp(mi*log(ffc))/m/rhow/rhow*drhowdt) +
00943 c2*rhoc*(exp(mi*log(0.96))*log(0.96)*dmidt*exp(mi*log(ffc))/rhow + exp(mi*log(0.96))*exp(mi*log(ffc))*log(ffc)
00944 *dmidt/rhow +
00945 -1.0*exp(mi*log(0.96))*exp(mi*log(ffc))/rhow/rhow*drhowdt);
00946
00947
00948 dddddt = -1.0*d1*ff/rhow/rhow*drhowdt + d2*d3*rhoc*(exp(mi*log(0.96))*log(0.96)*dmidt*exp(mi*log(ffc))/m/rhow +
00949 exp(mi*log(0.96))*exp(mi*log(ffc))*log(ffc)*dmidt/m/rhow +
00950 -1.0*exp(mi*log(0.96))*exp(mi*log(ffc))/m/m/rhow*dmdt +
00951 -1.0*exp(mi*log(0.96))*exp(mi*log(ffc))/m/rhow/rhow*drhowdt) +
00952 d2*rhoc*(exp(mi*log(0.96))*log(0.96)*dmidt*exp(mi*log(ffc))/rhow + exp(mi*log(0.96))*exp(mi*log(ffc))*log(ffc)
00953 *dmidt/rhow +
00954 -1.0*exp(mi*log(0.96))*exp(mi*log(ffc))/rhow/rhow*drhowdt);
00955
00956
00957
00958
00959
00960
00961
00962
00963
00964
00965
00966
00967
00968
00969
00970
00971
00972
00973
00974
00975
00976
00977
00978
00979
00980
00981
00982
00983
00984
00985
00986
00987
00988 help = (rhov*rv*psat - pv*dpsatdt)/psat/psat;
00989 dfracldt = daaadt*(pv/psat)*(pv/psat)*(pv/psat) + dbbbdt*(pv/psat)*(pv/psat) + dcccdt*(pv/psat) + dddddt;
00990 dfracldt = dfracldt + 3.0*aaa*(pv/psat)*(pv/psat)*help + 2.0*bbb*(pv/psat)*help + ccc*help;
00991
00992
00993 }
00994 if(h >= 1.04){
00995
00996
00997
00998
00999 help = 0.12*(rhov*rv*psat - pv*dpsatdt)/psat/psat;
01000 dfracldt = fracl0*rhow0*(help*rhow - (1.0+(0.12*(h+1.04)))*drhowdt)/rhow/rhow;
01001 }
01002
01003
01004
01005
01006
01007
01008
01009
01010
01011
01012
01013 return(dfracldt);
01014 }
01015
01016
01017
01018
01019
01020 double glasgowmat::f_dfracldrhov (double rhov, double t)
01021 {
01022 double dfracldrhov;
01023 double fracl0,rhow0,h;
01024 double rhow,m,pv,psat;
01025 double aaa,bbb,ccc;
01026
01027 h = f_h(rhov,t);
01028 fracl0 = f_fracl0(t);
01029 rhow0 = f_rhow0(t);
01030 rhow = f_rhow(t);
01031 m = f_m(t);
01032 pv = f_pv(rhov,t);
01033 psat = f_psat(t);
01034
01035 if(h <= 0.96){
01036
01037
01038
01039 dfracldrhov=(fracl0*rhow0*exp((-1.0 + 1.0/m)*log((fracl0*rhow0*pv)/(rhoc*psat)))*rv*t)/(m*psat*rhow);
01040 }
01041 if((h > 0.96) && (h < 1.04)){
01042
01043
01044
01045
01046
01047
01048
01049 */
01050
01051
01052
01053
01054
01055
01056
01057
01058
01059
01060
01061
01062
01063 aaa=(-3887.499999985135*fracl0*rhow0*m + 3906.249999985064*exp((1.0/m)*log(0.96))*rhoc*
01064 exp((1.0/m)*log((fracl0*rhow0)/rhoc))*(0.04166666666666362 + m))/(m*rhow);
01065
01066 bbb=(11663.249999955411*fracl0*rhow0*m - 11718.749999955196*exp((1.0/m)*log(0.96))*rhoc*
01067 exp((1.0/m)*log((fracl0*rhow0)/rhoc))*(0.042222222222221294 + m))/(m*rhow);
01068
01069 ccc=(-11645.279999955485*fracl0*rhow0*m + 11699.999999955271*exp((1.0/m)*log(0.96))*rhoc*
01070 exp((1.0/m)*log((fracl0*rhow0)/rhoc))*(0.042824074074075444 + m))/(m*rhow);
01071
01072
01073
01074 dfracldrhov = 3.0*aaa*(pv/psat)*(pv/psat/psat)*(rv*t) + 2.0*bbb*(pv/psat/psat)*(rv*t) + ccc*(rv*t/psat);
01075
01076
01077 }
01078 if(h >= 1.04){
01079
01080 dfracldrhov=(0.12*fracl0*rhow0*rv*t)/(psat*rhow);
01081 }
01082 return(dfracldrhov);
01083 }
01084
01085
01086
01087
01088
01089
01090 double glasgowmat::f_mcbwrel (double t)
01091 {
01092 double tc,mcbwrel,mcbw0;
01093
01094 tc=t-273.15;
01095
01096
01097
01098
01099
01100
01101
01102 mcbw0=0.09*rhoc;
01103 if(tc <= 200.0){
01104 mcbwrel=0.0;
01105 }
01106 if((tc > 200.0) && (tc <= 300.0)){
01107 mcbwrel=rhoc*(7.0e-4*(tc-200.0));
01108 }
01109 if((tc > 300.0) && (tc <= 800.0)){
01110 mcbwrel=rhoc*(0.4e-4*(tc-300.0)+0.07);
01111 }
01112 if(tc > 800.0){
01113 mcbwrel=rhoc*0.09;
01114 }
01115
01116
01117
01118
01119
01120
01121
01122
01123
01124
01125
01126
01127
01128
01129
01130
01131
01132
01133
01134
01135
01136
01137
01138
01139
01140
01141
01142
01143
01144
01145
01146
01147 return(mcbwrel);
01148 }
01149
01150
01151
01152
01153
01154 double glasgowmat::f_dmcbwreldt (double t)
01155 {
01156 double tc,dmcbwreldt,mcbw0;
01157
01158 tc=t-273.15;
01159
01160
01161
01162
01163
01164
01165 mcbw0=0.09*rhoc;
01166 if(tc <= 200.0){
01167 dmcbwreldt=0.0;
01168 }
01169 if((tc > 200.0) && (tc <= 300.0)){
01170 dmcbwreldt=rhoc*7.0e-4;
01171 }
01172 if((tc > 300.0) && (tc <= 800.0)){
01173 dmcbwreldt=rhoc*0.4e-4;
01174 }
01175 if(tc > 800.0){
01176 dmcbwreldt=0.0;
01177 }
01178
01179
01180
01181
01182
01183
01184
01185
01186
01187
01188
01189
01190
01191
01192
01193
01194
01195
01196
01197
01198
01199
01200
01201
01202
01203
01204
01205
01206
01207
01208
01209
01210 return(dmcbwreldt);
01211 }
01212
01213
01214
01215
01216
01217 double glasgowmat::f_fracd (double t)
01218 {
01219 double fracd,mcbwrel,rhow;
01220
01221 mcbwrel = f_mcbwrel(t);
01222 rhow = f_rhow(t);
01223 fracd=mcbwrel/rhow;
01224
01225 return(fracd);
01226 }
01227
01228
01229
01230
01231
01232 double glasgowmat::f_dfracddt(double t)
01233 {
01234 double dfracddt,dmcbwreldt,rhow;
01235
01236 dmcbwreldt = f_dmcbwreldt(t);
01237 rhow = f_rhow(t);
01238
01239 dfracddt=dmcbwreldt/rhow;
01240
01241 return(dfracddt);
01242 }
01243
01244
01245
01246
01247
01248 double glasgowmat::f_fracg(double rhov, double t)
01249 {
01250 double fracg,por,fracl;
01251
01252 por = f_por(t);
01253 fracl = f_fracl(rhov,t);
01254
01255 fracg=por-fracl;
01256
01257 return(fracg);
01258 }
01259
01260
01261
01262
01263
01264 double glasgowmat::f_s(double rhov, double t)
01265 {
01266 double s,fracl,por;
01267
01268 por = f_por(t);
01269 fracl = f_fracl(rhov,t);
01270
01271 s=fracl/por;
01272
01273 return(s);
01274 }
01275
01276
01277
01278
01279
01280 double glasgowmat::f_pl(double rhov, double pg, double t)
01281 {
01282 double pl,s,pc;
01283
01284 s = f_s(rhov,t);
01285 pc = f_pc(rhov,t);
01286
01287 switch (model){
01288 case 1:{
01289 pl=pg;
01290 break;
01291 }
01292 case 2:{
01293 if(s <= sssp){
01294 pl = 0.0;
01295 }
01296 else{
01297 pl = pg - pc;
01298 }
01299 break;
01300 }
01301 default:{
01302 fprintf (stderr,"\n\n unknown model type in glasgowmat.cpp (%s, line %d).\n",__FILE__,__LINE__);
01303 }
01304 }
01305 return(pl);
01306 }
01307
01308
01309
01310
01311
01312 double glasgowmat::f_pc (double rhov, double t)
01313 {
01314 double pc,s,h,rhow;
01315
01316 s = f_s(rhov,t);
01317 h = f_h(rhov,t);
01318 rhow = f_rhow(t);
01319
01320 switch (model){
01321 case 1:{
01322 pc=0.0;
01323 break;
01324 }
01325 case 2:{
01326 if(s <= sssp){
01327
01328 pc=-1.0*rv*t*rhow*log(h);
01329 }
01330 else{
01331 pc=-1.0*rv*t*rhow*log(h);
01332 if(pc < 0.0) {
01333 pc=0.0;
01334 }
01335 }
01336 break;
01337 }
01338 default:{
01339 fprintf (stderr,"\n\n unknown model type in glasgowmat.cpp (%s, line %d).\n",__FILE__,__LINE__);
01340 }
01341 }
01342 return(pc);
01343 }
01344
01345
01346
01347
01348
01349 double glasgowmat::f_dpcdt (double rhov, double t)
01350 {
01351 double dpcdt,pc,s,h,rhow,drhowdt,psat,dpsatdt;
01352
01353 s = f_s(rhov,t);
01354 h = f_h(rhov,t);
01355 rhow = f_rhow(t);
01356 drhowdt = f_drhowdt(t);
01357 psat = f_psat(t);
01358 dpsatdt = f_dpsatdt(t);
01359
01360 switch (model){
01361 case 1:{
01362 dpcdt=0.0;
01363 break;
01364 }
01365 case 2:{
01366 if(s <= sssp){
01367 dpcdt=0.0;
01368 }
01369 else{
01370 pc=-1.0*rv*t*rhow*log(h);
01371 if(pc < 0.0) {
01372 dpcdt=0.0;
01373 }
01374 else{
01375 dpcdt=rv*rhow*((t*dpsatdt/psat)-((1.0+t*drhowdt/rhow)*log(rhov*rv*t/psat))-1.0);
01376 }
01377 }
01378 break;
01379 }
01380 default:{
01381 fprintf (stderr,"\n\n unknown model type in glasgowmat.cpp (%s, line %d).\n",__FILE__,__LINE__);
01382 }
01383 }
01384
01385 return(dpcdt);
01386 }
01387
01388
01389
01390
01391
01392 double glasgowmat::f_dpcdrhov (double rhov, double t)
01393 {
01394 double dpcdrhov,pc,s,h,rhow;
01395
01396 s = f_s(rhov,t);
01397 h = f_h(rhov,t);
01398 rhow = f_rhow(t);
01399
01400 switch (model){
01401 case 1:{
01402 dpcdrhov=0.0;
01403 break;
01404 }
01405 case 2:{
01406 if(s <= sssp){
01407 dpcdrhov=0.0;
01408 }
01409 else{
01410 pc=-1.0*rv*t*rhow*log(h);
01411 if(pc < 0.0) {
01412 dpcdrhov=0.0;
01413 }
01414 else{
01415 dpcdrhov=-rv*t*rhow/rhov;
01416 }
01417 }
01418 break;
01419 }
01420 default:{
01421 fprintf (stderr,"\n\n unknown model type in glasgowmat.cpp (%s, line %d).\n",__FILE__,__LINE__);
01422 }
01423 }
01424
01425 return(dpcdrhov);
01426 }
01427
01428
01429
01430
01431
01432 double glasgowmat::f_sln (double rhov, double t)
01433 {
01434 double sln,s;
01435
01436 s = f_s(rhov,t);
01437
01438 switch (model){
01439 case 1:{
01440 sln=0.0;
01441 break;
01442 }
01443 case 2:{
01444 if(s <= sssp){
01445 sln=0.0;
01446 }
01447 else{
01448 sln=(s-sssp)/(1.0-sssp);
01449 }
01450 break;
01451 }
01452 default:{
01453 fprintf (stderr,"\n\n unknown model type in glasgowmat.cpp (%s, line %d).\n",__FILE__,__LINE__);
01454 }
01455 }
01456
01457 return(sln);
01458 }
01459
01460
01461
01462
01463
01464 double glasgowmat::f_ppore (double rhov, double pg, double t, double pginf)
01465 {
01466 double ppore;
01467 double sln,pc,s,h,rhow;
01468
01469 s = f_s(rhov,t);
01470 h = f_h(rhov,t);
01471 rhow = f_rhow(t);
01472
01473 switch (model){
01474 case 1:{
01475 ppore=pg-pginf;
01476 pc=0.0;
01477 sln=0.0;
01478 break;
01479 }
01480 case 2:{
01481 if(s <= sssp){
01482 pc=0.0;
01483 sln=0.0;
01484 ppore=pg-pginf;
01485 }
01486 else{
01487 pc=-1.0*rv*t*rhow*log(h);
01488 if(pc < 0.0) {
01489 pc=0.0;
01490 }
01491 else{}
01492 sln=(s-sssp)/(1.0-sssp);
01493 ppore=pg-(sln*pc)-pginf;
01494 }
01495 break;
01496 }
01497 default:{
01498 fprintf (stderr,"\n\n unknown model type in glasgowmat.cpp (%s, line %d).\n",__FILE__,__LINE__);
01499 }
01500 }
01501
01502 return(ppore);
01503 }
01504
01505
01506
01507
01508
01509 double glasgowmat::f_sb (double rhov, double t)
01510 {
01511 double sb,s;
01512
01513 s = f_s(rhov,t);
01514
01515 switch (model){
01516 case 1:{
01517 sb=0.0;
01518 break;
01519 }
01520 case 2:{
01521 if(s <= sssp){
01522 sb=s;
01523 }
01524 else{
01525 sb=sssp;
01526 }
01527 break;
01528 }
01529 default:{
01530 fprintf (stderr,"\n\n unknown model type in glasgowmat.cpp (%s, line %d).\n",__FILE__,__LINE__);
01531 }
01532 }
01533 return(sb);
01534 }
01535
01536
01537
01538
01539
01540 double glasgowmat::f_pa (double rhov, double pg, double t)
01541 {
01542 double pa,pv;
01543
01544 pv = f_pv(rhov,t);
01545
01546 pa=pg-pv;
01547
01548 return(pa);
01549 }
01550
01551
01552
01553
01554
01555 double glasgowmat::f_rhoa (double rhov, double pg, double t)
01556 {
01557 double rhoa,pa;
01558
01559 pa = f_pa(rhov,pg,t);
01560
01561 rhoa=pa/ra/t;
01562
01563 return(rhoa);
01564 }
01565
01566
01567
01568
01569
01570 double glasgowmat::f_rhog (double rhov, double pg, double t)
01571 {
01572 double rhog,rhoa;
01573
01574 rhoa = f_rhoa(rhov,pg,t);
01575
01576 rhog=rhoa+rhov;
01577
01578
01579
01580 return(rhog);
01581 }
01582
01583
01584
01585
01586
01587 double glasgowmat::f_dav (double rhov, double pg, double t)
01588 {
01589 double dif,delta,tor,dav;
01590
01591
01592
01593
01594
01595
01596
01597 dif=1.87*exp(2.072*log(t))/pg*1.0e-5;
01598
01599 delta=0.5;
01600 tor=3.0;
01601 dav=dif*delta/tor/tor;
01602
01603
01604
01605
01606
01607 return(dav);
01608 }
01609
01610
01611
01612
01613
01614 double glasgowmat::f_davex (double pg, double tinf)
01615 {
01616 double davex,tor,delta,difex;
01617
01618
01619
01620
01621
01622
01623
01624 delta=0.5;
01625 tor=3.0;
01626
01627
01628
01629 difex=1.87*exp(2.072*log(tinf))/pg*1.0e-5;
01630
01631 davex=difex*delta/tor/tor;
01632
01633 return(davex);
01634 }
01635
01636
01637
01638
01639
01640 double glasgowmat::f_db(double rhov, double t)
01641 {
01642 double db,db0,dbtref,s;
01643
01644 s = f_s(rhov,t);
01645
01646
01647 db0=1.57e-11;
01648 dbtref=295.0;
01649 if(s <= sssp)
01650 db=db0*exp(-2.08*(s/sssp)*(t/dbtref));
01651 else
01652 db=0.0;
01653
01654 return(db);
01655 }
01656
01657
01658
01659
01660
01661
01662 double glasgowmat::f_kk (double t)
01663 {
01664 double kk,por;
01665
01666
01667
01668
01669
01670
01671
01672
01673
01674
01675
01676
01677
01678
01679 por = f_por(t);
01680
01681
01682 kk=exp((2.0/3.0)*log(por/por0))*k0;
01683
01684
01685 kk = k0;
01686
01687
01688 return(kk);
01689 }
01690
01691
01692
01693
01694
01695 double glasgowmat::f_kg (double rhov, double t)
01696 {
01697 double kg,kb,km,s;
01698
01699
01700
01701
01702
01703
01704
01705
01706
01707
01708
01709
01710
01711
01712
01713
01714
01715
01716
01717
01718
01719 s = f_s(rhov,t);
01720
01721
01722 kb=2.2748;
01723 km=1.0/kb;
01724
01725 kg=(sqrt(1.0-s))*exp((2.0*km)*log(1.0-exp((1.0/km)*log(s))));
01726
01727
01728
01729
01730
01731
01732 return(kg);
01733 }
01734
01735
01736
01737
01738
01739 double glasgowmat::f_kl (double rhov, double t)
01740 {
01741 double kl,km,kb,s;
01742
01743
01744
01745
01746
01747
01748
01749
01750
01751
01752
01753
01754
01755
01756
01757
01758
01759
01760
01761
01762
01763 s = f_s(rhov,t);
01764
01765
01766 kb=2.2748;
01767 km=1.0/kb;
01768
01769 kl=sqrt(s)*exp(2.0*log(1.0-exp(km*log(1.0-exp((1.0/km)*log(s))))));
01770
01771
01772
01773
01774
01775
01776 return(kl);
01777 }
01778
01779
01780
01781
01782
01783 double glasgowmat::f_mul (double t)
01784 {
01785 double mul;
01786
01787
01788
01789
01790
01791 if(t <= 647.3)
01792
01793
01794 mul=0.6612*exp((-1.562)*log(t-229.0));
01795 else
01796 mul=0.000043;
01797
01798 return(mul);
01799 }
01800
01801
01802
01803
01804
01805
01806 double glasgowmat::f_muv (double t)
01807 {
01808 double muv,t0muv,amuv,muv0;
01809
01810
01811
01812
01813
01814 if(t <= 647.3){
01815 t0muv=273.15;
01816 muv0=8.85e-6;
01817 amuv=3.53e-8;
01818 muv=muv0+amuv*(t-t0muv);
01819 }
01820 else
01821 muv=0.00004313;
01822
01823 return(muv);
01824 }
01825
01826
01827
01828
01829
01830
01831 double glasgowmat::f_mua (double t)
01832 {
01833 double mua,t0mua,amua,bmua,mua0;
01834
01835
01836
01837
01838
01839 t0mua=273.15;
01840 mua0=17.17e-6;
01841 amua=4.73e-8;
01842 bmua=2.22e-11;
01843
01844
01845 mua=mua0+(amua*(t-t0mua))-(bmua*exp(2.0*log(t-t0mua)));
01846
01847 return(mua);
01848 }
01849
01850
01851
01852
01853
01854
01855 double glasgowmat::f_mug (double rhov, double pg, double t)
01856 {
01857 double mug,mua,muv,rhoa;
01858
01859 mua = f_mua(t);
01860 muv = f_muv(t);
01861 rhoa = f_rhoa(rhov,pg,t);
01862
01863
01864 if((rhoa == 0.0) && (rhov == 0.0))
01865 mug=0.0;
01866 else
01867 mug=((rhoa*mua)+(rhov*muv))/(rhoa+rhov);
01868
01869
01870
01871
01872 return(mug);
01873 }
01874
01875
01876
01877
01878
01879
01880 double glasgowmat::f_cl (double t)
01881 {
01882 double cl;
01883
01884
01885
01886
01887
01888 if(t <= 647.3)
01889
01890
01891
01892 cl=((2.4768*t)+3368.2)+exp(31.4447657616636*log((1.08542631988638*t)/513.15));
01893 else
01894 cl=24515.0;
01895
01896 return(cl);
01897 }
01898
01899
01900
01901
01902
01903
01904 double glasgowmat::f_cv (double t)
01905 {
01906 double cv;
01907
01908
01909
01910
01911
01912 if(t <= 647.3)
01913
01914
01915
01916 cv=((7.1399*t)-443.0)+exp(29.4435287521143*log((1.13771502228162*t)/513.15));
01917 else
01918 cv=45821.04;
01919
01920 return(cv);
01921 }
01922
01923
01924
01925
01926
01927
01928 double glasgowmat::f_ca (double t)
01929 {
01930 double ca;
01931
01932
01933
01934
01935 ca = -0.00000009849367018147350*t*t*t + 0.0003564362577698610*t*t - 0.1216179239877570*t +1012.502552163240;
01936
01937 return(ca);
01938 }
01939
01940
01941
01942
01943
01944
01945 double glasgowmat::f_cs (double t)
01946 {
01947 double cs,tc;
01948
01949 tc=t-273.15;
01950
01951
01952
01953
01954
01955 cs=900.0+80.0*(tc/120.0)-(4.0*tc*tc/120.0/120.0);
01956
01957
01958
01959 return(cs);
01960 }
01961
01962
01963
01964
01965
01966
01967 double glasgowmat::f_keff (double t)
01968 {
01969 double keff,tc;
01970
01971 tc=t-273.15;
01972
01973
01974
01975
01976
01977 keff=2.0-(0.24*(tc/120.0))+(0.012*tc*tc/120.0/120.0);
01978
01979
01980
01981
01982 return(keff);
01983 }
01984
01985
01986
01987
01988
01989
01990 double glasgowmat::f_fracs (double rhov, double t)
01991 {
01992 double fracs,fracg,fracl;
01993
01994 fracg = f_fracg(rhov,t);
01995 fracl = f_fracl(rhov,t);
01996
01997 fracs=1.0-fracg-fracl;
01998
01999 return(fracs);
02000 }
02001
02002
02003
02004
02005
02006
02007 double glasgowmat::f_crho (double rhov, double pg, double t)
02008 {
02009 double crho,fracs,cs,fracl,rhow,cl,cv,fracg,rhoa,ca;
02010
02011 fracs = f_fracs(rhov,t);
02012 cs = f_cs(t);
02013 fracl = f_fracl(rhov,t);
02014 rhow = f_rhow(t);
02015 cl = f_cl(t);
02016 fracg = f_fracg(rhov,t);
02017 cv = f_cv(t);
02018 rhoa = f_rhoa(rhov,pg,t);
02019 ca = f_ca(t);
02020
02021 crho=(rhos*fracs*cs)+(fracl*rhow*cl)+(fracg*rhov*cv)+(fracg*rhoa*ca);
02022
02023 return(crho);
02024 }
02025
02026
02027
02028
02029
02030
02031 double glasgowmat::f_le (double t)
02032 {
02033 double le;
02034
02035
02036
02037
02038
02039
02040
02041
02042
02043 if(t <= 647.3)
02044
02045
02046 le=2.672e5*exp(0.38*log(647.3-t));
02047 else
02048 le=0.0;
02049
02050 return(le);
02051 }
02052
02053
02054
02055
02056
02057
02058 double glasgowmat::f_ld ()
02059 {
02060 double ld;
02061
02062
02063 ld=2400.0e3;
02064
02065 return(ld);
02066 }
02067
02068
02069
02070
02071
02072
02073 double glasgowmat::f_hrad(double t, double tinf)
02074 {
02075 double hrad;
02076
02077
02078
02079
02080
02081 hrad=emmi*stef*((t*t)+(tinf*tinf))*(t+tinf);
02082
02083 return(hrad);
02084 }
02085
02086
02087
02088
02089
02090 double glasgowmat::f_hqr(double t, double tinf)
02091 {
02092 double hrad,hqr;
02093
02094 hrad=f_hrad(t,tinf);
02095
02096
02097
02098
02099 hqr=hq+hrad;
02100
02101 return(hqr);
02102 }
02103
02104
02105
02106
02107
02108
02109 double glasgowmat::f_beta (double pg, double tinf)
02110 {
02111 double beta,davex;
02112
02113 davex = f_davex(pg,tinf);
02114
02115
02116
02117
02118 beta=(hq/crhoair)*exp((2.0/3.0)*log(davex/alph));
02119
02120 return(beta);
02121 }
02122
02123
02124
02125
02126
02127
02128 double glasgowmat::ktt (double t,double pg,double rhov)
02129 {
02130 double ktt,keff;
02131
02132 keff = f_keff(t);
02133
02134 switch (model){
02135 case 1:{
02136 ktt = keff;
02137 break;
02138 }
02139 case 2:{
02140 double klt,fracl,rhow,sb,db,s,por,dfracldt,dpordt,kk,kl,dpcdt,mul,le;
02141
02142 fracl = f_fracl(rhov,t);
02143 rhow = f_rhow(t);
02144 sb = f_sb(rhov,t);
02145 db = f_db(rhov,t);
02146 s = f_s(rhov,t);
02147 por = f_por(t);
02148 dfracldt = f_dfracldt(rhov,t);
02149 dpordt = f_dpordt(t);
02150 kk = f_kk(t);
02151 kl = f_kl(rhov,t);
02152 dpcdt = f_dpcdt(rhov,t);
02153 mul = f_mul(t);
02154 le = f_le(t);
02155
02156 klt=fracl*rhow*(((sb*db/s/por)*(dfracldt-(fracl*dpordt/por)))-(1.0-(sb/s))*(kk*kl*dpcdt/mul));
02157 ktt = keff-(le*klt);
02158
02159
02160
02161
02162
02163
02164
02165
02166
02167 break;
02168 }
02169 default:{
02170 fprintf (stderr,"\n\n unknown model type in glasgowmat.cpp (%s, line %d).\n",__FILE__,__LINE__);
02171 }
02172 }
02173
02174 return (ktt);
02175 }
02176
02177 double glasgowmat::ktp (double t,double pg,double rhov)
02178 {
02179 double ktp,klp;
02180 double sb,s,fracl,rhow,kk,kl,mul,le;
02181
02182 sb = f_sb(rhov,t);
02183 s = f_s(rhov,t);
02184 fracl = f_fracl(rhov,t);
02185 rhow = f_rhow(t);
02186 kk = f_kk(t);
02187 kl = f_kl(rhov,t);
02188 mul = f_mul(t);
02189 le = f_le(t);
02190
02191 switch (model){
02192 case 1:{
02193 ktp= -le*fracl*rhow*kk*kl/mul;
02194 break;
02195 }
02196 case 2:{
02197 klp=(1.0-(sb/s))*fracl*rhow*kk*kl/mul;
02198 ktp= -le*klp;
02199 break;
02200 }
02201 default:{
02202 fprintf (stderr,"\n\n unknown model type in glasgowmat.cpp (%s, line %d).\n",__FILE__,__LINE__);
02203 }
02204 }
02205
02206 return (ktp);
02207 }
02208
02209 double glasgowmat::ktv (double t,double pg,double rhov)
02210 {
02211 double ktv;
02212 double klv,db,sb,s,fracl,rhow,kk,kl,mul,le,dfracldrhov,dpcdrhov,por;
02213
02214 sb = f_sb(rhov,t);
02215 db = f_db(rhov,t);
02216 s = f_s(rhov,t);
02217 fracl = f_fracl(rhov,t);
02218 dfracldrhov = f_dfracldrhov(rhov,t);
02219 dpcdrhov = f_dpcdrhov(rhov,t);
02220 rhow = f_rhow(t);
02221 kk = f_kk(t);
02222 por = f_por(t);
02223 kl = f_kl(rhov,t);
02224 mul = f_mul(t);
02225 le = f_le(t);
02226
02227 switch (model){
02228 case 1:{
02229 ktv=0.0;
02230 break;
02231 }
02232 case 2:{
02233 klv=fracl*rhow*((sb*db*dfracldrhov/s/por)-((1.0-(sb/s))*kk*kl*dpcdrhov/mul));
02234 ktv= -le*klv;
02235 break;
02236 }
02237 default:{
02238 fprintf (stderr,"\n\n unknown model type in glasgowmat.cpp (%s, line %d).\n",__FILE__,__LINE__);
02239 }
02240 }
02241
02242 return (ktv);
02243 }
02244
02245 double glasgowmat::kat (double t,double pg,double rhov)
02246 {
02247 double kat,dav,fracg,rhog;
02248
02249 dav = f_dav(rhov,pg,t);
02250 fracg = f_fracg(rhov,t);
02251 rhog = f_rhog(rhov,pg,t);
02252
02253
02254
02255
02256
02257
02258
02259
02260
02261 kat= -1.0*dav*fracg*rhov*pg/rhog/ra/t/t;
02262
02263 return (kat);
02264 }
02265
02266 double glasgowmat::kap (double t,double pg,double rhov)
02267 {
02268 double kap;
02269 double dav,fracg,rhoa,rhog,kk,kg,mug;
02270
02271 dav = f_dav(rhov,pg,t);
02272 fracg = f_fracg(rhov,t);
02273 rhoa = f_rhoa(rhov,pg,t);
02274 rhog = f_rhog(rhov,pg,t);
02275 kk = f_kk(t);
02276 kg = f_kg(rhov,t);
02277 mug = f_mug(rhov,pg,t);
02278
02279 kap= (kk*kg*fracg*rhoa/mug)+(dav*fracg*rhov/rhog/ra/t);
02280
02281 return (kap);
02282 }
02283
02284 double glasgowmat::kav (double t,double pg,double rhov)
02285 {
02286 double kav;
02287 double dav,fracg,rhoa,rhog;
02288
02289 dav = f_dav(rhov,pg,t);
02290 fracg = f_fracg(rhov,t);
02291 rhoa = f_rhoa(rhov,pg,t);
02292 rhog = f_rhog(rhov,pg,t);
02293
02294 kav= -1.0*(rhoa+(rv*rhov/ra))*dav*fracg/rhog;
02295
02296 return (kav);
02297 }
02298
02299 double glasgowmat::kmt (double t,double pg,double rhov)
02300 {
02301 double kmt,klt,kvt;
02302 double dav,fracg,rhoa,rhog;
02303
02304 dav = f_dav(rhov,pg,t);
02305 fracg = f_fracg(rhov,t);
02306 rhoa = f_rhoa(rhov,pg,t);
02307 rhog = f_rhog(rhov,pg,t);
02308
02309 switch (model){
02310 case 1:{
02311 kvt=fracg*dav*rhov*pg/rhog/ra/t/t;
02312 klt=0.0;
02313 kmt=kvt+klt;
02314 break;
02315 }
02316 case 2:{
02317 double db,sb,s,fracl,rhow,kk,kl,mul,le,dfracldt,dpcdt,por,dpordt;
02318
02319 sb = f_sb(rhov,t);
02320 db = f_db(rhov,t);
02321 s = f_s(rhov,t);
02322 fracl = f_fracl(rhov,t);
02323 dfracldt = f_dfracldt(rhov,t);
02324 dpcdt = f_dpcdt(rhov,t);
02325 rhow = f_rhow(t);
02326 kk = f_kk(t);
02327 por = f_por(t);
02328 dpordt = f_dpordt(t);
02329 kl = f_kl(rhov,t);
02330 mul = f_mul(t);
02331 le = f_le(t);
02332
02333 kvt=fracg*dav*rhov*pg/rhog/ra/t/t;
02334 klt=fracl*rhow*(((sb*db/s/por)*(dfracldt-(fracl*dpordt/por)))-(1.0-(sb/s))*(kk*kl*dpcdt/mul));
02335 kmt=kvt+klt;
02336 break;
02337 }
02338 default:{
02339 fprintf (stderr,"\n\n unknown model type in glasgowmat.cpp (%s, line %d).\n",__FILE__,__LINE__);
02340 }
02341 }
02342
02343 return (kmt);
02344 }
02345
02346 double glasgowmat::kmp (double t,double pg,double rhov)
02347 {
02348 double kmp,kvp,klp;
02349 double dav,fracg,rhoa,rhog;
02350 double sb,s,fracl,rhow,kk,kl,mul,le;
02351 double kg,mug;
02352
02353 dav = f_dav(rhov,pg,t);
02354 fracg = f_fracg(rhov,t);
02355 rhoa = f_rhoa(rhov,pg,t);
02356 rhog = f_rhog(rhov,pg,t);
02357 sb = f_sb(rhov,t);
02358 s = f_s(rhov,t);
02359 fracl = f_fracl(rhov,t);
02360 rhow = f_rhow(t);
02361 kk = f_kk(t);
02362 kl = f_kl(rhov,t);
02363 kg = f_kg(rhov,t);
02364 mul = f_mul(t);
02365 mug = f_mug(rhov,pg,t);
02366 le = f_le(t);
02367
02368 switch (model){
02369 case 1:{
02370 kvp=fracg*rhov*((kk*kg/mug)-(dav/rhog/ra/t));
02371 klp=fracl*rhow*(kk*kl/mul);
02372 kmp=kvp+klp;
02373 break;
02374 }
02375 case 2:{
02376 kvp=fracg*rhov*((kk*kg/mug)-(dav/rhog/ra/t));
02377 klp=(1.0-(sb/s))*fracl*rhow*kk*kl/mul;
02378 kmp=kvp+klp;
02379 break;
02380 }
02381 default:{
02382 fprintf (stderr,"\n\n unknown model type in glasgowmat.cpp (%s, line %d).\n",__FILE__,__LINE__);
02383 }
02384 }
02385 return(kmp);
02386 }
02387
02388 double glasgowmat::kmv (double t,double pg,double rhov)
02389 {
02390 double kmv,kvv,klv;
02391 double dav,fracg,fracl,rhoa,rhog,sb,db,s,dfracldrhov,rhow,kk,kl,mul,por,dpcdrhov;
02392
02393 dav = f_dav(rhov,pg,t);
02394 fracg = f_fracg(rhov,t);
02395 rhoa = f_rhoa(rhov,pg,t);
02396 rhog = f_rhog(rhov,pg,t);
02397 sb = f_sb(rhov,t);
02398 db = f_db(rhov,t);
02399 s = f_s(rhov,t);
02400 fracl = f_fracl(rhov,t);
02401 dfracldrhov = f_dfracldrhov(rhov,t);
02402 rhow = f_rhow(t);
02403 kk = f_kk(t);
02404 kl = f_kl(rhov,t);
02405 mul = f_mul(t);
02406 por = f_por(t);
02407 dpcdrhov = f_dpcdrhov(rhov,t);
02408
02409 switch (model){
02410 case 1:{
02411 kvv=(fracg*dav/rhog)*(rhoa+(rv*rhov/ra));
02412 klv=0.0;
02413 kmv=kvv+klv;
02414 break;
02415 }
02416 case 2:{
02417 kvv=(fracg*dav/rhog)*(rhoa+(rv*rhov/ra));
02418 klv=fracl*rhow*((sb*db*dfracldrhov/s/por)-((1.0-(sb/s))*kk*kl*dpcdrhov/mul));
02419 kmv=kvv+klv;
02420 break;
02421 }
02422 default:{
02423 fprintf (stderr,"\n\n unknown model type in glasgowmat.cpp (%s, line %d).\n",__FILE__,__LINE__);
02424 }
02425 }
02426 return (kmv);
02427 }
02428
02429 double glasgowmat::ctt (double t,double pg,double rhov)
02430 {
02431 double ctt;
02432 double crho,le,ld,rhow,drhowdt,fracl,fracd,dfracldt,dfracddt;
02433
02434 crho = f_crho(rhov,pg,t);
02435 rhow = f_rhow(t);
02436 drhowdt = f_drhowdt(t);
02437 fracl = f_fracl(rhov,t);
02438 fracd = f_fracd(t);
02439 dfracldt = f_dfracldt(rhov,t);
02440 dfracddt = f_dfracddt(t);
02441 ld = f_ld();
02442 le = f_le(t);
02443
02444 ctt = (crho+((ld+le)*(fracd*drhowdt+rhow*dfracddt))-(le*(fracl*drhowdt+rhow*dfracldt)));
02445
02446 return(ctt);
02447 }
02448
02449 double glasgowmat::ctp ()
02450 {
02451 return 0.0;
02452 }
02453
02454 double glasgowmat::ctv (double t,double pg,double rhov)
02455 {
02456 double ctv;
02457 double le,rhow,dfracldrhov;
02458
02459 le = f_le(t);
02460 rhow = f_rhow(t);
02461 dfracldrhov = f_dfracldrhov(rhov,t);
02462
02463 ctv = (-le*rhow*dfracldrhov);
02464
02465 return(ctv);
02466 }
02467
02468 double glasgowmat::cat (double t,double pg,double rhov)
02469 {
02470 double cat;
02471 double dpordt,rhoa,fracg,dfracldt;
02472
02473 dpordt = f_dpordt(t);
02474 rhoa = f_rhoa(rhov,pg,t);
02475 fracg = f_fracg(rhov,t);
02476 dfracldt = f_dfracldt(rhov,t);
02477
02478 cat = ((rhoa*(dpordt-dfracldt))-(fracg*pg/ra/t/t));
02479
02480 return(cat);
02481 }
02482
02483 double glasgowmat::cap (double t,double pg,double rhov)
02484 {
02485 double cap;
02486 double fracg;
02487
02488 fracg = f_fracg(rhov,t);
02489
02490 cap = (fracg/ra/t);
02491
02492 return(cap);
02493 }
02494
02495 double glasgowmat::cav (double t,double pg,double rhov)
02496 {
02497 double cav;
02498 double fracg,rhoa,dfracldrhov;
02499
02500 fracg = f_fracg(rhov,t);
02501 rhoa = f_rhoa(rhov,pg,t);
02502 dfracldrhov = f_dfracldrhov(rhov,t);
02503
02504 cav = (-1.0*(fracg*rv/ra)-(rhoa*dfracldrhov));
02505
02506 return(cav);
02507 }
02508
02509 double glasgowmat::cmt (double t,double pg,double rhov)
02510 {
02511 double cmt;
02512 double dpordt,rhow,drhowdt,fracl,fracd,dfracldt,dfracddt;
02513
02514 dpordt = f_dpordt(t);
02515 rhow = f_rhow(t);
02516 drhowdt = f_drhowdt(t);
02517 fracl = f_fracl(rhov,t);
02518 fracd = f_fracd(t);
02519 dfracldt = f_dfracldt(rhov,t);
02520 dfracddt = f_dfracddt(t);
02521
02522 cmt = ((rhov*dpordt)+((fracl-fracd)*drhowdt)+((rhow-rhov)*dfracldt)-(rhow*dfracddt));
02523
02524 return(cmt);
02525 }
02526
02527 double glasgowmat::cmp ()
02528 {
02529 return 0.0;
02530 }
02531
02532 double glasgowmat::cmv (double t,double pg,double rhov)
02533 {
02534 double cmv;
02535 double fracg,rhow,dfracldrhov;
02536
02537 fracg = f_fracg(rhov,t);
02538 rhow = f_rhow(t);
02539 dfracldrhov = f_dfracldrhov(rhov,t);
02540
02541 cmv = (fracg+((rhow-rhov)*dfracldrhov));
02542
02543 return(cmv);
02544 }
02545
02546
02547
02548
02549
02550
02551
02552
02553
02554
02555
02556
02557
02558
02559
02560 double glasgowmat::transmission_transcoeff(double trc,long ri,long ci,long nn,long bc,long ipp)
02561 {
02562 long k;
02563 double new_trc,rhov,pg,t;
02564 new_trc = 0.0;
02565
02566 k=Gtt->give_dof(nn,0);
02567 if (k>0) {rhov = Lsrst->lhs[k-1]+Lsrst->lhsi[k-1];}
02568 if (k==0) {rhov = 0.0;}
02569 if (k<0) {rhov = Tb->lc[0].pv[0-k-1].getval();}
02570 k=Gtt->give_dof(nn,1);
02571 if (k>0) {pg = Lsrst->lhs[k-1]+Lsrst->lhsi[k-1];}
02572 if (k==0) {pg = 0.0;}
02573 if (k<0) {pg = Tb->lc[0].pv[0-k-1].getval();}
02574 k=Gtt->give_dof(nn,2);
02575 if (k>0) {t = Lsrst->lhs[k-1]+Lsrst->lhsi[k-1];}
02576 if (k==0) {t = 0.0;}
02577 if (k<0) {t = Tb->lc[0].pv[0-k-1].getval();}
02578
02579 if((ri == 0) && (ci == 0))
02580 new_trc = get_transmission_transcoeff_mv(t,pg,rhov,bc,ipp);
02581 if((ri == 0) && (ci == 1))
02582 new_trc = 0.0;
02583 if((ri == 0) && (ci == 2))
02584 new_trc = get_transmission_transcoeff_mt(t,pg,rhov,bc,ipp);
02585
02586 if((ri == 1) && (ci == 0))
02587 new_trc = 0.0;
02588 if((ri == 1) && (ci == 1))
02589 new_trc = 0.0;
02590 if((ri == 1) && (ci == 2))
02591 new_trc = 0.0;
02592
02593 if((ri == 2) && (ci == 0))
02594 new_trc = 0.0;
02595 if((ri == 2) && (ci == 1))
02596 new_trc = 0.0;
02597 if((ri == 2) && (ci == 2))
02598 new_trc = get_transmission_transcoeff_tt(t,pg,rhov,bc,ipp);
02599
02600 new_trc = new_trc*trc;
02601
02602 return (new_trc);
02603 }
02604
02605
02606
02607
02608
02609
02610
02611
02612
02613
02614
02615
02616
02617
02618 double glasgowmat::transmission_nodval(double nodval,double trc2,long ri,long ci,long nn,long bc,long ipp)
02619 {
02620 long k;
02621 double new_nodval,rhov,pg,t;
02622 new_nodval = 0.0;
02623
02624 k=Gtt->give_dof(nn,0);
02625 if (k>0) {rhov = Lsrst->lhs[k-1]+Lsrst->lhsi[k-1];}
02626 if (k==0) {rhov = 0.0;}
02627 if (k<0) {rhov = Tb->lc[0].pv[0-k-1].getval();}
02628 k=Gtt->give_dof(nn,1);
02629 if (k>0) {pg = Lsrst->lhs[k-1]+Lsrst->lhsi[k-1];}
02630 if (k==0) {pg = 0.0;}
02631 if (k<0) {pg = Tb->lc[0].pv[0-k-1].getval();}
02632 k=Gtt->give_dof(nn,2);
02633 if (k>0) {t = Lsrst->lhs[k-1]+Lsrst->lhsi[k-1];}
02634 if (k==0) {t = 0.0;}
02635 if (k<0) {t = Tb->lc[0].pv[0-k-1].getval();}
02636
02637 if((ri == 0) && (ci == 0))
02638 new_nodval = get_transmission_nodval_mv(nodval,t,pg,rhov,bc,ipp);
02639 if((ri == 0) && (ci == 1))
02640 new_nodval = 0.0;
02641 if((ri == 0) && (ci == 2))
02642 new_nodval = get_transmission_nodval_mt(nodval,t,pg,rhov,bc,ipp);
02643
02644 if((ri == 1) && (ci == 0))
02645 new_nodval = 0.0;
02646 if((ri == 1) && (ci == 1))
02647 new_nodval = 0.0;
02648 if((ri == 1) && (ci == 2))
02649 new_nodval = 0.0;
02650
02651 if((ri == 2) && (ci == 0))
02652 new_nodval = 0.0;
02653 if((ri == 2) && (ci == 1))
02654 new_nodval = 0.0;
02655 if((ri == 2) && (ci == 2))
02656 new_nodval = get_transmission_nodval_tt(nodval,t,pg,rhov,bc,ipp);
02657
02658 return (new_nodval);
02659 }
02660
02661
02662
02663
02664
02665
02666
02667
02668
02669
02670
02671
02672
02673
02674
02675 double glasgowmat::transmission_flux(double nodval,double trc2,long ri,long ci,long nn,long bc,long ipp)
02676 {
02677 long k;
02678 double flux,rhov,pg,t;
02679 flux = 0.0;
02680
02681 k=Gtt->give_dof(nn,0);
02682 if (k>0) {rhov = Lsrst->lhs[k-1]+Lsrst->lhsi[k-1];}
02683 if (k==0) {rhov = 0.0;}
02684 if (k<0) {rhov = Tb->lc[0].pv[0-k-1].getval();}
02685 k=Gtt->give_dof(nn,1);
02686 if (k>0) {pg = Lsrst->lhs[k-1]+Lsrst->lhsi[k-1];}
02687 if (k==0) {pg = 0.0;}
02688 if (k<0) {pg = Tb->lc[0].pv[0-k-1].getval();}
02689 k=Gtt->give_dof(nn,2);
02690 if (k>0) {t = Lsrst->lhs[k-1]+Lsrst->lhsi[k-1];}
02691 if (k==0) {t = 0.0;}
02692 if (k<0) {t = Tb->lc[0].pv[0-k-1].getval();}
02693
02694 if((ri == 0) && (ci == 0))
02695 flux = get_transmission_flux_mv(nodval,t,pg,rhov,bc,ipp);
02696 if((ri == 0) && (ci == 1))
02697 flux = 0.0;
02698 if((ri == 0) && (ci == 2))
02699 flux = get_transmission_flux_mt(nodval,t,pg,rhov,bc,ipp);
02700
02701 if((ri == 1) && (ci == 0))
02702 flux = 0.0;
02703 if((ri == 1) && (ci == 1))
02704 flux = 0.0;
02705 if((ri == 1) && (ci == 2))
02706 flux = 0.0;
02707
02708 if((ri == 2) && (ci == 0))
02709 flux = 0.0;
02710 if((ri == 2) && (ci == 1))
02711 flux = 0.0;
02712 if((ri == 2) && (ci == 2))
02713 flux = get_transmission_flux_tt(nodval,t,pg,rhov,bc,ipp);
02714
02715 return (flux);
02716 }
02717
02718
02719
02720
02721
02722
02723
02724
02725
02726
02727
02728 double glasgowmat::get_transmission_transcoeff_mv(double t,double pg,double rhov,long bc,long ipp)
02729 {
02730 double trc;
02731
02732 switch (bc){
02733 case 30:{
02734 trc=1.0;
02735 break;
02736 }
02737 case 31:{
02738 double beta;
02739
02740 beta = f_beta(pg,t);
02741
02742 trc = beta;
02743 break;
02744 }
02745 case 32:{
02746 double beta,tinf;
02747
02748 tinf = f_tinf(Tp->time);
02749 beta = f_beta(pg,tinf);
02750
02751 trc=beta;
02752 break;
02753 }
02754 case 36:{
02755 trc=0.0;
02756 break;
02757 }
02758 default:{
02759 fprintf (stderr,"\n\n No real boundary condition is prescribed (%s, line %d).\n",__FILE__,__LINE__);
02760 exit(0);
02761 }
02762 }
02763
02764 return(trc);
02765 }
02766
02767
02768
02769
02770
02771
02772
02773
02774
02775
02776
02777
02778 double glasgowmat::get_transmission_nodval_mv(double bv,double t,double pg,double rhov,long bc,long ipp)
02779 {
02780 double nodval;
02781
02782 switch (bc){
02783 case 30:{
02784 nodval = bv;
02785 break;
02786 }
02787 case 31:{
02788 nodval = bv;
02789 break;
02790 }
02791 case 32:{
02792 double beta,tinf;
02793
02794 tinf = f_tinf(Tp->time);
02795 beta = f_beta(pg,tinf);
02796
02797 nodval = beta*bv;
02798 break;
02799 }
02800 case 36:{
02801 double psat;
02802
02803 psat = f_psat(t);
02804
02805 nodval = bv*psat/rv/t;
02806
02807 nodval = nodval - rhov;
02808 break;
02809 }
02810 default:{
02811 fprintf (stderr,"\n\n No real boundary condition is prescribed (%s, line %d).\n",__FILE__,__LINE__);
02812 exit(0);
02813 }
02814 }
02815 return(nodval);
02816 }
02817
02818
02819
02820
02821
02822
02823
02824
02825
02826
02827
02828
02829 double glasgowmat::get_transmission_flux_mv(double bv,double t,double pg,double rhov,long bc,long ipp)
02830 {
02831 double flux;
02832
02833 switch (bc){
02834 case 30:{
02835 flux = (bv - rhov);
02836 break;
02837 }
02838 case 31:{
02839 double beta;
02840
02841 beta = f_beta(pg,t);
02842
02843 flux = beta*(bv - rhov);
02844 break;
02845 }
02846 case 32:{
02847 double beta,tinf;
02848
02849 tinf = f_tinf(Tp->time);
02850 beta = f_beta(pg,tinf);
02851
02852 flux = beta*(bv - rhov);
02853 break;
02854 }
02855 case 36:{
02856 double psat;
02857
02858 psat = f_psat(t);
02859
02860 flux = bv*psat/rv/t;
02861
02862 flux = flux - rhov;
02863 break;
02864 }
02865 default:{
02866 fprintf (stderr,"\n\n No real boundary condition is prescribed (%s, line %d).\n",__FILE__,__LINE__);
02867 exit(0);
02868 }
02869 }
02870 return(flux);
02871 }
02872
02873
02874
02875
02876
02877
02878
02879
02880
02881
02882 double glasgowmat::get_transmission_transcoeff_tt(double t,double pg,double rhov,long bc,long ipp)
02883 {
02884 double trc;
02885
02886 switch (bc){
02887 case 30:{
02888 trc=1.0;
02889 break;
02890 }
02891 case 31:{
02892 double hqr;
02893
02894 hqr = f_hqr(t,t);
02895
02896 trc=hqr;
02897 break;
02898 }
02899 case 32:{
02900 double hqr,tinf;
02901
02902 tinf = f_tinf(Tp->time);
02903 hqr = f_hqr(t,tinf);
02904
02905 trc=hqr;
02906 break;
02907 }
02908 default:{
02909 fprintf (stderr,"\n\n No real boundary condition is prescribed (%s, line %d).\n",__FILE__,__LINE__);
02910 exit(0);
02911 }
02912 }
02913 return(trc);
02914 }
02915
02916
02917
02918
02919
02920
02921
02922
02923
02924
02925
02926
02927
02928 double glasgowmat::get_transmission_nodval_tt(double bv,double t,double pg,double rhov,long bc,long ipp)
02929 {
02930 double nodval;
02931
02932 switch (bc){
02933 case 30:{
02934 nodval = bv;
02935 break;
02936 }
02937 case 31:{
02938 nodval = bv;
02939 break;
02940 }
02941 case 32:{
02942 double hqr,tinf;
02943
02944 tinf = f_tinf(Tp->time);
02945 hqr = f_hqr(t,tinf);
02946
02947 nodval = hqr*tinf;
02948 break;
02949 }
02950 default:{
02951 fprintf (stderr,"\n\n No real boundary condition is prescribed (%s, line %d).\n",__FILE__,__LINE__);
02952 exit(0);
02953 }
02954 }
02955 return(nodval);
02956 }
02957
02958
02959
02960
02961
02962
02963
02964
02965
02966
02967
02968
02969
02970 double glasgowmat::get_transmission_flux_tt(double bv,double t,double pg,double rhov,long bc,long ipp)
02971 {
02972 double flux;
02973
02974 switch (bc){
02975 case 30:{
02976 flux = (bv - t);
02977 break;
02978 }
02979 case 31:{
02980 double hqr;
02981
02982 hqr = f_hqr(t,bv);
02983
02984 flux = hqr*(bv - t);
02985 break;
02986 }
02987 case 32:{
02988 double hqr,tinf;
02989
02990 tinf = f_tinf(Tp->time);
02991 hqr = f_hqr(t,tinf);
02992
02993 flux = hqr*(tinf - t);
02994 break;
02995 }
02996 default:{
02997 fprintf (stderr,"\n\n No real boundary condition is prescribed (%s, line %d).\n",__FILE__,__LINE__);
02998 exit(0);
02999 }
03000 }
03001 return(flux);
03002 }
03003
03004
03005
03006
03007
03008
03009
03010
03011
03012
03013
03014 double glasgowmat::get_transmission_transcoeff_mt(double t,double pg,double rhov,long bc,long ipp)
03015 {
03016 double trc;
03017 double kvt,keff,hqr,tinf;
03018 double dav,fracg,rhoa,rhog;
03019
03020 dav = f_dav(rhov,pg,t);
03021 fracg = f_fracg(rhov,t);
03022 rhoa = f_rhoa(rhov,pg,t);
03023 rhog = f_rhog(rhov,pg,t);
03024
03025 kvt=fracg*dav*rhov*pg/rhog/ra/t/t;
03026 keff = f_keff (t);
03027
03028 switch (bc){
03029 case 30:{
03030
03031 trc=-1.0*kvt/keff;
03032 break;
03033 }
03034 case 31:{
03035 hqr = f_hqr(t,t);
03036
03037 trc=-1.0*kvt*hqr/keff;
03038 break;
03039 }
03040 case 32:{
03041 tinf = f_tinf(Tp->time);
03042 hqr = f_hqr(t,tinf);
03043
03044 trc=-1.0*kvt*hqr/keff;
03045 break;
03046 }
03047 case 36:{
03048 trc=0.0;
03049 break;
03050 }
03051 default:{
03052 fprintf (stderr,"\n\n No real boundary condition is prescribed (%s, line %d).\n",__FILE__,__LINE__);
03053 exit(0);
03054 }
03055 }
03056
03057 return(trc);
03058 }
03059
03060
03061
03062
03063
03064
03065
03066
03067
03068
03069
03070
03071 double glasgowmat::get_transmission_nodval_mt(double bv,double t,double pg,double rhov,long bc,long ipp)
03072 {
03073 double nodval;
03074 double kvt,keff,hqr,tinf;
03075 double dav,fracg,rhoa,rhog;
03076
03077 dav = f_dav(rhov,pg,t);
03078 fracg = f_fracg(rhov,t);
03079 rhoa = f_rhoa(rhov,pg,t);
03080 rhog = f_rhog(rhov,pg,t);
03081
03082 kvt=fracg*dav*rhov*pg/rhog/ra/t/t;
03083
03084 keff = f_keff (t);
03085
03086 switch (bc){
03087 case 30:{
03088
03089 nodval = -1.0*kvt*bv/keff;
03090 break;
03091 }
03092 case 31:{
03093 nodval = -1.0*kvt*bv/keff;
03094 break;
03095 }
03096 case 32:{
03097 tinf = f_tinf(Tp->time);
03098 hqr = f_hqr(t,tinf);
03099
03100 nodval = -1.0*kvt*hqr*tinf/keff;
03101 break;
03102 }
03103 case 36:{
03104 nodval=0.0;
03105 break;
03106 }
03107 default:{
03108 fprintf (stderr,"\n\n No real boundary condition is prescribed (%s, line %d).\n",__FILE__,__LINE__);
03109 exit(0);
03110 }
03111 }
03112 return(nodval);
03113 }
03114
03115
03116
03117
03118
03119
03120
03121
03122
03123
03124
03125 double glasgowmat::get_transmission_flux_mt(double bv,double t,double pg,double rhov,long bc,long ipp)
03126 {
03127 double flux;
03128 double kvt,keff,hqr,tinf;
03129 double dav,fracg,rhoa,rhog;
03130
03131 dav = f_dav(rhov,pg,t);
03132 fracg = f_fracg(rhov,t);
03133 rhoa = f_rhoa(rhov,pg,t);
03134 rhog = f_rhog(rhov,pg,t);
03135
03136 kvt=fracg*dav*rhov*pg/rhog/ra/t/t;
03137
03138 keff = f_keff (t);
03139
03140 switch (bc){
03141 case 30:{
03142
03143 flux = -1.0*kvt*(bv - t)/keff;
03144 break;
03145 }
03146 case 31:{
03147 hqr = f_hqr(t,bv);
03148
03149 flux = -1.0*kvt*hqr*(bv - t)/keff;
03150 break;
03151 }
03152 case 32:{
03153 tinf = f_tinf(Tp->time);
03154 hqr = f_hqr(t,tinf);
03155
03156 flux = -1.0*kvt*hqr*(tinf - t)/keff;
03157 break;
03158 }
03159 case 36:{
03160 flux=0.0;
03161 break;
03162 }
03163 default:{
03164 fprintf (stderr,"\n\n No real boundary condition is prescribed (%s, line %d).\n",__FILE__,__LINE__);
03165 exit(0);
03166 }
03167 }
03168 return(flux);
03169 }
03170
03171
03172
03173
03174
03175
03176
03177
03178
03179
03180
03181 double glasgowmat::get_othervalue(long compother,long ipp, double *r)
03182 {
03183 double other;
03184
03185 switch (compother){
03186 case 0:{
03187 other = r[0];
03188 break;
03189 }
03190 case 1:{
03191 other = r[1];
03192 break;
03193 }
03194 case 2:{
03195 other = r[2];
03196 break;
03197 }
03198 case 3:{
03199 double h;
03200 h = f_h(r[0],r[2]);
03201
03202 other = h;
03203 break;
03204 }
03205 case 4:{
03206 double s;
03207 s = f_s(r[0],r[2]);
03208
03209 other = s;
03210 break;
03211 }
03212 case 5:{
03213 double pv;
03214 pv = f_pv(r[0],r[2]);
03215
03216 other = pv;
03217 break;
03218 }
03219 case 6:{
03220 double pl;
03221 pl = f_pl(r[0],r[1],r[2]);
03222
03223 other = pl;
03224 break;
03225 }
03226 case 7:{
03227 double pc;
03228 pc = f_pc(r[0],r[2]);
03229
03230 other = pc;
03231 break;
03232 }
03233 case 8:{
03234 double pc;
03235 pc = f_pc(r[0],r[2]);
03236
03237 other = pc;
03238 break;
03239 }
03240 default:{
03241 fprintf (stderr,"\n\n unknown type of component is required in function (%s, line %d).\n",__FILE__,__LINE__);
03242 }
03243 }
03244 return (other);
03245 }
03246
03247
03248
03249
03250
03251
03252 void glasgowmat::print_othervalue_name(FILE *out,long compother)
03253 {
03254 switch (compother){
03255 case 0:{
03256 fprintf (out,"Vapour content (kg/m3) ");
03257 break;
03258 }
03259 case 1:{
03260 fprintf (out,"Gas pressure (Pa) ");
03261 break;
03262 }
03263 case 2:{
03264 fprintf (out,"Temperature (K) ");
03265 break;
03266 }
03267 case 3:{
03268 fprintf (out,"Relative humidity () ");
03269 break;
03270 }
03271 case 4:{
03272 fprintf (out,"Degree of saturation () ");
03273 break;
03274 }
03275 case 5:{
03276 fprintf (out,"Pore water vapor pressure (Pa)");
03277 break;
03278 }
03279 case 6:{
03280 fprintf (out,"Pore water pressure (Pa) ");
03281 break;
03282 }
03283 case 7:{
03284 fprintf (out,"Capilary pressure (Pa) ");
03285 break;
03286 }
03287 case 8:{
03288 fprintf (out,"Capilary pressure (Pa) ");
03289 break;
03290 }
03291 default:{
03292 fprintf (stderr,"\n\n unknown type of component is required in function (%s, line %d).\n",__FILE__,__LINE__);
03293 }
03294 }
03295 }
03296
03297
03298
03299
03300
03301
03302
03303
03304
03305
03306
03307
03308 void glasgowmat::give_dof_names(namevart *dofname, long ntm)
03309 {
03310 if (ntm < 1)
03311 {
03312 print_err("the model defines %ld unknowns but number of transported media is %ld",
03313 __FILE__, __LINE__, __func__, 1, ntm);
03314 abort();
03315 }
03316 dofname[0] = trf_vapor_content;
03317 dofname[1] = trf_gas_press;
03318 dofname[2] = trf_temperature;
03319 }