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