00001
00002
00003
00004
00005
00006
00007
00008
00009
00010 #include <stdio.h>
00011 #include <stdlib.h>
00012 #include <math.h>
00013 #include "pedersen.h"
00014 #include "globalt.h"
00015
00016 pedmat::pedmat()
00017 {
00018 por = 0.0;
00019 rho = 0.0; ceff = 0.0; chieff = 0.0;
00020 a_0 = 0.0; nn = 0.0; phi_c = 0.0;
00021 delta_dry = 0.0; delta_wet = 0.0;
00022 w_h_sorp = 0.0; n_sorp = 0.0; a_sorp = 0.0;
00023
00024
00025 mw = 18.01528e-3;
00026 gasr = 8.31441;
00027 rhow = 998.0;
00028 awet = 0.0; bwet = 0.0;
00029 k_wg = 0.0; ak = 0.0; bk = 0.0; nk = 0.0;
00030 w_cr = 0.0; w_cap = 0.0; w_vac = 0.0; w_98 = 0.0;
00031 dhvap = 2.7e+5;
00032
00033
00034
00035
00036 gf=NULL;
00037 s_type = 0;
00038 }
00039
00040 pedmat::~pedmat()
00041 {}
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056 void pedmat::matcond (matrix &d,long ri,long ci,long ipp)
00057 {
00058 long n;
00059 n = d.n;
00060
00061 switch (n){
00062 case 1:{
00063 matcond1d (d,ri,ci,ipp);
00064 break;
00065 }
00066 case 2:{
00067 matcond2d (d,ri,ci,ipp);
00068 break;
00069 }
00070 case 3:{
00071 matcond3d (d,ri,ci,ipp);
00072 break;
00073 }
00074 default:{
00075 fprintf (stderr,"\n unknown number of components of stress tensor is required");
00076 fprintf (stderr,"\n in function (%s, line %d).\n",__FILE__,__LINE__);
00077 }
00078 }
00079 }
00080
00081
00082
00083
00084
00085
00086
00087
00088
00089
00090 void pedmat::matcond1d (matrix &d,long ri,long ci,long ipp)
00091 {
00092 double k;
00093 double w,t;
00094 k = 0.0;
00095
00096 w = Tm->ip[ipp].av[0];
00097 t = Tm->ip[ipp].av[1];
00098
00099 if((ri == 0) && (ci == 0))
00100 k = perm_ww(w,t);
00101 if((ri == 0) && (ci == 1))
00102 k = perm_wt(w,t);
00103 if((ri == 1) && (ci == 0))
00104 k = perm_tw(w,t);
00105 if((ri == 1) && (ci == 1))
00106 k = perm_tt(w,t);
00107
00108 d[0][0] = k;
00109 }
00110
00111
00112
00113
00114
00115
00116
00117
00118
00119 void pedmat::matcond2d (matrix &d,long ri,long ci,long ipp)
00120 {
00121 double k;
00122 double w,t;
00123 k = 0.0;
00124
00125 w = Tm->ip[ipp].av[0];
00126 t = Tm->ip[ipp].av[1];
00127
00128 if((ri == 0) && (ci == 0))
00129 k = perm_ww(w,t);
00130 if((ri == 0) && (ci == 1))
00131 k = perm_wt(w,t);
00132 if((ri == 1) && (ci == 0))
00133 k = perm_tw(w,t);
00134 if((ri == 1) && (ci == 1))
00135 k = perm_tt(w,t);
00136
00137 fillm(0.0,d);
00138
00139 d[0][0] = k; d[0][1] = 0.0;
00140 d[1][0] = 0.0; d[1][1] = k;
00141 }
00142
00143
00144
00145
00146
00147
00148
00149
00150
00151
00152 void pedmat::matcond3d (matrix &d,long ri,long ci,long ipp)
00153 {
00154 double k;
00155 double w,t;
00156 k = 0.0;
00157
00158 w = Tm->ip[ipp].av[0];
00159 t = Tm->ip[ipp].av[1];
00160
00161 if((ri == 0) && (ci == 0))
00162 k = perm_ww(w,t);
00163 if((ri == 0) && (ci == 1))
00164 k = perm_wt(w,t);
00165 if((ri == 1) && (ci == 0))
00166 k = perm_tw(w,t);
00167 if((ri == 1) && (ci == 1))
00168 k = perm_tt(w,t);
00169
00170 fillm(0.0,d);
00171
00172 d[0][0]=k; d[0][1]=0.0; d[0][2]=0.0;
00173 d[1][0]=0.0; d[1][1]=k; d[1][2]=0.0;
00174 d[2][0]=0.0; d[2][1]=0.0; d[2][2]=k;
00175 }
00176
00177
00178
00179
00180
00181
00182
00183
00184
00185
00186
00187
00188
00189
00190 void pedmat::matcond2 (matrix &d,long ri,long ci,long ipp)
00191 {
00192 long n;
00193 n = d.n;
00194
00195 switch (n){
00196 case 1:{
00197 matcond1d_2 (d,ri,ci,ipp);
00198 break;
00199 }
00200 case 2:{
00201 matcond2d_2 (d,ri,ci,ipp);
00202 break;
00203 }
00204 case 3:{
00205 matcond3d_2 (d,ri,ci,ipp);
00206 break;
00207 }
00208 default:{
00209 fprintf (stderr,"\n unknown number of components of stress tensor is required");
00210 fprintf (stderr,"\n in function (%s, line %d).\n",__FILE__,__LINE__);
00211 }
00212 }
00213 }
00214
00215
00216
00217
00218
00219
00220
00221
00222
00223
00224
00225
00226 void pedmat::matcond1d_2 (matrix &d,long ri,long ci,long ipp)
00227 {
00228 double g,w,t,phi,delta_gw,kwg,pgw,rhogw;
00229
00230 w = Tm->ip[ipp].av[0];
00231 t = Tm->ip[ipp].av[1];
00232
00233 phi = inverse_sorption_isotherm(w);
00234 delta_gw = get_delta_gw(phi,w);
00235 kwg = get_kw_g(w);
00236 pgw = phi*exp(23.5771 - 4042.9/(t-37.58));
00237 rhogw = pgw/t/gasr*mw;
00238
00239
00240 fillm(0.0,d);
00241
00242 if((ri == 0) && (ci == 0)){
00243 g = Tp->gr1;
00244
00245
00246 d[0][0] = -1.0*kwg*rhow*g - delta_gw*rhogw*g;
00247 }
00248 }
00249
00250
00251
00252
00253
00254
00255
00256
00257
00258
00259 void pedmat::matcond2d_2 (matrix &d,long ri,long ci,long ipp)
00260 {
00261 double w,t,phi,delta_gw,kwg,pgw,rhogw;
00262 double *g;
00263
00264 w = Tm->ip[ipp].av[0];
00265 t = Tm->ip[ipp].av[1];
00266
00267 phi = inverse_sorption_isotherm(w);
00268 delta_gw = get_delta_gw(phi,w);
00269 kwg = get_kw_g(w);
00270 pgw = phi*exp(23.5771 - 4042.9/(t-37.58));
00271 rhogw = pgw/t/gasr*mw;
00272
00273 fillm(0.0,d);
00274
00275 if((ri == 0) && (ci == 0)){
00276 g = new double [2];
00277 g[0] = Tp->gr1;
00278 g[1] = Tp->gr2;
00279
00280 d[0][0] = -1.0*kwg*rhow*g[0] - delta_gw*rhogw*g[0];
00281 d[0][1] = -1.0*kwg*rhow*g[1] - delta_gw*rhogw*g[1];
00282
00283 delete [] g;
00284 }
00285 }
00286
00287
00288
00289
00290
00291
00292
00293
00294
00295
00296
00297 void pedmat::matcond3d_2 (matrix &d,long ri,long ci,long ipp)
00298 {
00299 double w,t,phi,delta_gw,kwg,pgw,rhogw;
00300 double *g;
00301
00302 w = Tm->ip[ipp].av[0];
00303 t = Tm->ip[ipp].av[1];
00304
00305 phi = inverse_sorption_isotherm(w);
00306 delta_gw = get_delta_gw(phi,w);
00307 kwg = get_kw_g(w);
00308 pgw = phi*exp(23.5771 - 4042.9/(t-37.58));
00309 rhogw = pgw/t/gasr*mw;
00310
00311 fillm(0.0,d);
00312
00313 if((ri == 0) && (ci == 0)){
00314 g = new double [3];
00315 g[0] = Tp->gr1;
00316 g[1] = Tp->gr2;
00317 g[2] = Tp->gr3;
00318
00319 d[0][0] = -1.0*kwg*rhow*g[0] - delta_gw*rhogw*g[0];
00320 d[0][1] = -1.0*kwg*rhow*g[1] - delta_gw*rhogw*g[1];
00321 d[0][2] = -1.0*kwg*rhow*g[2] - delta_gw*rhogw*g[2];
00322
00323 delete [] g;
00324 }
00325 }
00326
00327
00328
00329
00330
00331
00332
00333
00334
00335
00336
00337
00338
00339 void pedmat::matcap (double &c,long ri,long ci,long ipp)
00340 {
00341 double w,t;
00342 c = 0.0;
00343
00344 w = Tm->ip[ipp].av[0];
00345 t = Tm->ip[ipp].av[1];
00346
00347 if((ri == 0) && (ci == 0))
00348 c = c_ww(w,t);
00349 if((ri == 0) && (ci == 1))
00350 c = c_wt(w,t);
00351 if((ri == 1) && (ci == 0))
00352 c = c_tw(w,t);
00353 if((ri == 1) && (ci == 1))
00354 c = c_tt(w,t);
00355 }
00356
00357
00358
00359
00360
00361
00362
00363
00364
00365
00366
00367 double pedmat::perm_ww(double w,double t)
00368 {
00369 double k_ww,phi,delta_gw,dphi_dw,p_gws;
00370 double kwg;
00371 double help1,help2,help3;
00372
00373 phi = inverse_sorption_isotherm(w);
00374 delta_gw = get_delta_gw(phi,w);
00375 dphi_dw = get_dphi_dw(w);
00376 p_gws = get_p_gws(t);
00377
00378 k_ww = delta_gw*p_gws*dphi_dw;
00379
00380
00381
00382
00383
00384
00385
00386
00387
00388
00389 if (w < w_cap){
00390
00391 kwg = get_kw_g(w);
00392 help1 = pow((w_cap-w)/awet,(1.0/bwet));
00393 help2 = exp(help1);
00394 help3 = bwet/(w-w_cap);
00395
00396
00397
00398
00399 k_ww = k_ww - kwg*help1*help2/help3;
00400 }
00401 else{
00402 w = w_cap;
00403 }
00404
00405 return(k_ww);
00406 }
00407
00408
00409
00410
00411
00412
00413
00414
00415
00416 double pedmat::perm_wt(double w,double t)
00417 {
00418 double k_wt,phi,delta_gw,dpgw_dt;
00419
00420 phi = inverse_sorption_isotherm(w);
00421 delta_gw = get_delta_gw(phi,w);
00422 dpgw_dt = get_dpgw_dt(t,phi);
00423
00424 k_wt = delta_gw*dpgw_dt;
00425
00426
00427
00428
00429
00430
00431
00432
00433
00434
00435 if (w < w_cap){
00436
00437 }
00438 else{
00439 w = w_cap;
00440 }
00441
00442 return(k_wt);
00443 }
00444
00445
00446
00447
00448
00449
00450
00451
00452
00453 double pedmat::perm_tw(double w,double t)
00454 {
00455 double k_tw,phi,delta_gw,dphi_dw,p_gws;
00456
00457
00458 phi = inverse_sorption_isotherm(w);
00459 delta_gw = get_delta_gw(phi,w);
00460 dphi_dw = get_dphi_dw(w);
00461 p_gws = get_p_gws(t);
00462
00463 k_tw = dhvap*delta_gw*p_gws*dphi_dw;
00464
00465 return(k_tw);
00466 }
00467
00468
00469
00470
00471
00472
00473
00474
00475
00476
00477 double pedmat::perm_tt(double w,double t)
00478 {
00479 double k_tt,phi,delta_gw,dpgw_dt;
00480
00481
00482 phi = inverse_sorption_isotherm(w);
00483 delta_gw = get_delta_gw(phi,w);
00484 dpgw_dt = get_dpgw_dt(t,phi);
00485
00486 k_tt = chieff + dhvap*delta_gw*dpgw_dt;
00487
00488 return(k_tt);
00489 }
00490
00491
00492
00493
00494
00495
00496
00497
00498
00499 double pedmat::c_ww (double w,double t)
00500 {
00501 double c;
00502
00503 c=1.0*rho;
00504
00505 return(c);
00506 }
00507
00508
00509
00510
00511
00512
00513
00514
00515
00516 double pedmat::c_wt (double w,double t)
00517 {
00518 double c;
00519
00520 c=0.0;
00521
00522 return(c);
00523 }
00524
00525
00526
00527
00528
00529
00530
00531
00532
00533 double pedmat::c_tw (double w,double t)
00534 {
00535 double phi,dphi_dw,pgw,s,rhogw,ds_dw,c;
00536
00537
00538
00539 phi = inverse_sorption_isotherm(w);
00540 dphi_dw = get_dphi_dw(w);
00541 pgw = phi*exp(23.5771 - 4042.9/(t-37.58));
00542 rhogw = pgw/t/gasr*mw;
00543
00544 s = ((1.0 - por)*rho*w - por*rhogw)/(por*rhow - por*rhogw);
00545 ds_dw = ((1.0 - por)*rho)/(por*rhow - por*rhogw);
00546 c = dhvap*por*rhogw*(phi*ds_dw - (1.0-s)*dphi_dw);
00547
00548
00549
00550 return(c);
00551 }
00552
00553
00554
00555
00556
00557
00558
00559
00560
00561 double pedmat::c_tt (double w,double t)
00562 {
00563 double c;
00564
00565 c = ceff * rho;
00566
00567 return(c);
00568 }
00569
00570
00571 void pedmat::values_correction (vector &nv)
00572 {
00573
00574 moisture_check(nv[0],nv[1],0);
00575 }
00576
00577
00578
00579
00580
00581
00582
00583
00584
00585
00586
00587 void pedmat::moisture_check(double &w,double t,long ipp)
00588 {
00589
00590 if (w > w_cap){
00591 w = w_cap;
00592
00593 }
00594
00595
00596 Tm->ip[ipp].av[0] = w;
00597 }
00598
00599
00600
00601
00602
00603
00604
00605
00606
00607
00608
00609
00610
00611
00612 double pedmat::get_delta_gw(double phi,double w)
00613 {
00614 double delta_gw;
00615
00616
00617
00618
00619
00620 if(phi <= 0.60){
00621 delta_gw = delta_dry;
00622 }
00623 if(phi > 0.60 && phi < 0.98){
00624 delta_gw = delta_dry + ((phi-0.60)/(0.98-phi))*(delta_wet-delta_dry);
00625 }
00626 if(phi >= 0.98 && w <= w_vac){
00627 delta_gw = delta_wet * ((w_cap - w)/(w_vac - w_98));
00628 }
00629 if(w > w_vac){
00630 delta_gw = 0.0;
00631 }
00632
00633 delta_gw = fabs(delta_gw);
00634
00635 return(delta_gw);
00636 }
00637
00638
00639
00640
00641
00642
00643
00644
00645 double pedmat::get_kw_g(double w)
00646 {
00647 double kw_g,a0_k,wc_k;
00648
00649
00650
00651
00652
00653
00654
00655
00656
00657
00658 a0_k = ak;
00659 wc_k = bk;
00660 kw_g = k_wg*(a0_k + (1.0 - a0_k)/(1 + pow((wc_k/w),nk)));
00661
00662 return(kw_g);
00663 }
00664
00665
00666
00667
00668
00669
00670 double pedmat::sorption_isotherm(double phi)
00671 {
00672 double w;
00673
00674 if (s_type == 0)
00675 w = sorption_isotherm_data_get_w(phi);
00676 if (s_type == 1)
00677 w = sorption_isotherm_hansen(w_h_sorp,a_sorp,n_sorp,phi);
00678
00679 return(w);
00680 }
00681
00682
00683
00684
00685
00686
00687 double pedmat::inverse_sorption_isotherm(double w)
00688 {
00689 double phi;
00690
00691 if (s_type == 0)
00692 phi = inverse_sorption_isotherm_data_get_phi(w);
00693 if (s_type == 1)
00694 phi = inverse_sorption_isotherm_hansen(w_h_sorp,a_sorp,n_sorp,w);
00695
00696 return(phi);
00697 }
00698
00699
00700
00701
00702
00703
00704 double pedmat::get_dphi_dw(double w)
00705 {
00706 double dphi_dw;
00707
00708 if (s_type == 0)
00709 dphi_dw = inverse_sorption_isotherm_data_get_derivative(w);
00710 if (s_type == 1)
00711 dphi_dw = get_dphi_dw_hansen(w);
00712
00713 return(dphi_dw);
00714 }
00715
00716
00717
00718
00719
00720
00721
00722
00723
00724
00725
00726
00727
00728
00729 double pedmat::sorption_isotherm_hansen(double w_h_s, double a_s, double n_s, double phi)
00730 {
00731 double w;
00732
00733
00734 w = w_h_s*pow((1.0-log(phi)/a_s),(-1.0/n_s));
00735
00736 return(w);
00737 }
00738
00739 double pedmat::inverse_sorption_isotherm_hansen(double w)
00740 {
00741 double phi;
00742
00743 if (w > w_cap)
00744 w = w_h_sorp;
00745
00746
00747 phi = inverse_sorption_isotherm_hansen(w_h_sorp,a_sorp,n_sorp,w);
00748
00749 return(phi);
00750 }
00751
00752
00753 double pedmat::suction_curve(double s)
00754 {
00755 double w;
00756
00757 if (s <= 0)
00758 w = w_cap;
00759 else
00760 w = w_cap - awet*pow(log(s),bwet);
00761
00762 return(w);
00763 }
00764
00765
00766 double pedmat::inverse_suction_curve(double w,double t)
00767 {
00768 double s,help;
00769
00770 if (w > w_cap)
00771 w = w_cap;
00772
00773 help = pow(((w_cap-w)/awet),(1.0/bwet));
00774 s = exp(help);
00775
00776 return(s);
00777 }
00778
00779
00780
00781
00782
00783
00784
00785
00786
00787
00788
00789
00790
00791 double pedmat::inverse_sorption_isotherm_hansen(double w_h_s, double a_s, double n_s, double w)
00792 {
00793 double phi;
00794
00795 if (w > w_cap)
00796 w = w_h_sorp;
00797
00798
00799 phi = exp(a_s*(1.0-pow((w_h_s/w),(n_s))));
00800
00801 return(phi);
00802 }
00803
00804
00805
00806
00807
00808
00809
00810
00811
00812
00813
00814
00815
00816 double pedmat::get_dphi_dw_hansen(double w)
00817 {
00818
00819 double dphi_dw;
00820
00821 if (w > w_cap)
00822 w = w_h_sorp;
00823
00824 dphi_dw = exp(a_sorp*(1.0-pow((w_h_sorp/w),n_sorp)))*a_sorp*n_sorp*pow(w_h_sorp,n_sorp)*pow(w,(-1.0-n_sorp));
00825
00826 return(dphi_dw);
00827 }
00828
00829
00830
00831
00832
00833
00834
00835
00836
00837
00838
00839
00840 double pedmat::get_p_gws(double t)
00841 {
00842 double p_gws;
00843
00844
00845 p_gws = exp(23.5771 - 4042.9/(t-37.58));
00846
00847 return(p_gws);
00848 }
00849
00850
00851
00852
00853
00854
00855
00856
00857
00858
00859
00860
00861
00862 double pedmat::get_dpgw_dt(double t, double phi)
00863 {
00864 double dp_gw_dt,dp_gws_dt;
00865
00866 dp_gws_dt = exp(23.5771 - 4042.9/(t-37.58))*4042.9/(t-37.58)/(t-37.58);
00867
00868 dp_gw_dt = phi*dp_gws_dt;
00869
00870 return(dp_gw_dt);
00871 }
00872
00873
00874
00875
00876
00877
00878
00879
00880
00881 double pedmat::get_dpc_dw(double w, double t)
00882 {
00883 double dphi_dw,rh,dpc_drh,dpc_dw;
00884
00885 dphi_dw = get_dphi_dw(w);
00886 rh = inverse_sorption_isotherm(w);
00887
00888 dpc_drh = -1.0*rhow/gasr/t*mw/rh;
00889
00890 dpc_dw = dpc_drh*dphi_dw;
00891
00892 return(dpc_dw);
00893 }
00894
00895
00896
00897
00898
00899
00900
00901
00902
00903 double pedmat::get_dpc_dt(double w, double t)
00904 {
00905 double rh,dpc_dt;
00906
00907 rh = inverse_sorption_isotherm(w);
00908
00909 dpc_dt = -1.0*log(rh)/mw*rhow*gasr;
00910
00911 return(dpc_dt);
00912 }
00913
00914
00915
00916
00917
00918
00919
00920
00921
00922
00923
00924
00925
00926
00927
00928 double pedmat::transmission_transcoeff(double trc,long ri,long ci,long nn,long bc,long ipp)
00929 {
00930 long k;
00931 double new_trc,w,t;
00932 new_trc = 0.0;
00933
00934 k=Gtt->give_dof(nn,0);
00935 if (k>0) {w = Lsrst->lhs[k-1]+Lsrst->lhsi[k-1];}
00936 if (k==0) {w = 0.0;}
00937 if (k<0) {w = Tb->lc[0].pv[0-k-1].getval();}
00938 k=Gtt->give_dof(nn,1);
00939 if (k>0) {t = Lsrst->lhs[k-1]+Lsrst->lhsi[k-1];}
00940 if (k==0) {t = 0.0;}
00941 if (k<0) {t = Tb->lc[0].pv[0-k-1].getval();}
00942
00943 if((ri == 0) && (ci == 0))
00944 new_trc = get_transmission_transcoeff_ww(w,t,bc,ipp);
00945 if((ri == 0) && (ci == 1))
00946 new_trc = 0.0;
00947
00948 if((ri == 1) && (ci == 0))
00949 new_trc = 0.0;
00950 if((ri == 1) && (ci == 1))
00951 new_trc = get_transmission_transcoeff_tt(w,t,bc,ipp);
00952
00953 new_trc = new_trc*trc;
00954
00955 return (new_trc);
00956 }
00957
00958
00959
00960
00961
00962
00963
00964
00965
00966
00967
00968
00969
00970
00971 double pedmat::transmission_nodval(double nodval,double trc2,long ri,long ci,long nn,long bc,long ipp)
00972 {
00973 long k;
00974 double new_nodval,w,t;
00975 new_nodval = 0.0;
00976
00977 k=Gtt->give_dof(nn,0);
00978 if (k>0) {w = Lsrst->lhs[k-1]+Lsrst->lhsi[k-1];}
00979 if (k==0) {w = 0.0;}
00980 if (k<0) {w = Tb->lc[0].pv[0-k-1].getval();}
00981 k=Gtt->give_dof(nn,1);
00982 if (k>0) {t = Lsrst->lhs[k-1]+Lsrst->lhsi[k-1];}
00983 if (k==0) {t = 0.0;}
00984 if (k<0) {t = Tb->lc[0].pv[0-k-1].getval();}
00985
00986 if((ri == 0) && (ci == 0))
00987 new_nodval = get_transmission_nodval_ww(nodval,w,t,bc,ipp);
00988 if((ri == 0) && (ci == 1))
00989 new_nodval = 0.0;
00990
00991 if((ri == 1) && (ci == 0))
00992 new_nodval = 0.0;
00993 if((ri == 1) && (ci == 1))
00994 new_nodval = get_transmission_nodval_tt(nodval,w,t,bc,ipp);
00995
00996 return (new_nodval);
00997 }
00998
00999
01000
01001
01002
01003
01004
01005
01006
01007
01008
01009
01010
01011
01012
01013 double pedmat::transmission_flux(double nodval,double trc2,long ri,long ci,long nn,long bc,long ipp)
01014 {
01015 long k;
01016 double flux,w,t;
01017 flux = 0.0;
01018
01019 k=Gtt->give_dof(nn,0);
01020 if (k>0) {w = Lsrst->lhs[k-1]+Lsrst->lhsi[k-1];}
01021 if (k==0) {w = 0.0;}
01022 if (k<0) {w = Tb->lc[0].pv[0-k-1].getval();}
01023 k=Gtt->give_dof(nn,1);
01024 if (k>0) {t = Lsrst->lhs[k-1]+Lsrst->lhsi[k-1];}
01025 if (k==0) {t = 0.0;}
01026 if (k<0) {t = Tb->lc[0].pv[0-k-1].getval();}
01027
01028 if((ri == 0) && (ci == 0))
01029 flux = get_transmission_flux_ww(nodval,w,t,bc,ipp);
01030 if((ri == 0) && (ci == 1))
01031 flux = 0.0;
01032
01033 if((ri == 1) && (ci == 0))
01034 flux = 0.0;
01035 if((ri == 1) && (ci == 1))
01036 flux = get_transmission_flux_tt(nodval,w,t,bc,ipp);
01037
01038 return (flux);
01039 }
01040
01041
01042
01043
01044
01045
01046
01047
01048
01049 double pedmat::get_transmission_transcoeff_ww(double w,double t,long bc,long ipp)
01050 {
01051 double trc;
01052 double pgws;
01053
01054 switch (bc){
01055 case 30:{
01056
01057 pgws = get_p_gws(t);
01058
01059 trc = pgws*get_dphi_dw(w);
01060 break;
01061 }
01062 case 31:{
01063
01064 pgws = get_p_gws(t);
01065
01066 trc = pgws*get_dphi_dw(w);
01067 break;
01068 }
01069 case 32:{
01070 trc = 0.0;
01071 break;
01072 }
01073 case 33:{
01074 trc = 0.0;
01075 break;
01076 }
01077 default:{
01078 fprintf (stderr,"\n\n No real boundary condition is prescribed (%s, line %d).\n",__FILE__,__LINE__);
01079 exit(0);
01080 }
01081 }
01082
01083 return(trc);
01084 }
01085
01086
01087
01088
01089
01090
01091
01092
01093
01094 double pedmat::get_transmission_nodval_ww(double bv,double w,double t,long bc,long ipp)
01095 {
01096 double nodval;
01097 double pgws,rh;
01098
01099 switch (bc){
01100 case 30:{
01101
01102 pgws = get_p_gws(t);
01103 nodval = bv*pgws;
01104 nodval = nodval - pgws*inverse_sorption_isotherm(w) + pgws*w*get_dphi_dw(w);
01105 break;
01106 }
01107 case 31:{
01108 pgws = get_p_gws(t);
01109 nodval = bv;
01110 nodval = nodval - pgws*inverse_sorption_isotherm(w) + pgws*w*get_dphi_dw(w);
01111 break;
01112 }
01113 case 32:{
01114
01115 pgws = get_p_gws(t);
01116 rh = inverse_sorption_isotherm(w);
01117
01118 nodval = pgws*rh;
01119 bv = pgws*bv;
01120 nodval = bv - nodval;
01121 break;
01122 }
01123 case 33:{
01124
01125 pgws = get_p_gws(t);
01126 rh = inverse_sorption_isotherm(w);
01127
01128 nodval = pgws*rh;
01129 nodval = bv - nodval;
01130 break;
01131 }
01132 default:{
01133 fprintf (stderr,"\n\n No real boundary condition is prescribed (%s, line %d).\n",__FILE__,__LINE__);
01134 exit(0);
01135 }
01136 }
01137 return(nodval);
01138 }
01139
01140
01141
01142
01143
01144
01145
01146
01147
01148
01149 double pedmat::get_transmission_flux_ww(double bv,double w,double t,long bc,long ipp)
01150 {
01151 double flux;
01152 double pgws,rh;
01153
01154 switch (bc){
01155 case 30:{
01156
01157 pgws = get_p_gws(t);
01158 rh = inverse_sorption_isotherm(w);
01159
01160 flux = pgws*rh;
01161 bv = pgws*bv;
01162 flux = bv - flux;
01163 break;
01164 }
01165 case 31:{
01166
01167 pgws = get_p_gws(t);
01168 rh = inverse_sorption_isotherm(w);
01169
01170 flux = pgws*rh;
01171 flux = bv - flux;
01172 break;
01173 }
01174 case 32:{
01175
01176 pgws = get_p_gws(t);
01177 rh = inverse_sorption_isotherm(w);
01178
01179 flux = pgws*rh;
01180 bv = pgws*bv;
01181 flux = bv - flux;
01182 break;
01183 }
01184 case 33:{
01185
01186 pgws = get_p_gws(t);
01187 rh = inverse_sorption_isotherm(w);
01188
01189 flux = pgws*rh;
01190 flux = bv - flux;
01191 break;
01192 }
01193 default:{
01194 fprintf (stderr,"\n\n No real boundary condition is prescribed (%s, line %d).\n",__FILE__,__LINE__);
01195 exit(0);
01196 }
01197 }
01198 return(flux);
01199 }
01200
01201
01202
01203
01204
01205
01206
01207
01208
01209 double pedmat::get_transmission_transcoeff_tt(double w,double t,long bc,long ipp)
01210 {
01211 double trc;
01212
01213 switch (bc){
01214 case 30:{
01215 trc = 1.0;
01216 break;
01217 }
01218 case 31:{
01219 trc = 0.0;
01220 break;
01221 }
01222 default:{
01223 fprintf (stderr,"\n\n No real boundary condition is prescribed (%s, line %d).\n",__FILE__,__LINE__);
01224 exit(0);
01225 }
01226 }
01227
01228 return(trc);
01229 }
01230
01231
01232
01233
01234
01235
01236
01237
01238
01239 double pedmat::get_transmission_nodval_tt(double bv,double w,double t,long bc,long ipp)
01240 {
01241 double nodval;
01242
01243 switch (bc){
01244 case 30:{
01245 nodval = bv;
01246 break;
01247 }
01248 case 31:{
01249 nodval = (bv - t);
01250 break;
01251 }
01252 default:{
01253 fprintf (stderr,"\n\n No real boundary condition is prescribed (%s, line %d).\n",__FILE__,__LINE__);
01254 exit(0);
01255 }
01256 }
01257 return(nodval);
01258 }
01259
01260
01261
01262
01263
01264
01265
01266
01267
01268
01269 double pedmat::get_transmission_flux_tt(double bv,double w,double t,long bc,long ipp)
01270 {
01271 double flux;
01272
01273 switch (bc){
01274 case 30:{
01275 flux = (bv - t);
01276 break;
01277 }
01278 case 31:{
01279 flux = (bv - t);
01280 break;
01281 }
01282 default:{
01283 fprintf (stderr,"\n\n No real boundary condition is prescribed (%s, line %d).\n",__FILE__,__LINE__);
01284 exit(0);
01285 }
01286 }
01287 return(flux);
01288 }
01289
01290
01291
01292
01293
01294
01295
01296
01297
01298
01299
01300 double pedmat::get_othervalue(long compother,double w,double t)
01301 {
01302 double other;
01303
01304 switch (compother){
01305 case 0:{
01306 other = w;
01307 break;
01308 }
01309 case 1:{
01310 other = inverse_sorption_isotherm(w);
01311 break;
01312 }
01313 case 2:{
01314 other = t;
01315 break;
01316 }
01317 case 3:{
01318 other = inverse_sorption_isotherm(w) * get_p_gws(t);
01319 break;
01320 }
01321 case 4:{
01322 double pc,pgw,rh;
01323
01324 pgw = inverse_sorption_isotherm(w) * get_p_gws(t);
01325 rh = inverse_sorption_isotherm(w);
01326 pc = -rhow*gasr*t/mw*log(rh);
01327 other = pgw - pc;
01328 break;
01329 }
01330 case 5:{
01331 double rh;
01332
01333 rh = inverse_sorption_isotherm(w);
01334 other = -rhow*gasr*t/mw*log(rh);
01335 break;
01336 }
01337 case 6:{
01338 other = inverse_suction_curve(w,t);
01339 break;
01340 }
01341 default:{
01342 fprintf (stderr,"\n\n unknown type of component is required in function (%s, line %d).\n",__FILE__,__LINE__);
01343 }
01344 }
01345 return (other);
01346
01347 }
01348
01349
01350
01351
01352
01353
01354 void pedmat::print_othervalue_name(FILE *out,long compother)
01355 {
01356 switch (compother){
01357 case 0:{
01358 fprintf (out,"Moisture content u (kg/kg)");
01359 break;
01360 }
01361 case 1:{
01362 fprintf (out,"Relative humidity () ");
01363 break;
01364 }
01365 case 2:{
01366 fprintf (out,"Temperature (K) ");
01367 break;
01368 }
01369 case 3:{
01370 fprintf (out,"Water vapour pressure (Pa)");
01371 break;
01372 }
01373 case 4:{
01374 fprintf (out,"Water pressure (Pa) ");
01375 break;
01376 }
01377 case 5:{
01378 fprintf (out,"Capillary pressure (Pa) ");
01379 break;
01380 }
01381 case 6:{
01382 fprintf (out,"Suction pressure (Pa) ");
01383 break;
01384 }
01385 default:{
01386 fprintf (stderr,"\n\n unknown type of component is required in function (%s, line %d).\n",__FILE__,__LINE__);
01387 }
01388 }
01389 }
01390
01391
01392
01393
01394
01395
01396
01397 void pedmat::read(XFILE *in)
01398 {
01399
01400
01401
01402
01403 xfscanf (in,"%lf %lf %lf %lf %lf %lf %d", &por, &rho, &ceff, &chieff, &delta_dry, &delta_wet, &s_type);
01404
01405 if(s_type == 0)
01406 sorption_isotherm_data_read(in);
01407 else
01408 xfscanf (in,"%lf %lf %lf", &w_h_sorp, &a_sorp, &n_sorp);
01409
01410 xfscanf (in," %lf %lf %lf %lf %lf %lf %lf %lf %lf", &awet, &bwet, &w_cr, &w_cap, &w_vac, &k_wg, &ak, &bk, &nk);
01411
01412 w_98 = sorption_isotherm(0.98);
01413 }
01414
01415
01416
01417
01418
01419
01420
01421 void pedmat::print(FILE *out)
01422 {
01423
01424
01425
01426
01427 fprintf (out," %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e ", por,rho, ceff, chieff, delta_dry, delta_wet, w_h_sorp, a_sorp, n_sorp, awet, bwet, w_cr, w_cap, w_vac, k_wg, ak, bk, nk);
01428
01429 }
01430
01431
01432
01433
01434
01435
01436
01437
01438
01439
01440
01441 double pedmat::sorption_isotherm_data_get_w(double phi)
01442 {
01443 double help,w;
01444
01445 w = gf->getval3(phi,help);
01446
01447 return w;
01448 }
01449
01450
01451
01452
01453
01454
01455
01456 double pedmat::inverse_sorption_isotherm_data_get_phi(double w)
01457 {
01458 double phi;
01459
01460 phi = gf->getval (w);
01461
01462 return phi;
01463 }
01464
01465
01466
01467
01468
01469
01470
01471
01472 double pedmat::inverse_sorption_isotherm_data_get_derivative(double w)
01473 {
01474 double dphi_dw,phi;
01475
01476 phi = gf->getval2(w, dphi_dw);
01477
01478 return dphi_dw;
01479 }
01480
01481
01482
01483
01484
01485
01486
01487 void pedmat::sorption_isotherm_data_read(XFILE *in)
01488 {
01489 gf = new tablefunct;
01490 gf->read(in);
01491 }
01492
01493
01494
01495
01496
01497
01498
01499
01500 void pedmat::sorption_isotherm_data_print(FILE *out)
01501 {
01502 if (gf != NULL)
01503 gf->print (out);
01504 }
01505