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