00001 #include "globalt.h"
00002 #include "millymat.h"
00003 #include "kunzel.h"
00004 #include "math.h"
00005
00006
00007
00008
00009
00010
00011 static double T0 = 293.15;
00012 static double M = 0.018;
00013 static double R = 8.314;
00014 static double rl = 1000;
00015 static double p = 101325;
00016 static double g = 9.81;
00017
00018
00019 static double Rvapd = 462.0;
00020 static double cwatd = 4180;
00021 static double cvapd = 1617;
00022
00023
00024 static kunmat *kun;
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044 millymat::millymat ()
00045 {}
00046
00047 millymat::~millymat ()
00048 {}
00049
00050 void millymat::values_correction (vector &nv)
00051 {
00052
00053 nv[0] = nv[0];
00054 nv[1] = nv[1];
00055
00056
00057 }
00058
00059
00060
00061
00062
00063
00064
00065
00066
00067
00068
00069 void millymat::matcond (matrix &d,long ri,long ci,long ipp)
00070 {
00071 long n;
00072 n = d.n;
00073
00074 switch (n){
00075 case 1:{
00076 matcond1d (d,ri,ci,ipp);
00077 break;
00078 }
00079 case 2:{
00080 matcond2d (d,ri,ci,ipp);
00081 break;
00082 }
00083 case 3:{
00084 matcond3d (d,ri,ci,ipp);
00085 break;
00086 }
00087 default:{
00088 fprintf (stderr,"\n unknown number of components of stress tensor is required");
00089 fprintf (stderr,"\n in function (%s, line %d).\n",__FILE__,__LINE__);
00090 }
00091 }
00092 }
00093
00094
00095
00096
00097
00098
00099
00100
00101
00102
00103 void millymat::matcond1d (matrix &d,long ri,long ci,long ipp)
00104 {
00105 double k;
00106 double x1,x2;
00107 k = 0.0;
00108
00109 x1 = Tm->ip[ipp].av[0];
00110 x2 = Tm->ip[ipp].av[1];
00111
00112 if((ri == 0) && (ci == 0))
00113 k = k11 (x1,x2, ipp);
00114 if((ri == 0) && (ci == 1))
00115 k = k12 (x1,x2, ipp);
00116
00117 if((ri == 1) && (ci == 0))
00118 k = k21 (x1,x2, ipp);
00119 if((ri == 1) && (ci == 1))
00120 k = k22 (x1,x2, ipp);
00121
00122 d[0][0] = k;
00123 }
00124
00125
00126
00127
00128
00129
00130
00131
00132
00133 void millymat::matcond2d (matrix &d,long ri,long ci,long ipp)
00134 {
00135 double k;
00136 double x1,x2;
00137 k = 0.0;
00138
00139 x1 = Tm->ip[ipp].av[0];
00140 x2 = Tm->ip[ipp].av[1];
00141
00142 if((ri == 0) && (ci == 0))
00143 k = k11 (x1,x2, ipp);
00144 if((ri == 0) && (ci == 1))
00145 k = k12 (x1,x2, ipp);
00146
00147 if((ri == 1) && (ci == 0))
00148 k = k21 (x1,x2, ipp);
00149 if((ri == 1) && (ci == 1))
00150 k = k22 (x1,x2, ipp);
00151
00152 fillm(0.0,d);
00153
00154 d[0][0] = k; d[0][1] = 0.0;
00155 d[1][0] = 0.0; d[1][1] = k;
00156 }
00157
00158
00159
00160
00161
00162
00163
00164
00165
00166
00167 void millymat::matcond3d (matrix &d,long ri,long ci,long ipp)
00168 {
00169 double k;
00170 double x1,x2;
00171 k = 0.0;
00172
00173 x1 = Tm->ip[ipp].av[0];
00174 x2 = Tm->ip[ipp].av[1];
00175
00176 if((ri == 0) && (ci == 0))
00177 k = k11 (x1,x2, ipp);
00178 if((ri == 0) && (ci == 1))
00179 k = k12 (x1,x2, ipp);
00180
00181 if((ri == 1) && (ci == 0))
00182 k = k21 (x1,x2, ipp);
00183 if((ri == 1) && (ci == 1))
00184 k = k22 (x1,x2, ipp);
00185
00186 fillm(0.0,d);
00187
00188 d[0][0]=k; d[0][1]=0.0; d[0][2]=0.0;
00189 d[1][0]=0.0; d[1][1]=k; d[1][2]=0.0;
00190 d[2][0]=0.0; d[2][1]=0.0; d[2][2]=k;
00191 }
00192
00193
00194
00195
00196
00197
00198
00199
00200
00201
00202 void millymat::matcap (double &c,long ri,long ci,long ipp)
00203 {
00204 double x1,x2;
00205 c = 0.0;
00206
00207 x1 = Tm->ip[ipp].av[0];
00208 x2 = Tm->ip[ipp].av[1];
00209
00210 if((ri == 0) && (ci == 0))
00211 c = c11 (x1,x2, ipp);
00212 if((ri == 0) && (ci == 1))
00213 c = c12 (x1,x2, ipp);
00214
00215 if((ri == 1) && (ci == 0))
00216 c = c21 (x1,x2, ipp);
00217 if((ri == 1) && (ci == 1))
00218 c = c22 (x1,x2, ipp);
00219
00220 }
00221
00222
00223
00224
00225
00226
00227
00228
00229 void millymat::read (XFILE *in)
00230 {
00231 int i, j;
00232
00233 for (i = 2; i<19; i++)
00234 {
00235 xfscanf(in, "%d",&MatChar[i]);
00236 }
00237 xfscanf(in, "%d",&MatChar[19]);
00238 for (i = 2;i <= 19; i++)
00239 {
00240
00241 switch (MatChar[i]){
00242 case 0:{
00243 break;
00244 }
00245 case 1:{
00246 xfscanf(in, "%lf",&MatConst[i]);
00247
00248 break;
00249 }
00250 case 2:{
00251 if( MatChar[i-1] ==2)
00252 {
00253 xfscanf(in, "%lf %lf",&MatData[i][0][0],&MatData[i][1][0]);
00254 }
00255 else
00256 {
00257 xfscanf(in, "%lf %lf",&MatData[i][0][0],&MatData[i][1][0]);
00258 }
00259
00260 switch (int(MatData[i][1][0])){
00261 case 2:{
00262 for (j=1; j<=MatData[i][0][0];j++)
00263 {
00264 xfscanf(in, "%lf %lf",&MatData[i][0][j],&MatData[i][1][j] );
00265 }
00266 break;
00267 }
00268 case 3:{
00269 for (j=1; j<=MatData[i][0][0];j++)
00270 {
00271 xfscanf(in, "%lf %lf %lf",&MatData[i][0][j],&MatData[i][1][j],&MatData[i][2][j] );
00272 }
00273 break;
00274 }
00275 default:{
00276 fprintf (stderr,"\n\n unknown definition of Material KUNZEL-DATA is required");
00277 fprintf (stderr,"\n in function read (file %s, line %d)\n",__FILE__,__LINE__);
00278 }
00279 }
00280 xfscanf(in, "%d",&MatData[i][0][0]);
00281
00282 break;
00283 }
00284 case 3:{
00285
00286 break;
00287 }
00288 case 30:{
00289 if( MatChar[i-1] ==1)
00290 {
00291 xfscanf(in, "%lf %lf",&MatFunce[i][0],&MatFunce[i][1]);
00292 }
00293 else
00294 {
00295 xfscanf(in, "%lf %lf",&MatFunce[i][0],&MatFunce[i][1]);
00296 }
00297 break;
00298 }
00299 case 31:{
00300 if( MatChar[i-1] ==1)
00301 {
00302 if (i == 6)
00303 {
00304 xfscanf(in, "%lf %lf %lf",&MatFunce[i][0],&MatFunce[i][1], &MatFunce[i][2]);
00305 }
00306 else
00307 {
00308 xfscanf(in, "%lf %lf",&MatFunce[i][0],&MatFunce[i][1]);
00309 }
00310 }
00311 else
00312 {
00313 if (i == 6)
00314 {
00315 xfscanf(in, "%lf %lf %lf",&MatFunce[i][0],&MatFunce[i][1], &MatFunce[i][2]);
00316 }
00317 else
00318 {
00319 xfscanf(in, "%lf %lf",&MatFunce[i][0],&MatFunce[i][1]);
00320 }
00321 }
00322
00323 break;
00324 }
00325 case 32:{
00326 if( MatChar[i-1] ==1)
00327 {
00328 xfscanf(in, "%lf %lf",&MatFunce[i][0],&MatFunce[i][1]);
00329 }
00330 else
00331 {
00332 xfscanf(in, "%lf %lf",&MatFunce[i][0],&MatFunce[i][1]);
00333 }
00334 break;
00335 }
00336 default:{
00337 fprintf (stderr,"\n\n unknown definition of Material KUNZEL is required");
00338 fprintf (stderr,"\n in function read (file %s, line %d)\n",__FILE__,__LINE__);
00339 }
00340 }
00341 }
00342
00343
00344
00345 }
00346
00347
00348
00349
00350
00351
00352 double millymat::k11(double x1,double x2,long ipp)
00353 {
00354 double k11, rs, Kl, Kv;
00355
00356 CorD(2,kd,0,0,x1,rs,a2,a3);
00357
00358 Kl = Tm->ip[ipp].eqother[5];
00359 Kv =coeff_Kv(x1,x2,ipp);
00360
00361
00362 k11 = rl*(Kl+Kv);
00363 return (k11);
00364 }
00365
00366 double millymat::k12(double x1,double x2,long ipp)
00367 {
00368 double k12, rs,Dtv;
00369
00370 CorD(2,kd,0,0,x1,rs,a2,a3);
00371 Dtv =coeff_Dtv(x1,x2,ipp);
00372
00373
00374 k12 = rl*Dtv;
00375
00376 return (k12);
00377 }
00378
00379
00380 double millymat::k21(double x1,double x2,long ipp)
00381 {
00382 double k21, rs,Kl, Kv, hl, hv;
00383
00384 CorD(2,kd,0,0,x1,rs,a2,a3);
00385 Kl = Tm->ip[ipp].eqother[5];
00386 Kv =coeff_Kv(x1,x2,ipp);
00387 hl =specific_enthalpy_of_liquid_phase(x1,x2);
00388 hv = specific_enthalpy_of_gaseous_phase(x1,x2);
00389
00390
00391 k21 = (-hl*rl*Kl-hv*rl*Kv);
00392
00393 return (k21);
00394 }
00395
00396 double millymat::k22(double x1,double x2,long ipp)
00397 {
00398 double k22, rs, lambda, Dtv, hv;
00399 CorD(2,kd,0,0,x1,rs,a2,a3);
00400 lambda = Tm->ip[ipp].eqother[6];
00401 hv = specific_enthalpy_of_gaseous_phase(x1,x2);
00402 Dtv =coeff_Dtv(x1,x2,ipp);
00403
00404
00405 k22 = (lambda-Dtv*hv*rl);
00406
00407 fprintf(Outt,"k22 je %3.6e pro rh= %lf a t= %lf \n",k22,x1,x2);
00408 return (k22);
00409 }
00410
00411 double millymat::c11(double x1,double x2,long ipp)
00412 {
00413 double c11, rs,A, B, duldpsi;
00414 CorD(2,kd,0,0,x1,rs,a2,a3);
00415 A =coeef_zaporneA(x1,x2, ipp);
00416 B = coeef_B(x1,x2, ipp);
00417 duldpsi = Tm->ip[ipp].eqother[8];
00418
00419 c11 = (1+A)*duldpsi + B;
00420 c11 = c11 *rs;
00421 return(c11);
00422 }
00423
00424 double millymat::c12(double x1,double x2,long ipp)
00425 {
00426 double c12,rs, A, C, duldt;
00427
00428 CorD(2,kd,0,0,x1,rs,a2,a3);
00429 A =coeef_zaporneA(x1,x2, ipp);
00430 C = coeef_C(x1,x2, ipp);
00431 duldt = Tm->ip[ipp].eqother[4];
00432
00433 c12 = (1+A)*duldt + C;
00434 c12 = c12 *rs;
00435
00436 return(c12);
00437 }
00438
00439
00440 double millymat::c21(double x1,double x2,long ipp)
00441 {
00442 double c21, A, B, duldpsi, hl, hv, rs;
00443 A =coeef_zaporneA(x1,x2, ipp);
00444 B = coeef_B(x1,x2, ipp);
00445 duldpsi = Tm->ip[ipp].eqother[8];
00446 hv = specific_enthalpy_of_gaseous_phase(x1,x2);
00447 hl = specific_enthalpy_of_liquid_phase(x1,x2);
00448
00449 c21 = (hl+hv*A)*duldpsi+B*hv;
00450
00451 CorD(2,kd,0,0,x1,rs,a2,a3);
00452 c21 = c21 * rs;
00453 return(c21);
00454 }
00455
00456 double millymat::c22(double x1,double x2,long ipp)
00457 {
00458 double c22, A, C, duldt, hl, hv, cs, ul, uv, rs;
00459
00460 A =coeef_zaporneA(x1,x2, ipp);
00461 C = coeef_C(x1,x2, ipp);
00462 duldt = Tm->ip[ipp].eqother[4];
00463 ul = Tm->ip[ipp].eqother[0];
00464 hv = specific_enthalpy_of_gaseous_phase(x1,x2);
00465 hl = specific_enthalpy_of_liquid_phase(x1,x2);
00466 CorD(9,kd,0,0,x1,cs,a2,a3);
00467 uv = gaseous_moisture_content_by_mass(x1,x2,ipp);
00468
00469 c22 = (hl+A*hv)*duldt+(cs+cwatd*ul+cvapd*uv+hv*C);
00470 CorD(2,kd,0,0,x1,rs,a2,a3);
00471 c22 = c22*rs;
00472 fprintf(Outt,"c22 je %3.6e pro rh= %lf a t= %lf \n",c22,x1,x2);
00473 return(c22);
00474 }
00475
00476
00477
00478
00479 void millymat::auxiliarydata (double x1,double x2)
00480 {}
00481
00482
00483
00484
00485
00486
00487
00488
00489
00490
00491
00492
00493
00494
00495
00496
00497
00498 double millymat::transmission_transcoeff(double trc,long ri,long ci,long nn,long bc,long ipp)
00499 {
00500 long k;
00501 double new_trc,x1,x2;
00502 new_trc = 0.0;
00503
00504 k=Gtt->give_dof(nn,0);
00505 if (k>0) {x1 = Lsrst->lhs[k-1]+Lsrst->lhsi[k-1];}
00506 if (k==0) {x1 = 0.0;}
00507 if (k<0) {x1 = Tb->lc[0].pv[0-k-1].getval();}
00508 k=Gtt->give_dof(nn,1);
00509 if (k>0) {x2 = Lsrst->lhs[k-1]+Lsrst->lhsi[k-1];}
00510 if (k==0) {x2 = 0.0;}
00511 if (k<0) {x2 = Tb->lc[0].pv[0-k-1].getval();}
00512
00513
00514 if((ri == 0) && (ci == 0))
00515 new_trc = get_transmission_transcoeff_11(x1,x2,bc,ipp);
00516 if((ri == 0) && (ci == 1))
00517 new_trc = 0.0;
00518
00519
00520 if((ri == 1) && (ci == 0))
00521 new_trc = 0.0;
00522 if((ri == 1) && (ci == 1))
00523 new_trc = 0.0;
00524
00525
00526 new_trc = new_trc*trc;
00527
00528 return (new_trc);
00529 }
00530
00531
00532
00533
00534
00535
00536
00537
00538
00539
00540
00541
00542 double millymat::transmission_nodval(double nodval,long ri,long ci,long nn,long bc,long ipp)
00543 {
00544 long k;
00545 double new_nodval,x1,x2;
00546 new_nodval = 0.0;
00547
00548 k=Gtt->give_dof(nn,0);
00549 if (k>0) {x1 = Lsrst->lhs[k-1]+Lsrst->lhsi[k-1];}
00550 if (k==0) {x1 = 0.0;}
00551 if (k<0) {x1 = Tb->lc[0].pv[0-k-1].getval();}
00552 k=Gtt->give_dof(nn,1);
00553 if (k>0) {x2 = Lsrst->lhs[k-1]+Lsrst->lhsi[k-1];}
00554 if (k==0) {x2 = 0.0;}
00555 if (k<0) {x2 = Tb->lc[0].pv[0-k-1].getval();}
00556
00557
00558 if((ri == 0) && (ci == 0))
00559 new_nodval = get_transmission_nodval_11(nodval,x1,x2,bc,ipp);
00560 if((ri == 0) && (ci == 1))
00561 new_nodval = 0.0;
00562
00563
00564 if((ri == 1) && (ci == 0))
00565 new_nodval = 0.0;
00566 if((ri == 1) && (ci == 1))
00567 new_nodval = 0.0;
00568
00569
00570
00571
00572 return (new_nodval);
00573 }
00574
00575
00576
00577
00578
00579
00580
00581
00582
00583
00584
00585
00586
00587 double millymat::transmission_flux(double nodval,long ri,long ci,long nn,long bc,long ipp)
00588 {
00589 long k;
00590 double flux,x1,x2;
00591 flux = 0.0;
00592
00593 k=Gtt->give_dof(nn,0);
00594 if (k>0) {x1 = Lsrst->lhs[k-1]+Lsrst->lhsi[k-1];}
00595 if (k==0) {x1 = 0.0;}
00596 if (k<0) {x1 = Tb->lc[0].pv[0-k-1].getval();}
00597 k=Gtt->give_dof(nn,1);
00598 if (k>0) {x2 = Lsrst->lhs[k-1]+Lsrst->lhsi[k-1];}
00599 if (k==0) {x2 = 0.0;}
00600 if (k<0) {x2 = Tb->lc[0].pv[0-k-1].getval();}
00601
00602
00603 if((ri == 0) && (ci == 0))
00604 flux = get_transmission_flux_11(nodval,x1,x2,bc,ipp);
00605 if((ri == 0) && (ci == 1))
00606 flux = 0.0;
00607
00608
00609 if((ri == 1) && (ci == 0))
00610 flux = 0.0;
00611 if((ri == 1) && (ci == 1))
00612 flux = 0.0;
00613
00614
00615 return (flux);
00616 }
00617
00618
00619
00620
00621
00622
00623
00624
00625
00626
00627 double millymat::get_transmission_nodval_11(double bv,double x1,double x2,long bc,long ipp)
00628 {
00629 double nodval;
00630
00631 switch (bc){
00632 case 30:{
00633 nodval = bv;
00634 break;
00635 }
00636 default:{
00637 fprintf (stderr,"\n\n No real boundary condition is prescribed (%s, line %d).\n",__FILE__,__LINE__);
00638 exit(0);
00639 }
00640 }
00641
00642 return(nodval);
00643 }
00644
00645
00646
00647
00648
00649
00650
00651
00652
00653 double millymat::get_transmission_transcoeff_11(double x1,double x2,long bc,long ipp)
00654 {
00655 double trc;
00656
00657 switch (bc){
00658 case 30:{
00659 trc=1.0;
00660 break;
00661 }
00662 default:{
00663 fprintf (stderr,"\n\n No real boundary condition is prescribed (%s, line %d).\n",__FILE__,__LINE__);
00664 exit(0);
00665 }
00666 }
00667
00668 return(trc);
00669 }
00670
00671
00672
00673
00674
00675
00676
00677
00678
00679
00680 double millymat::get_transmission_flux_11(double bv,double x1,double x2,long bc,long ipp)
00681 {
00682 double flux;
00683
00684 switch (bc){
00685 case 30:{
00686 flux = (bv - x1);
00687 break;
00688 }
00689 default:{
00690 fprintf (stderr,"\n\n No real boundary condition is prescribed (%s, line %d).\n",__FILE__,__LINE__);
00691 exit(0);
00692 }
00693 }
00694
00695 return(flux);
00696 }
00697
00698
00699
00700
00701
00702
00703
00704
00705 double millymat::get_othervalue(long compother,long ipp, double x1,double x2)
00706 {
00707 double other;
00708
00709 switch (compother){
00710 case 0:{
00711 other = relative_humidity_psi(x1,293.15);
00712 break;
00713 }
00714 case 1:
00715 case 2:
00716 case 3:
00717 case 4:
00718 case 5:
00719 case 6:
00720 case 7:
00721 case 8:
00722 case 9:
00723 case 10:
00724 case 11:
00725 case 12:
00726 break;
00727 default:{
00728
00729 }
00730 }
00731 return (other);
00732 }
00733
00734
00735
00736
00737
00738
00739 void millymat::print_othervalue_name(FILE *out,long compother)
00740 {
00741 switch (compother){
00742 case 0:{
00743 fprintf (out,"First unknown (Units) ");
00744 break;
00745 }
00746 default:{
00747 fprintf (stderr,"\n\n unknown type of component is required in function (%s, line %d).\n",__FILE__,__LINE__);
00748 }
00749 }
00750 }
00751
00752 double millymat::saturation_water_vapour_pressure(double x1, double x2)
00753 {
00754 return (exp(23.5771 - 4042.9/(x2 - 37.58)));
00755 }
00756
00757 double millymat::derivation_saturation_water_vapour_pressure_temperature(double x1, double x2)
00758 {
00759 return ((4042.9 * exp(23.5771-4042.9/(x2 - 37.58)))/((x2-37.58)*(x2-37.58)));
00760 }
00761
00762 double millymat::relative_volume_ration_a(double x1, double x2, double ul)
00763 {
00764 double pi, rs;
00765 CorD(3,kd,0,0,x1,pi,a2,a3);
00766 CorD(2,kd,0,0,x1,rs,a2,a3);
00767
00768
00769
00770 return (pi-rs/rl*ul);
00771 }
00772 double millymat::relative_humidity_pc(double x1, double x2)
00773 {
00774 double pc, fi;
00775
00776
00777
00778 rl = 1000;
00779 fi = exp(-pc*M/(rl*R*x2));
00780
00781 if (fi>1.0) {
00782 fi = 1.0;
00783 }
00784 return (fi);
00785 }
00786 double millymat::relative_humidity_psi(double x1, double x2)
00787 {
00788 double fi;
00789
00790
00791
00792 fi = exp(x1*M*g/(R*x2));
00793 return (fi);
00794 }
00795 double millymat::coeef_zaporneA (double x1, double x2, long ipp)
00796 {
00797 double fi, pvs;
00798 fi = Tm->ip[ipp].eqother[2];
00799 pvs = saturation_water_vapour_pressure(x1,x2);
00800 return (-1/rl*fi*pvs*M/(R*x2));
00801
00802 }
00803
00804 double millymat::coeef_B (double x1, double x2, long ipp)
00805 {
00806 double fi, pvs, a, rs;
00807 fi = Tm->ip[ipp].eqother[2];
00808 pvs = saturation_water_vapour_pressure(x1,x2);
00809 a = Tm->ip[ipp].eqother[7];
00810 CorD(2,kd,0,0,x1,rs,a2,a3);
00811 return (fi*pvs*a/rs*(M/(R*x2))*(M/(R*x2)));
00812
00813 }
00814 double millymat::coeef_C (double x1, double x2, long ipp)
00815 {
00816 double fi, pvs, a, rs, dpvsdt;
00817
00818 fi = Tm->ip[ipp].eqother[2];
00819 pvs = saturation_water_vapour_pressure(x1,x2);
00820 dpvsdt = derivation_saturation_water_vapour_pressure_temperature(x1,x2);
00821 a = Tm->ip[ipp].eqother[7];
00822
00823 CorD(2,kd,0,0,x1,rs,a2,a3);
00824
00825
00826 return (1/rs*M/(R*x2)*a*pvs*fi*(-x1*g*M/(R*x2*x2)+1/pvs*dpvsdt-1/x2));
00827
00828
00829 }
00830 double millymat::diffusion_coefficient_of_water_vapor_in_air(double x1, double x2)
00831 {
00832 double mi;
00833 CorD(4,kd,0,0,x1,mi,a2,a3);
00834 return(2.3e-5 * pow(x2/273.15,1.81));
00835 }
00836
00837 double millymat::partial_water_vapour_pressure_function(double x1, double x2, long ipp)
00838 {
00839 double fi, pvs;
00840
00841 pvs = saturation_water_vapour_pressure(x1,x2);
00842
00843 fi = Tm->ip[ipp].eqother[2];
00844 return(pvs*fi);
00845
00846 }
00847
00848 double millymat::coeff_Kv (double x1, double x2, long ipp)
00849 {
00850 double alf,a, mi, D, pv, fi,pvs;
00851
00852 CorD(4,kd,0,0,x1,mi,a2,a3);
00853 alf =1/mi;
00854
00855 a = Tm->ip[ipp].eqother[7];
00856 D = diffusion_coefficient_of_water_vapor_in_air(x1,x2);
00857 fi = Tm->ip[ipp].eqother[2];
00858 pv = partial_water_vapour_pressure_function(x1,x2,ipp);
00859 pvs = saturation_water_vapour_pressure(x1,x2);
00860
00861 return(1/rl*alf*a*D*p/(p-pv)*(M/(R*x2))*(M/(R*x2))*g*pvs*fi );
00862 }
00863
00864 double millymat::coeff_Dtv (double x1, double x2, long ipp)
00865 {
00866 double alf,a, mi, D, pv, fi,pvs,dpvsdt;
00867
00868 CorD(4,kd,0,0,x1,mi,a2,a3);
00869 alf =1/mi;
00870
00871 a = Tm->ip[ipp].eqother[7];
00872 D = diffusion_coefficient_of_water_vapor_in_air(x1,x2);
00873 fi = Tm->ip[ipp].eqother[2];
00874 pv = partial_water_vapour_pressure_function(x1,x2,ipp);
00875 pvs = saturation_water_vapour_pressure(x1,x2);
00876 dpvsdt = derivation_saturation_water_vapour_pressure_temperature(x1,x2);
00877
00878 return(1/rl*alf*a*D*p/(p-pv)*(M/(R*x2))*fi*(-pvs/x2+dpvsdt-M/(R*x2*x2)*g*x1*pvs));
00879 }
00880 double millymat::gaseous_moisture_content_by_mass (double x1, double x2, long ipp)
00881 {
00882
00883 double pvs, ul, pi, rs;
00884
00885 pvs = saturation_water_vapour_pressure(x1,x2);
00886 CorD(3,kd,0,0,x1,pi,a2,a3);
00887 CorD(2,kd,0,0,x1,rs,a2,a3);
00888 ul = Tm->ip[ipp].eqother[0];
00889
00890 return( (pvs*M/(R*x2))*(pi/rs-ul/rl)*exp(x1*g*M/(R*x2)));
00891
00892 }
00893
00894 double millymat:: specific_enthalpy_of_porous_matrix (double x1, double x2)
00895 {
00896 double cs;
00897 CorD(9,kd,0,0,x1,cs,a2,a3);
00898
00899 return(cs*(x2-T0));
00900 }
00901
00902 double millymat:: specific_enthalpy_of_liquid_phase (double x1, double x2)
00903 {
00904 return(cwatd*(x2-T0));
00905 }
00906 double millymat:: specific_enthalpy_of_gaseous_phase (double x1, double x2)
00907 {
00908 double Lv;
00909 Lv =latent_heat_of_evaporation_of_water(x1,x2);
00910 return(cvapd*(x2-T0)+Lv);
00911 }
00912 double millymat::latent_heat_of_evaporation_of_water(double x1, double x2)
00913 {
00914 return (2.5008e6)*pow((273.15/x2),(0.167+x2*3.67e-4));
00915 }
00916
00917 double millymat:: derivation_water_retention_curve_by_pressure_head (double x1, double x2)
00918 {
00919 double dulpsi;
00920
00921 return (dulpsi);
00922 }
00923
00924 double millymat:: derivation_water_retention_curve_by_temperature (double x1, double x2)
00925 {
00926 double ul, dut;
00927 dut = 0.0;
00928 return (dut);
00929 }
00930 double millymat::give_ul (double x1, double x2, double &fi, double &dmoistdrh )
00931 {
00932 double ul;
00933
00934
00935
00936 fi = relative_humidity_psi(x1,x2);
00937 if (fi>1.0) {
00938 fi = 1.0;
00939 }
00940 if (fi>0.976) {
00941 ul =give_ul_retention_curve(x1,x2,fi, dmoistdrh);
00942
00943 }
00944 else {
00945 ul = give_ul_sorption_izothemrs(x1,x2,fi, dmoistdrh);
00946
00947 }
00948
00949 return(ul);
00950 }
00951 double millymat::give_ul_sorption_izothemrs (double x1, double x2,double rh, double &dmoistdrh )
00952 {
00953 double rs;
00954
00955
00956 int s;
00957 double moist,hmrh, mhmc, smc, u, a, n;
00958
00959 if (rh < 0) rh = 0;
00960
00961 CorD(6,s,0,0,0,a1,a2,a3);
00962
00963 switch (s)
00964 {
00965 case 2: {
00966 CorD(6,s,1,0,rh,moist,a1,a3);
00967 dmoistdrh =derivation_sorption_izoterm_data(rh,x2, 6);
00968 }
00969 break;
00970 case 30: CorD(6,s,1,0,rh,a1,hmrh,a3);
00971 CorD(7,s,0,0,rh,smc,a2,a3);
00972 CorD(6,s,0,0,rh,mhmc,hmrh,a3);
00973
00974 dmoistdrh = DerivativeOfTheSorptionIsotherm (x1,x2,rh,moist);
00975
00976 break;
00977 case 31: CorD(6,s,0,0,0,u,a,n);
00978
00979
00980 break ;
00981 }
00982 CorD(2,kd,0,0,moist,rs,a2,a3);
00983
00984 return (moist*rl/rs);
00985 }
00986
00987 double millymat::give_ul_retention_curve (double x1, double x2, double rh, double &dmoistdrh )
00988 {
00989 double moist, rs, pc;
00990 int s;
00991
00992 pc = capilar_pressure (x1,x2);
00993
00994
00995 CorD(11,s,0,0,0,a1,a2,a3);
00996
00997
00998
00999 switch (s)
01000 {
01001 case 2: {
01002 CorD(11,s,1,0,pc,moist,a1,a3);
01003 dmoistdrh =derivation_sorption_izoterm_data(pc,x2, 11);
01004 }
01005 break;
01006 }
01007
01008 CorD(2,kd,0,0,moist,rs,a2,a3);
01009
01010 return (moist*rl/rs);
01011
01012 }
01013 double millymat::capilar_pressure (double x1, double x2)
01014 {
01015 return(-x1*g*rl);
01016 }
01017
01018
01019 double millymat::ul_to_w (double ul)
01020 {
01021 double rs;
01022 CorD(2,kd,0,0,ul,rs,a2,a3);
01023 return(ul*rs/rl);
01024 }
01025
01026 double millymat::derivation_sorption_izoterm_data(double x1, double x2, int kod)
01027 {
01028
01029 int i = 1;
01030 double derbi;
01031 if ( x1 < 0)
01032 {
01033 x1 = 0;
01034 }
01035
01036 while (x1 > MatData[kod][0][i]){
01037 i++;
01038 }
01039 derbi = (MatData [kod][1][i] - MatData [kod][1][i-1])/(MatData [kod][0][i] - MatData [kod][0][i-1]);
01040 return(derbi);
01041 }
01042
01043 double millymat::DerivativeOfTheSorptionIsotherm(double x1, double x2,double rh,double &moistakt1)
01044 {
01045
01046 double mhmc, mhrh, rs,rh2,moistakt2, ul1, ul2, psi1, psi2;
01047 double hvezdicka;
01048
01049 if (rh < 0) rh = 0;
01050
01051 CorD(6,kd,0,0,rh,mhmc,mhrh,a3);
01052
01053 hvezdicka = mhmc/(1-sqrt(1-mhrh));
01054 moistakt1 = (1-sqrt(1-rh)) * hvezdicka ;
01055 rh2 = rh +.001;
01056 moistakt2 = (1-sqrt(1-rh2)) * hvezdicka ;
01057
01058 CorD(2,kd,0,0,moistakt1,rs,a2,a3);
01059
01060 ul1 = (moistakt1*rl/rs);
01061 ul2 = (moistakt2*rl/rs);
01062
01063
01064
01065 psi1 = x1;
01066 psi2 = R*x2/(g*M)*log(rh2);
01067
01068 return ((ul2-ul1)/(psi2-psi1));
01069
01070 }
01071
01072 void millymat::CorD(int cislochar, int &kvyhl,double in,int rhw, double x, double & y, double & z, double &z2)
01073 {
01074
01075 switch (MatChar[cislochar])
01076 {
01077 case 0:
01078 break;
01079 case 1: y = MatConst[cislochar];
01080 break;
01081 case 2: {
01082 int i = 1;
01083
01084 int psl, prad;
01085 prad = int(MatData[cislochar][0][0]);
01086 psl = int(MatData[cislochar][1][0]);
01087 if(rhw == 1)
01088 {
01089 if (x <=0)
01090 {
01091 x = 0;
01092 y = MatData[cislochar][0][1];
01093 if (psl == 3)
01094 {
01095 z = MatData[cislochar][2][1];
01096 }
01097
01098 }
01099 else {
01100 if (x > in)
01101 {
01102 x = in;
01103 y = MatData[cislochar][0][prad];
01104 if (psl == 3)
01105 {
01106 z = MatData[cislochar][2][prad];
01107 }
01108 }
01109 else
01110 {
01111 while (x > MatData[cislochar][1][i]){
01112 i++;
01113 }
01114
01115 y = MatData[cislochar][0][i-1] + ((x - MatData[cislochar][1][i-1])*(MatData[cislochar][0][i] - MatData[cislochar][0][i-1]))/(MatData[cislochar][1][i] - MatData[cislochar][1][i-1]);
01116 if (psl == 3)
01117 {
01118 z = MatData[cislochar][2][i-1] + ((x - MatData[cislochar][1][i-1])*(MatData[cislochar][2][i] - MatData[cislochar][2][i-1]))/(MatData[cislochar][1][i] - MatData[cislochar][1][i-1]);
01119 }
01120 }
01121 }
01122 }
01123 else {
01124 if (x <=0)
01125 {
01126 x = 0;
01127 y = MatData[cislochar][1][1];
01128 if (psl == 3)
01129 {
01130 z = MatData[cislochar][2][1];
01131 }
01132
01133 }
01134 else {
01135 if (x > in)
01136 {
01137 x = in;
01138 y = MatData[cislochar][1][prad];
01139 if (psl == 3)
01140 {
01141 z = MatData[cislochar][2][prad];
01142 }
01143 }
01144 else
01145 {
01146 while (x > MatData[cislochar][0][i]){
01147 i++;
01148 }
01149
01150 y = MatData[cislochar][1][i-1] + ((x - MatData[cislochar][0][i-1])*(MatData[cislochar][1][i] - MatData[cislochar][1][i-1]))/(MatData[cislochar][0][i] - MatData[cislochar][0][i-1]);
01151 if (psl == 3)
01152 {
01153 z = MatData[cislochar][2][i-1] + ((x - MatData[cislochar][0][i-1])*(MatData[cislochar][2][i] - MatData[cislochar][2][i-1]))/(MatData[cislochar][0][i] - MatData[cislochar][0][i-1]);
01154 }
01155 }
01156 }
01157 }
01158 kvyhl = 2;
01159 }
01160 break;
01161 case 30: y = MatFunce[cislochar][0];
01162 z = MatFunce[cislochar][1];
01163 kvyhl = 30;
01164 break;
01165 case 31: y = MatFunce[cislochar][0];
01166 z = MatFunce[cislochar][1];
01167 if (cislochar == 6)
01168 {
01169 z2 = MatFunce[cislochar][2];
01170 }
01171 kvyhl = 31;
01172 break;
01173 case 32: y = MatFunce[cislochar][0];
01174 z = MatFunce[cislochar][1];
01175 kvyhl = 32;
01176 break;
01177 case 33:
01178 break;
01179 }
01180 }
01181
01182
01183 void millymat::initvalues (long ipp,long ido)
01184 {
01185 double x1,x2;
01186
01187 x1 = Tm->ip[ipp].av[0];
01188 x2 = Tm->ip[ipp].av[1];
01189
01190 }
01191
01192 void millymat::save_values (long ipp,double *out)
01193 {
01194 Tm->ip[ipp].eqother[0]=out[0];
01195 Tm->ip[ipp].eqother[1]=out[1];
01196 Tm->ip[ipp].eqother[2]=out[2];
01197 Tm->ip[ipp].eqother[3]=out[3];
01198 Tm->ip[ipp].eqother[4]=out[4];
01199 Tm->ip[ipp].eqother[5]=out[5];
01200 Tm->ip[ipp].eqother[6]=out[6];
01201 Tm->ip[ipp].eqother[7]=out[7];
01202 Tm->ip[ipp].eqother[8]=out[8];
01203
01204 }
01205
01206 void millymat::give_values (long ipp,double *av, double *pv, double *eq)
01207 {
01208 av[0] = Tm->ip[ipp].av[0];
01209 av[1] = Tm->ip[ipp].av[1];
01210
01211
01212 eq[0] = Tm->ip[ipp].eqother[0];
01213 eq[1] = Tm->ip[ipp].eqother[1];
01214 eq[2] = Tm->ip[ipp].eqother[2];
01215 eq[3] = Tm->ip[ipp].eqother[3];
01216 eq[4] = Tm->ip[ipp].eqother[4];
01217 eq[5] = Tm->ip[ipp].eqother[5];
01218 eq[6] = Tm->ip[ipp].eqother[6];
01219 eq[7] = Tm->ip[ipp].eqother[7];
01220 eq[8] = Tm->ip[ipp].eqother[8];
01221
01222 }
01223
01224 void millymat::aux_values (long ipp,double *in,double *inp, double *ineq,double *out)
01225 {
01226
01227 double x1,x2;
01228
01229 double ul, w, rh, dudpsi, dudT, kl, lambda,a;
01230
01231 x1 = in[0];
01232 x2 = in[1];
01233
01234
01235
01236
01237
01238
01239
01240
01241
01242
01243
01244
01245 ul = give_ul(x1,x2, rh, dudpsi);
01246 w = ul_to_w(ul);
01247 dudT = derivation_water_retention_curve_by_temperature (x1, x2);
01248
01249 CorD(8,kd,0,0,w,kl,a2,a3);
01250 CorD(10,kd,0,0,w,lambda,a2,a3);
01251 a = relative_volume_ration_a(x1,x2,ul);
01252 if(rh<0.97)
01253 {
01254 ul = 0.0;
01255 }
01256
01257 out[0] = ul;
01258 out[1] = w;
01259 out[2] = rh;
01260 out[8] = dudpsi;
01261 out[4] = dudT;
01262 out[5] = kl;
01263 out[6] = lambda;
01264 out[7] = a;
01265 out[3] = MatConst[7];
01266
01267
01268
01269
01270
01271
01272
01273
01274
01275
01276
01277
01278
01279
01280
01281
01282
01283
01284
01285
01286
01287
01288
01289
01290
01291 }