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