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