00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029 #include "kunzel.h"
00030 #include "globalt.h"
00031
00032
00033 kunmat::kunmat()
00034 {
00035
00036 kunzeltype = 0;
00037
00038
00039 rho_m = 0.0;
00040
00041 rho_w = 1000.0;
00042
00043 M = 0.01801528;
00044
00045 R = 8.314472;
00046 }
00047
00048 kunmat::~kunmat()
00049 {
00050
00051 }
00052
00053
00054
00055
00056
00057
00058
00059
00060 void kunmat::read(XFILE *in)
00061 {
00062 long i,nr,nc;
00063
00064 for (i = 2; i<20; i++){
00065 xfscanf(in, "%d",&MatChar[i]);
00066 }
00067
00068 for (i = 2;i < 19; i++){
00069
00070 switch (MatChar[i]){
00071 case 0:{
00072
00073 break;
00074 }
00075 case 1:{
00076
00077
00078 data[i]=new gfunct ();
00079
00080
00081 data[i]->tfunc=stat;
00082
00083 xfscanf(in, "%lf",&data[i]->f);
00084
00085 break;
00086 }
00087 case 2:{
00088
00089
00090
00091
00092 xfscanf(in, "%ld %ld",&nr,&nc);
00093
00094
00095 data[i]=new gfunct (tab,nr);
00096
00097 data[i]->tabf[0].itype = piecewiselin;
00098
00099 data[i]->tabf[0].readval (in);
00100
00101
00102 xfscanf(in, "%d",&nr);
00103
00104 break;
00105 }
00106 case 3:{
00107
00108 break;
00109 }
00110 case 30:{
00111 if( MatChar[i-1] ==1)
00112 {
00113 xfscanf(in, "%lf %lf",&MatFunce[i][0],&MatFunce[i][1]);
00114 }
00115 else
00116 {
00117 xfscanf(in, "%lf %lf",&MatFunce[i][0],&MatFunce[i][1]);
00118 }
00119 break;
00120 }
00121 case 31:{
00122 if( MatChar[i-1] ==1)
00123 {
00124 if (i == 6)
00125 {
00126 xfscanf(in, "%lf %lf %lf",&MatFunce[i][0],&MatFunce[i][1], &MatFunce[i][2]);
00127 }
00128 else
00129 {
00130 xfscanf(in, "%lf %lf",&MatFunce[i][0],&MatFunce[i][1]);
00131 }
00132 }
00133 else
00134 {
00135 if (i == 6)
00136 {
00137 xfscanf(in, "%lf %lf %lf",&MatFunce[i][0],&MatFunce[i][1], &MatFunce[i][2]);
00138 }
00139 else
00140 {
00141 xfscanf(in, "%lf %lf",&MatFunce[i][0],&MatFunce[i][1]);
00142 }
00143 }
00144
00145 break;
00146 }
00147 case 32:{
00148 if( MatChar[i-1] ==1)
00149 {
00150 xfscanf(in, "%lf %lf",&MatFunce[i][0],&MatFunce[i][1]);
00151 }
00152 else
00153 {
00154 xfscanf(in, "%lf %lf",&MatFunce[i][0],&MatFunce[i][1]);
00155 }
00156 break;
00157 }
00158 default:{
00159 print_err("unknown definition of Material KUNZEL is required",__FILE__,__LINE__,__func__);
00160 }
00161 }
00162 }
00163
00164
00165 kunzeltype = MatChar[19];
00166
00167 CorD(2,kd,0,rho_m,a2,a3);
00168 }
00169
00170
00171
00172
00173
00174
00175
00176
00177 void kunmat::print(FILE *out)
00178 {
00179 long i,j;
00180
00181 fprintf (out,"\n");
00182 for (i = 2; i<20; i++)
00183 {
00184 fprintf(out, " %ld ",MatChar[i]);
00185 }
00186 fprintf(out, " %ld ",MatChar[19]);
00187 for (i = 2;i < 19; i++)
00188 {
00189
00190 switch (MatChar[i]){
00191 case 0:{
00192
00193 break;
00194 }
00195 case 1:{
00196
00197 fprintf(out, "\n %e ",data[i]->f);
00198 break;
00199 }
00200 case 2:{
00201
00202
00203 fprintf(out, "\n %ld %d ",data[i]->tabf[0].asize,2);
00204
00205 for (j=0;j<data[i]->tabf[0].asize;j++){
00206 fprintf(out, "\n %e %e ",data[i]->tabf[0].x[j],data[i]->tabf[0].y[j]);
00207 }
00208
00209 fprintf(out, "\n %ld ",data[i]->tabf[0].asize);
00210 break;
00211 }
00212 case 3:{
00213
00214 break;
00215 }
00216 case 30:{
00217 if( MatChar[i-1] ==1)
00218 {
00219 fprintf(out, "\n %e %e ",MatFunce[i][0],MatFunce[i][1]);
00220 }
00221 else
00222 {
00223 fprintf(out, "\n %e %e ",MatFunce[i][0],MatFunce[i][1]);
00224 }
00225 break;
00226 }
00227 case 31:{
00228 if( MatChar[i-1] ==1)
00229 {
00230 if (i == 6)
00231 {
00232 fprintf(out, "\n %e %e %e ",MatFunce[i][0],MatFunce[i][1], MatFunce[i][2]);
00233 }
00234 else
00235 {
00236 fprintf(out, "\n %e %e ",MatFunce[i][0],MatFunce[i][1]);
00237 }
00238 }
00239 else
00240 {
00241 if (i == 6)
00242 {
00243 fprintf(out, "\n %e %e %e",MatFunce[i][0],MatFunce[i][1], MatFunce[i][2]);
00244 }
00245 else
00246 {
00247 fprintf(out, "\n %e %e ",MatFunce[i][0],MatFunce[i][1]);
00248 }
00249 }
00250
00251 break;
00252 }
00253 case 32:{
00254 if( MatChar[i-1] ==1)
00255 {
00256 fprintf(out, "\n %e %e ",MatFunce[i][0],MatFunce[i][1]);
00257 }
00258 else
00259 {
00260 fprintf(out, "\n %e %e ",MatFunce[i][0],MatFunce[i][1]);
00261 }
00262 break;
00263 }
00264 default:{
00265 print_err("unknown definition of Material KUNZEL is required",__FILE__,__LINE__,__func__);
00266 fprintf (out,"\n");
00267 }
00268 }
00269 }
00270 fprintf (out,"\n");
00271 }
00272
00273
00274
00275
00276
00277
00278
00279
00280
00281
00282
00283
00284
00285 void kunmat::matcond (matrix &d,long ri,long ci,long ipp)
00286 {
00287 long n;
00288 n = d.n;
00289
00290 switch (n){
00291 case 1:{
00292 matcond1d (d,ri,ci,ipp);
00293 break;
00294 }
00295 case 2:{
00296 matcond2d (d,ri,ci,ipp);
00297 break;
00298 }
00299 case 3:{
00300 matcond3d (d,ri,ci,ipp);
00301 break;
00302 }
00303 default:{
00304 print_err("unknown number of components of conductivity tensor is required",__FILE__,__LINE__,__func__);
00305 }
00306 }
00307 }
00308
00309
00310
00311
00312
00313
00314
00315
00316
00317
00318
00319
00320 void kunmat::matcond1d (matrix &d,long ri,long ci,long ipp)
00321 {
00322 double k;
00323 k = 0.0;
00324
00325 if((ri == 0) && (ci == 0))
00326 k = kmm(ipp);
00327 if((ri == 0) && (ci == 1))
00328 k = kmt(ipp);
00329 if((ri == 1) && (ci == 0))
00330 k = khm(ipp);
00331 if((ri == 1) && (ci == 1))
00332 k = kht(ipp);
00333
00334 d[0][0] = k;
00335 }
00336
00337
00338
00339
00340
00341
00342
00343
00344
00345
00346
00347 void kunmat::matcond2d (matrix &d,long ri,long ci,long ipp)
00348 {
00349 double k;
00350 k = 0.0;
00351
00352
00353 if((ri == 0) && (ci == 0))
00354 k = kmm(ipp);
00355 if((ri == 0) && (ci == 1))
00356 k = kmt(ipp);
00357 if((ri == 1) && (ci == 0))
00358 k = khm(ipp);
00359 if((ri == 1) && (ci == 1))
00360 k = kht(ipp);
00361
00362 fillm(0.0,d);
00363
00364
00365
00366 d[0][0] = k; d[0][1] = 0.0;
00367 d[1][0] = 0.0; d[1][1] = k;
00368 }
00369
00370
00371
00372
00373
00374
00375
00376
00377
00378
00379
00380 void kunmat::matcond3d (matrix &d,long ri,long ci,long ipp)
00381 {
00382 double k;
00383 k = 0.0;
00384
00385
00386 if((ri == 0) && (ci == 0))
00387 k = kmm(ipp);
00388 if((ri == 0) && (ci == 1))
00389 k = kmt(ipp);
00390 if((ri == 1) && (ci == 0))
00391 k = khm(ipp);
00392 if((ri == 1) && (ci == 1))
00393 k = kht(ipp);
00394
00395 fillm(0.0,d);
00396
00397 d[0][0]=k; d[0][1]=0.0; d[0][2]=0.0;
00398 d[1][0]=0.0; d[1][1]=k; d[1][2]=0.0;
00399 d[2][0]=0.0; d[2][1]=0.0; d[2][2]=k;
00400 }
00401
00402
00403
00404
00405
00406
00407
00408
00409
00410
00411
00412 void kunmat::matcap (double &c,long ri,long ci,long ipp)
00413 {
00414 c=0.0;
00415
00416 if((ri == 0) && (ci == 0))
00417 c = cmm(ipp);
00418 if((ri == 0) && (ci == 1))
00419 c = cmt(ipp);
00420 if((ri == 1) && (ci == 0))
00421 c = chm(ipp);
00422 if((ri == 1) && (ci == 1)){
00423 c = cht(ipp);
00424 }
00425 }
00426
00427
00428
00429
00430
00431
00432
00433
00434
00435 void kunmat::values_correction (vector &nv)
00436 {
00437
00438 switch (kunzeltype) {
00439 case 0:{
00440 double rh, tk;
00441
00442 rh=nv[0];
00443
00444 tk=nv[1];
00445
00446 if (rh >= 1.0)
00447 rh = 1.0;
00448 if (rh <= 0.0)
00449 rh = 0.0;
00450
00451 if (tk >= 350.0)
00452 tk = 350.0;
00453 if (tk <= 240.0)
00454 tk = 240.0;
00455
00456 nv[0]=rh;
00457 nv[1]=tk;
00458 break;
00459 }
00460 case 1:{
00461 double pv, tk;
00462
00463 pv=nv[0];
00464
00465 tk=nv[1];
00466
00467
00468
00469
00470
00471
00472 if (tk >= 350.0)
00473 tk = 350.0;
00474 if (tk <= 240.0)
00475 tk = 240.0;
00476
00477 nv[0]=pv;
00478 nv[1]=tk;
00479 break;
00480 }
00481 default:{
00482 print_err("unknown type of Kunzel is required",__FILE__,__LINE__,__func__);
00483
00484 }
00485 }
00486
00487
00488
00489 double rh,t;
00490
00491
00492 }
00493
00494
00495
00496
00497
00498
00499
00500
00501
00502 double kunmat::kmm(long ipp)
00503 {
00504 double kk;
00505
00506
00507 switch (kunzeltype) {
00508 case 0:{
00509 double rh, tk;
00510 double dw, mi, dp, dmretc, ps, dphi, moist, smc;
00511
00512
00513 rh = Tm->ip[ipp].av[0];
00514
00515
00516 tk = Tm->ip[ipp].av[1];
00517
00518 if (rh>1.0) rh = 1.0;
00519 if (rh<0.0) rh = 0.0;
00520
00521
00522 moist = Tm->ip[ipp].eqother[0];
00523
00524
00525 CorD(4,kd,moist,mi,a2,a3);
00526
00527
00528 dw = Tm->ip[ipp].eqother[3];
00529
00530 dp = water_vapour_permeability(tk,ipp);
00531 ps = saturated_water_vapor_pressure( tk);
00532
00533
00534 dmretc = Tm->ip[ipp].eqother[1];
00535
00536 dphi = dw * dmretc;
00537
00538 kk = dphi + dp * ps;
00539
00540 break;
00541 }
00542 case 1:{
00543 double pv,tk,dw, dp;
00544
00545
00546 pv = Tm->ip[ipp].av[0];
00547
00548
00549 tk = Tm->ip[ipp].av[1];
00550
00551
00552 dw = Tm->ip[ipp].eqother[3];
00553
00554 dp = water_vapour_permeability(tk,ipp);
00555
00556 kk = dp+dw*M/(R*tk);
00557 break;
00558 }
00559 default:{
00560 print_err("unknown type of Kunzel is required",__FILE__,__LINE__,__func__);
00561
00562 }
00563 }
00564
00565
00566
00567
00568 return (kk);
00569 }
00570
00571
00572
00573
00574
00575
00576
00577
00578
00579
00580
00581 double kunmat::kmt(long ipp)
00582 {
00583 double kk;
00584
00585 switch (kunzeltype) {
00586 case 0:{
00587 double rh,tk;
00588 double dp, dpvs;
00589
00590
00591 rh = Tm->ip[ipp].av[0];
00592
00593
00594 tk = Tm->ip[ipp].av[1];
00595
00596 if (rh>1.0) rh = 1.0;
00597 if (rh<0.0) rh = 0.0;
00598
00599 dpvs = derivative_of_saturated_water_vapor_pressure(tk);
00600 dp = water_vapour_permeability(tk,ipp);
00601
00602
00603 kk = dp * rh * dpvs;
00604
00605 break;
00606 }
00607 case 1:{
00608 double pv,tk, dw;
00609
00610 pv = Tm->ip[ipp].av[0];
00611
00612 tk = Tm->ip[ipp].av[1];
00613
00614
00615 dw = Tm->ip[ipp].eqother[3];
00616
00617 kk = -dw*M*pv/(R*tk*tk);
00618
00619 break;
00620 }
00621
00622 default:{
00623 print_err("unknown type of Kunzel is required",__FILE__,__LINE__,__func__);
00624
00625 }
00626
00627 }
00628
00629
00630
00631 return (kk);
00632 }
00633
00634
00635
00636
00637
00638
00639
00640
00641
00642
00643 double kunmat::kht(long ipp)
00644 {
00645 double kk;
00646
00647 switch (kunzeltype) {
00648 case 0:{
00649 double rh, tk;
00650 double lv1, lambda, dpvs, dp, moist;
00651
00652
00653 rh = Tm->ip[ipp].av[0];
00654
00655
00656 tk = Tm->ip[ipp].av[1];
00657
00658 if (rh>1.0) rh = 1.0;
00659 if (rh<0.0) rh = 0.0;
00660
00661
00662 moist = Tm->ip[ipp].eqother[0];
00663
00664 CorD(10,kd,moist,lambda,a2,a3);
00665
00666 dpvs = derivative_of_saturated_water_vapor_pressure(tk);
00667 dp = water_vapour_permeability(tk,ipp);
00668
00669 lv1 = latent_heat_of_evaporation_of_water (tk);
00670
00671
00672 kk = lambda + lv1 * dp * rh * dpvs;
00673
00674 break;
00675 }
00676 case 1:{
00677 double pv, tk, moist, lambda;
00678
00679
00680 pv = Tm->ip[ipp].av[0];
00681
00682
00683 tk = Tm->ip[ipp].av[1];
00684
00685
00686 moist = Tm->ip[ipp].eqother[0];
00687
00688
00689 CorD(10,kd,moist,lambda,a2,a3);
00690
00691 kk = lambda;
00692 break;
00693 }
00694 default:{
00695 print_err("unknown type of Kunzel is required",__FILE__,__LINE__,__func__);
00696
00697 }
00698 }
00699
00700 return (kk);
00701 }
00702
00703
00704
00705
00706
00707
00708
00709
00710
00711
00712
00713 double kunmat::khm(long ipp)
00714 {
00715 double kk;
00716
00717 switch (kunzeltype) {
00718 case 0:{
00719 double rh, tk;
00720
00721 double lv, dp, ps;
00722
00723 rh = Tm->ip[ipp].av[0];
00724
00725
00726 tk = Tm->ip[ipp].av[1];
00727
00728 if (rh>1.0) rh = 1.0;
00729 if (rh<0.0) rh = 0.0;
00730
00731 dp = water_vapour_permeability(tk,ipp);
00732 ps = saturated_water_vapor_pressure (tk);
00733 lv = latent_heat_of_evaporation_of_water (tk);
00734
00735
00736 kk = lv * dp * ps;
00737
00738 break;
00739 }
00740 case 1:{
00741 double pv,tk;
00742 double lv, dp, ps;
00743
00744
00745 pv = Tm->ip[ipp].av[0];
00746
00747
00748 tk = Tm->ip[ipp].av[1];
00749
00750 dp = water_vapour_permeability(tk,ipp);
00751 lv = latent_heat_of_evaporation_of_water (tk);
00752
00753 kk = lv*dp;
00754 break;
00755 }
00756 default:{
00757 print_err("unknown type of Kunzel is required",__FILE__,__LINE__,__func__);
00758
00759 }
00760 }
00761
00762
00763
00764 return (kk);
00765 }
00766
00767
00768
00769
00770
00771
00772
00773
00774
00775
00776 double kunmat::cmm(long ipp)
00777 {
00778 double cc;
00779
00780 switch (kunzeltype) {
00781 case 0:{
00782
00783 cc = Tm->ip[ipp].eqother[1];
00784
00785 break;
00786 }
00787 case 1:{
00788 double pv, tk;
00789
00790 pv = Tm->ip[ipp].av[0];
00791
00792 tk = Tm->ip[ipp].av[1];
00793
00794 cc = M/(R*tk);
00795
00796 break;
00797 }
00798
00799 default:{
00800 print_err("unknown type of Kunzel is required",__FILE__,__LINE__,__func__);
00801
00802 }
00803
00804 }
00805
00806
00807
00808 return (cc);
00809 }
00810
00811
00812
00813
00814
00815
00816
00817
00818
00819
00820 double kunmat::cmt(long ipp)
00821 {
00822 double cc;
00823
00824 switch (kunzeltype) {
00825 case 0:{
00826 double rh,tk;
00827 double dp, dpvs;
00828
00829
00830 rh = Tm->ip[ipp].av[0];
00831
00832
00833 tk = Tm->ip[ipp].av[1];
00834
00835 cc = 0.0;
00836
00837 break;
00838 }
00839 case 1:{
00840 double pv,tk, dw;
00841
00842 pv = Tm->ip[ipp].av[0];
00843
00844 tk = Tm->ip[ipp].av[1];
00845
00846 cc = -M*pv/(R*tk*tk);
00847 break;
00848 }
00849
00850 default:{
00851 print_err("unknown type of Kunzel is required",__FILE__,__LINE__,__func__);
00852 }
00853
00854 }
00855
00856
00857
00858 return (cc);
00859 }
00860
00861
00862
00863
00864
00865
00866
00867
00868
00869
00870 double kunmat::chm(long ipp)
00871 {
00872 double cc;
00873
00874 switch (kunzeltype) {
00875 case 0:{
00876 double rh,tk;
00877 double dp, dpvs;
00878
00879
00880 rh = Tm->ip[ipp].av[0];
00881
00882
00883 tk = Tm->ip[ipp].av[1];
00884
00885 cc = 0.0;
00886
00887 break;
00888 }
00889 case 1:{
00890 double pv,tk;
00891
00892 pv = Tm->ip[ipp].av[0];
00893
00894 tk = Tm->ip[ipp].av[1];
00895
00896 cc = 0.0;
00897
00898 break;
00899 }
00900
00901 default:{
00902 print_err("unknown type of Kunzel is required",__FILE__,__LINE__,__func__);
00903
00904 }
00905
00906 }
00907
00908
00909
00910 return (cc);
00911 }
00912
00913
00914
00915
00916
00917
00918
00919
00920
00921
00922 double kunmat::cht(long ipp)
00923 {
00924 double cc;
00925
00926 switch (kunzeltype) {
00927 case 0:{
00928 double rh,tk;
00929 double dp, dpvs;
00930
00931
00932 rh = Tm->ip[ipp].av[0];
00933
00934
00935 tk = Tm->ip[ipp].av[1];
00936
00937 if (rh>1.0) rh = 1.0;
00938 if (rh<0.0) rh = 0.0;
00939
00940 cc = derivative_of_the_enthalpy_density(ipp);
00941
00942 break;
00943 }
00944 case 1:{
00945 double pv,tk;
00946
00947 pv = Tm->ip[ipp].av[0];
00948
00949 tk = Tm->ip[ipp].av[1];
00950
00951 cc = derivative_of_the_enthalpy_density(ipp);
00952
00953 break;
00954 }
00955
00956 default:{
00957 print_err("unknown type of Kunzel is required",__FILE__,__LINE__,__func__);
00958 }
00959
00960 }
00961
00962
00963
00964 return (cc);
00965 }
00966
00967
00968
00969
00970
00971
00972
00973
00974 double kunmat::derivative_of_saturated_water_vapor_pressure(double tk)
00975 {
00976
00977 return ((4042.9 * exp(23.5771-4042.9/(tk - 37.58)))/((tk-37.58)*(tk-37.58)));
00978 }
00979
00980
00981
00982
00983
00984
00985
00986
00987
00988
00989
00990
00991
00992
00993 double kunmat::water_vapour_permeability( double tk,long ipp)
00994 {
00995 double P, Rv, Pa;
00996 double da, dp, mi;
00997
00998
00999 moist = Tm->ip[ipp].eqother[0];
01000
01001
01002 CorD(4,kd,moist,mi,a2,a3);
01003
01004 Rv = R/M;
01005
01006
01007 Pa = 101325;
01008 P = 101325;
01009
01010 da = (2.306e-5 * Pa)/(Rv * tk * P)*pow((tk/273.15),1.81);
01011
01012 if (mi==0.0)
01013 dp=0.0;
01014 else
01015 dp = da/mi;
01016
01017 return (dp);
01018 }
01019
01020
01021
01022
01023
01024
01025 double kunmat::latent_heat_of_evaporation_of_water(double tk)
01026 {
01027 return (2.5008e6)*pow((273.15/tk),(0.167+tk*3.67e-4));
01028 }
01029
01030
01031
01032
01033
01034
01035
01036
01037
01038
01039
01040
01041
01042 double kunmat::saturated_water_vapor_pressure(double tk)
01043 {
01044 return (exp(23.5771 - 4042.9/(tk - 37.58)));
01045 }
01046
01047
01048
01049
01050
01051
01052
01053
01054
01055
01056
01057 double kunmat::derivative_of_the_enthalpy_density(long ipp)
01058 {
01059 double rom, c, moist;
01060
01061
01062 moist = Tm->ip[ipp].eqother[0];
01063
01064
01065 CorD(2,kd,0,rom,a2,a3);
01066
01067
01068 CorD(9,kd,moist,c,a2,a3);
01069
01070 rho_m = rom;
01071
01072 return (rom * c);
01073 }
01074
01075
01076
01077
01078
01079
01080
01081 double kunmat::hygroscopic_moisture(double rh, double tk)
01082 {
01083
01084 double mhmc, mhrh,a3;
01085 CorD(6,kd,rh,mhmc,mhrh,a3);
01086
01087 return (mhmc/(1-sqrt(1-mhrh)));
01088
01089 }
01090
01091
01092
01093
01094
01095
01096
01097
01098 double kunmat::derivative_of_the_sorption_isotherm_root(double rh, double tk)
01099 {
01100
01101 double mhv;
01102 if (rh >1) rh = 1;
01103 if (rh < 0) rh = 0;
01104
01105 mhv = hygroscopic_moisture (rh,tk);
01106
01107 return (1000 * mhv/(2*sqrt(1-rh)));
01108
01109 }
01110
01111
01112
01113
01114
01115
01116
01117 double kunmat::derivative_of_retention_curve_root(double rh, double tk, long ipp)
01118 {
01119
01120 double smc, mhmc, mhrh;
01121
01122 smc = Tm->ip[ipp].eqother[2];
01123 CorD(6,kd,rh,mhmc,mhrh,a3);
01124
01125
01126 return (1000 * (smc - mhmc)/(1 - mhrh));
01127 }
01128
01129
01130
01131
01132
01133
01134
01135
01136
01137
01138 double kunmat::sorption_izotherm_derivative_kk_hansen(double rh, double tk)
01139 {
01140 double a, b, n;
01141
01142 CorD(6,kd,rh,a,b,n);
01143
01144
01145 if (rh >1) rh = 1;
01146 if (rh < 0) rh = 0;
01147
01148
01149 return (1000 * (a/(b*rh*n))*pow((1-(log(rh))/b),(-1-(1/n))));
01150
01151 }
01152
01153
01154
01155
01156
01157
01158
01159
01160
01161
01162
01163
01164
01165
01166
01167 double kunmat::transmission_transcoeff(double trc,long ri,long ci,long nn,long bc,long ipp)
01168 {
01169 long k;
01170 double new_trc,h,t;
01171 new_trc = 0.0;
01172
01173 k=Gtt->give_dof(nn,0);
01174 if (k>0) {h = Lsrst->lhs[k-1]+Lsrst->lhsi[k-1];}
01175 if (k==0) {h = 0.0;}
01176 if (k<0) {h = Tb->lc[0].pv[0-k-1].getval();}
01177 k=Gtt->give_dof(nn,1);
01178 if (k>0) {t = Lsrst->lhs[k-1]+Lsrst->lhsi[k-1];}
01179 if (k==0) {t = 0.0;}
01180 if (k<0) {t = Tb->lc[0].pv[0-k-1].getval();}
01181
01182 if((ri == 0) && (ci == 0))
01183 new_trc = get_transmission_transcoeff_hh(h,t,bc,ipp);
01184 if((ri == 0) && (ci == 1))
01185 new_trc = 0.0;
01186
01187 if((ri == 1) && (ci == 0))
01188 new_trc = 0.0;
01189 if((ri == 1) && (ci == 1))
01190 new_trc = get_transmission_transcoeff_tt(h,t,bc,ipp);
01191
01192 new_trc = new_trc*trc;
01193
01194 return (new_trc);
01195 }
01196
01197
01198
01199
01200
01201
01202
01203
01204
01205
01206
01207
01208
01209
01210 double kunmat::transmission_nodval(double nodval,double trc2,long ri,long ci,long nn,long bc,long ipp)
01211 {
01212 long k;
01213 double new_nodval,h,t;
01214 new_nodval = 0.0;
01215
01216 k=Gtt->give_dof(nn,0);
01217 if (k>0) {h = Lsrst->lhs[k-1]+Lsrst->lhsi[k-1];}
01218 if (k==0) {h = 0.0;}
01219 if (k<0) {h = Tb->lc[0].pv[0-k-1].getval();}
01220 k=Gtt->give_dof(nn,1);
01221 if (k>0) {t = Lsrst->lhs[k-1]+Lsrst->lhsi[k-1];}
01222 if (k==0) {t = 0.0;}
01223 if (k<0) {t = Tb->lc[0].pv[0-k-1].getval();}
01224
01225 if((ri == 0) && (ci == 0))
01226 new_nodval = get_transmission_nodval_hh(nodval,h,t,bc,ipp);
01227 if((ri == 0) && (ci == 1))
01228 new_nodval = 0.0;
01229
01230 if((ri == 1) && (ci == 0))
01231 new_nodval = 0.0;
01232 if((ri == 1) && (ci == 1))
01233 new_nodval = get_transmission_nodval_tt(nodval,h,t,bc,ipp);
01234
01235 return (new_nodval);
01236 }
01237
01238
01239
01240
01241
01242
01243
01244
01245
01246
01247
01248
01249
01250
01251
01252 double kunmat::transmission_flux(double nodval,double trc2,long ri,long ci,long nn,long bc,long ipp)
01253 {
01254 long k;
01255 double flux,h,t;
01256 flux = 0.0;
01257
01258 k=Gtt->give_dof(nn,0);
01259 if (k>0) {h = Lsrst->lhs[k-1]+Lsrst->lhsi[k-1];}
01260 if (k==0) {h = 0.0;}
01261 if (k<0) {h = Tb->lc[0].pv[0-k-1].getval();}
01262 k=Gtt->give_dof(nn,1);
01263 if (k>0) {t = Lsrst->lhs[k-1]+Lsrst->lhsi[k-1];}
01264 if (k==0) {t = 0.0;}
01265 if (k<0) {t = Tb->lc[0].pv[0-k-1].getval();}
01266
01267 if((ri == 0) && (ci == 0))
01268 flux = get_transmission_flux_hh(nodval,h,t,bc,ipp);
01269 if((ri == 0) && (ci == 1))
01270 flux = 0.0;
01271
01272 if((ri == 1) && (ci == 0))
01273 flux = 0.0;
01274 if((ri == 1) && (ci == 1))
01275 flux = get_transmission_flux_tt(nodval,h,t,bc,ipp);
01276
01277 return (flux);
01278 }
01279
01280
01281
01282
01283
01284
01285
01286
01287
01288 double kunmat::get_transmission_transcoeff_hh(double rh,double t,long bc,long ipp)
01289 {
01290 double trc,pgws;
01291
01292 switch (bc){
01293 case 30:{
01294 pgws = saturated_water_vapor_pressure(t);
01295 trc = pgws;
01296 break;
01297 }
01298 case 31:{
01299 pgws = saturated_water_vapor_pressure(t);
01300 trc = pgws;
01301 break;
01302 }
01303 case 32:{
01304 trc = 0.0;
01305 break;
01306 }
01307 default:{
01308 print_err("no real boundary condition is prescribed",__FILE__,__LINE__,__func__);
01309 exit(0);
01310 }
01311 }
01312 return(trc);
01313 }
01314
01315
01316
01317
01318
01319
01320
01321
01322
01323 double kunmat::get_transmission_nodval_hh(double bv,double rh,double t,long bc,long ipp)
01324 {
01325 double nodval,pgws, pgwspred;
01326
01327 switch (bc){
01328 case 30:{
01329
01330 pgws = saturated_water_vapor_pressure(t);
01331 pgwspred = saturated_water_vapor_pressure(t);
01332
01333 nodval = bv*pgwspred;
01334
01335 break;
01336 }
01337 case 31:{
01338 nodval = bv;
01339 break;
01340 }
01341 case 32:{
01342 pgws = saturated_water_vapor_pressure(t);
01343
01344 nodval = pgws*rh;
01345 bv = pgws*bv;
01346 nodval = bv - nodval;
01347
01348 break;
01349 }
01350 default:{
01351 print_err("no real boundary condition is prescribed",__FILE__,__LINE__,__func__);
01352 exit(0);
01353 }
01354 }
01355 return(nodval);
01356 }
01357
01358
01359
01360
01361
01362
01363
01364
01365
01366
01367 double kunmat::get_transmission_flux_hh(double bv,double rh,double t,long bc,long ipp)
01368 {
01369 double flux,pgws;
01370
01371 switch (bc){
01372 case 30:{
01373 pgws = saturated_water_vapor_pressure( t);
01374
01375 flux = pgws*rh;
01376 bv = pgws*bv;
01377 flux = bv - flux;
01378
01379 break;
01380 }
01381 case 31:{
01382 pgws = saturated_water_vapor_pressure(t);
01383 flux = pgws*rh;
01384
01385 flux = bv - flux;
01386 break;
01387 }
01388 case 32:{
01389 pgws = saturated_water_vapor_pressure(t);
01390
01391 flux = pgws*rh;
01392 bv = pgws*bv;
01393 flux = bv - flux;
01394
01395 break;
01396 }
01397 default:{
01398 print_err("no real boundary condition is prescribed",__FILE__,__LINE__,__func__);
01399 exit(0);
01400 }
01401 }
01402 return(flux);
01403 }
01404
01405
01406
01407
01408
01409
01410
01411
01412
01413 double kunmat::get_transmission_transcoeff_tt(double h,double t,long bc,long ipp)
01414 {
01415 double trc;
01416
01417 switch (bc){
01418 case 30:{
01419 trc = 1.0;
01420 break;
01421 }
01422 default:{
01423 print_err("no real boundary condition is prescribed",__FILE__,__LINE__,__func__);
01424 exit(0);
01425 }
01426 }
01427
01428 return(trc);
01429 }
01430
01431
01432
01433
01434
01435
01436
01437
01438
01439 double kunmat::get_transmission_nodval_tt(double bv,double h,double t,long bc,long ipp)
01440 {
01441 double nodval;
01442
01443 switch (bc){
01444 case 30:{
01445 nodval = bv;
01446 break;
01447 }
01448 default:{
01449 print_err("no real boundary condition is prescribed",__FILE__,__LINE__,__func__);
01450 exit(0);
01451 }
01452 }
01453 return(nodval);
01454 }
01455
01456
01457
01458
01459
01460
01461
01462
01463
01464
01465 double kunmat::get_transmission_flux_tt(double bv,double h,double t,long bc,long ipp)
01466 {
01467 double flux;
01468
01469 switch (bc){
01470 case 30:{
01471 flux = (bv - t);
01472 break;
01473 }
01474 default:{
01475 print_err("no real boundary condition is prescribed",__FILE__,__LINE__,__func__);
01476 exit(0);
01477 }
01478 }
01479 return(flux);
01480 }
01481
01482
01483
01484
01485
01486
01487
01488
01489
01490
01491
01492 double kunmat::get_othervalue(long compother,double rh,double t, long ipp)
01493 {
01494 double other;
01495
01496 switch (compother){
01497 case 0:{
01498 other = rh;
01499 break;
01500 }
01501 case 1:{
01502 other = t;
01503 break;
01504 }
01505 case 2:{
01506 double w, a3;
01507 sorption_izoterms_giva_data(0,rh,t,w,a3, ipp);
01508 other = w;
01509 break;
01510 }
01511 case 3:{
01512 other = rh * saturated_water_vapor_pressure (t);
01513 break;
01514 }
01515 case 4:{
01516 double w, a3;
01517 sorption_izoterms_giva_data(0,rh,t,w,a3, ipp);
01518 other = w*rho_w/rho_m;
01519 break;
01520 }
01521 case 5:{
01522 other = rh/saturated_water_vapor_pressure (t);
01523 break;
01524 }
01525 case 6:{
01526 other = 0.0;
01527 break;
01528 }
01529 case 7:{
01530 other = 0.0;
01531 break;
01532 }
01533 case 8:{
01534 other = 0.0;
01535 break;
01536 }
01537 case 9:{
01538 other = 0.0;
01539 break;
01540 }
01541 case 10:{
01542 other = 0.0;
01543 break;
01544 }
01545 default:{
01546 print_err("unknown type of component is required",__FILE__,__LINE__,__func__);
01547 }
01548 }
01549 return (other);
01550
01551 }
01552
01553
01554
01555
01556
01557
01558 void kunmat::print_othervalue_name(FILE *out,long compother)
01559 {
01560 switch (compother){
01561 case 0:{
01562 fprintf (out,"Relative humidity () ");
01563 break;
01564 }
01565 case 1:{
01566 fprintf (out,"Temperature (K) ");
01567 break;
01568 }
01569 case 2:{
01570 fprintf (out,"Moisture content w (m3/m3) ");
01571 break;
01572 }
01573 case 3:{
01574 fprintf (out,"Water vapour pressure (Pa) ");
01575 break;
01576 }
01577 case 4:{
01578 fprintf (out,"Moisture content u (kg/kg) ");
01579 break;
01580 }
01581 default:{
01582 print_err("unknown type of component is required",__FILE__,__LINE__,__func__);
01583 }
01584 }
01585 }
01586
01587
01588 void kunmat::give_data(double rh,double Mhmc, double Smc, double Mhrh, double & moistakt)
01589 {
01590 double hvezdicka;
01591 hvezdicka = Mhmc/(1-sqrt(1-Mhrh));
01592
01593 if (rh < Mhrh)
01594 {
01595 moistakt = (1-sqrt(1-rh)) * hvezdicka ;
01596 }
01597 else {
01598 if (rh > 1)
01599 {
01600 moistakt = Smc;
01601 }
01602 else
01603 {
01604 moistakt = Mhmc + (rh - Mhrh )/(1-Mhrh )* (Smc - Mhmc );
01605 }
01606 }
01607
01608 }
01609
01610
01611
01612
01613
01614
01615
01616
01617
01618
01619 double kunmat::kapa_exp(double a, double b, double rh)
01620 {
01621 return (a * exp(b * rh));
01622 }
01623
01624
01625
01626
01627
01628
01629
01630
01631
01632
01633
01634
01635
01636
01637
01638
01639
01640
01641
01642
01643
01644
01645
01646
01647
01648
01649
01650
01651
01652
01653
01654
01655
01656 void kunmat::CorD(long charid,long &kvyhl, double x, double &y, double &z, double &z2)
01657 {
01658
01659 switch (MatChar[charid]){
01660 case 0:{
01661 kvyhl = 0;
01662 break;
01663 }
01664 case 1:{
01665
01666
01667 y = data[charid]->getval(0.0);
01668 kvyhl = 1;
01669 break;
01670 }
01671 case 2:{
01672
01673
01674 y = data[charid]->getval(x);
01675 kvyhl = 2;
01676
01677
01678 long nr=data[charid]->tabf[0].asize;
01679
01680
01681 z2 = data[charid]->tabf[0].y[nr-1];
01682
01683 break;
01684 }
01685 case 30:{
01686 y = MatFunce[charid][0];
01687 z = MatFunce[charid][1];
01688 kvyhl = 30;
01689 break;
01690 }
01691 case 31:{
01692 y = MatFunce[charid][0];
01693 z = MatFunce[charid][1];
01694 if (charid == 6)
01695 {
01696 z2 = MatFunce[charid][2];
01697 }
01698 kvyhl = 31;
01699 break;
01700 }
01701 case 32:{
01702 y = MatFunce[charid][0];
01703 z = MatFunce[charid][1];
01704 kvyhl = 32;
01705 break;
01706 }
01707 case 33:{
01708 kvyhl = 33;
01709
01710 break;
01711 }
01712 default:{
01713 print_err("unknown definition of material parameter is required",__FILE__,__LINE__,__func__);
01714 }
01715 }
01716 }
01717
01718 void kunmat::sorption_izoterms_giva_data(long kod,double rh, double tk, double & moist, double & dmoistdrh, long ipp)
01719 {
01720
01721
01722
01723 long s;
01724 double hmrh, mhmc, smc, u, a, n;
01725
01726 if (rh < 0) rh = 0;
01727
01728 CorD(6,s,0,a1,a2,a3);
01729
01730 switch (s){
01731 case 2:{
01732 if(rh>0.976)
01733 {
01734
01735 CorD(7,kd,rh,smc,a2,a3);
01736 CorD(6,kd,rh,mhmc,a2,mhmc);
01737 give_data(rh,mhmc, smc, 0.976,moist);
01738 dmoistdrh =1000 * (smc - mhmc)/(1 - 0.976);
01739 }
01740 else
01741 {
01742 CorD(6,s,rh,moist,dmoistdrh,a3);
01743 dmoistdrh =1000* derivative_sorption_izoterm_data(rh);
01744 }
01745 break;
01746 }
01747 case 30:{
01748 CorD(6,s,rh,a1,hmrh,a3);
01749 CorD(7,s,rh,smc,a2,a3);
01750 CorD(6,s,rh,mhmc,hmrh,a3);
01751 give_data(rh,mhmc, smc, hmrh,moist);
01752 if (rh < hmrh)
01753 {
01754 dmoistdrh = derivative_of_the_sorption_isotherm_root (rh, tk);
01755 }
01756 else {
01757 dmoistdrh = derivative_of_retention_curve_root(rh, tk, ipp);
01758 }
01759
01760 break;
01761 }
01762 case 31:{
01763 switch (kod){
01764 case 0:{
01765 CorD(6,kod,0,u,a,n);
01766 moist = si_kk_hansen(rh, 0,u,a,n);
01767 dmoistdrh = sorption_izotherm_derivative_kk_hansen(rh, tk);
01768 break;
01769 }
01770 case 1:{
01771 dmoistdrh = sorption_izotherm_derivative_kk_hansen(rh, tk);
01772 break;
01773 }
01774 default:{
01775 print_err("unknown definition of Sorption isotherm HANSEN is required",__FILE__,__LINE__,__func__);
01776 }
01777 }
01778 break ;
01779 }
01780 default:{
01781 print_err("unknown definition of Sorption isotherm is required",__FILE__,__LINE__,__func__);
01782 }
01783 }
01784
01785
01786 }
01787
01788 double kunmat::si_kk_hansen(double rh, double tk,double u, double a, double n)
01789 {
01790 return (u*pow((1-(log(rh))/a),(-1/n)));
01791 }
01792
01793
01794
01795
01796
01797
01798
01799
01800 double kunmat::derivative_sorption_izoterm_data(double x1)
01801 {
01802 long charid;
01803 double derbi;
01804
01805
01806 charid=6;
01807 derbi = data[charid]->getderiv (x1);
01808
01809 return derbi;
01810 }
01811
01812
01813
01814
01815
01816
01817
01818
01819
01820 void kunmat::aux_values (long ipp,double *in,double *inp, double *ineq,double *out)
01821 {
01822
01823
01824 switch (kunzeltype) {
01825 case 0:{
01826 double rh,tk;
01827 double w,dwdf, kapaA, smc;
01828
01829
01830 rh = in[0];
01831
01832 tk = in[1];
01833
01834 smc = data[7]->getval(0.0);
01835
01836 sorption_izoterms_giva_data(0,rh,tk,w,dwdf, ipp);
01837 if (w>smc) w = smc;
01838 if(w<0) w = 0.0;
01839
01840 kapa_values(MatChar[5],ipp,w, kapaA);
01841
01842
01843 out[3] = kapaA;
01844
01845 out[0] = w;
01846
01847 out[1] = dwdf;
01848
01849 out[2] = smc;
01850 break;
01851 }
01852 case 1:{
01853 double pv,tk, pvs;
01854 double w,dwdf, kapaA, smc, rh;
01855
01856 pv = Tm->ip[ipp].av[0];
01857
01858 tk = Tm->ip[ipp].av[1];
01859
01860 pvs = saturated_water_vapor_pressure(tk);
01861
01862 rh = pv/pvs;
01863
01864
01865 smc = data[7]->getval(0.0);
01866
01867 sorption_izoterms_giva_data(0,rh,tk,w,dwdf, ipp);
01868 if (w>smc) w = smc;
01869 if(w<0) w = 0.0;
01870
01871 kapa_values(MatChar[5],ipp,w, kapaA);
01872
01873
01874 out[3] = kapaA;
01875
01876 out[0] = w;
01877
01878 out[1] = -10000.0;
01879
01880 out[2] = smc;
01881
01882 break;
01883 }
01884 default:
01885 ;
01886 }
01887 }
01888
01889
01890
01891
01892
01893
01894
01895
01896
01897 void kunmat::save_values (long ipp,double *out)
01898 {
01899 Tm->ip[ipp].eqother[0]=out[0];
01900 Tm->ip[ipp].eqother[1]=out[1];
01901 Tm->ip[ipp].eqother[2]=out[2];
01902 Tm->ip[ipp].eqother[3]=out[3];
01903 }
01904
01905
01906
01907
01908
01909
01910
01911
01912
01913
01914 void kunmat::give_values (long ipp,double *av,double *eq)
01915 {
01916 av[0] = Tm->ip[ipp].av[0];
01917 av[1] = Tm->ip[ipp].av[1];
01918
01919 eq[0] = Tm->ip[ipp].eqother[0];
01920 eq[1] = Tm->ip[ipp].eqother[1];
01921 eq[2] = Tm->ip[ipp].eqother[2];
01922 eq[3] = Tm->ip[ipp].eqother[3];
01923 }
01924
01925
01926
01927
01928
01929
01930 void kunmat::kapa_values (long kod, long ipp,double x1, double &kapa)
01931 {
01932 double kapa1, kapa2, kapa3;
01933 long kd1;
01934
01935 switch (kod){
01936 case 0: kapa = 0.0;
01937 break;
01938 case 1: CorD(5,kod,x1,kapa,kapa1,kapa2);
01939 break;
01940 case 2: CorD(5,kod,x1,kapa,kapa1,kapa2);
01941 break;
01942 case 30:{
01943 CorD(5,kod,x1,kapa1,kapa2,kapa3);
01944 kapa = kapa_exp(kapa1,kapa2,x1);
01945 break;
01946 }
01947 default:{
01948 print_err("unknown kod is required",__FILE__,__LINE__,__func__);
01949 }
01950 }
01951 }
01952
01953 void kunmat::initvalues (long ipp,long ido)
01954 {
01955 double x1,x2;
01956
01957 x1 = Tm->ip[ipp].av[0];
01958 x2 = Tm->ip[ipp].av[1];
01959
01960 }
01961
01962