00001
00002
00003
00004
00005
00006
00007
00008
00009 #include <stdio.h>
00010 #include <stdlib.h>
00011 #include <math.h>
00012 #include "constrel.h"
00013 #include "consol_awf2.h"
00014 #include "globalt.h"
00015
00016 con_awf2mat::con_awf2mat()
00017 {
00018 alpha = 1.0;
00019
00020 ks = 2.167e6;
00021
00022 kw = 2.0e9;
00023
00024 phi0 = 0.297;
00025
00026 kintr = 4.5e-13;
00027
00028 rhow = 1000.0;
00029
00030 muw0 = 1.0e-3;
00031
00032 mug = 1.8e-5;
00033
00034 }
00035
00036 con_awf2mat::~con_awf2mat()
00037 {}
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047 void con_awf2mat::read(XFILE *in)
00048 {
00049 xfscanf (in,"%lf %lf %lf %lf %lf", &alpha, &ks, &phi0, &kw, &kintr);
00050 }
00051
00052
00053
00054
00055
00056
00057
00058
00059
00060 void con_awf2mat::print(FILE *out)
00061 {
00062 fprintf (out," %lf %lf %lf %lf %lf \n",alpha,ks,phi0,kw,kintr);
00063 }
00064
00065
00066
00067
00068
00069
00070
00071
00072
00073
00074
00075
00076
00077
00078 void con_awf2mat::matcond (matrix &d,long ri,long ci,long ipp)
00079 {
00080 long n;
00081 n = d.n;
00082
00083 switch (n){
00084 case 1:{
00085 matcond1d (d,ri,ci,ipp);
00086 break;
00087 }
00088 case 2:{
00089 matcond2d (d,ri,ci,ipp);
00090 break;
00091 }
00092 case 3:{
00093 matcond3d (d,ri,ci,ipp);
00094 break;
00095 }
00096 default:{
00097 fprintf (stderr,"\n unknown number of components of stress tensor is required");
00098 fprintf (stderr,"\n in function (%s, line %d).\n",__FILE__,__LINE__);
00099 }
00100 }
00101 }
00102
00103
00104
00105
00106
00107
00108
00109
00110
00111
00112 void con_awf2mat::matcond1d (matrix &d,long ri,long ci,long ipp)
00113 {
00114 double k;
00115 double pw,pg;
00116 k = 0.0;
00117
00118 pw = Tm->ip[ipp].av[0];
00119 pg = Tm->ip[ipp].av[1];
00120
00121 if((ri == 0) && (ci == 0))
00122 k = get_kww(pw,pg);
00123 if((ri == 0) && (ci == 1))
00124 k = get_kwg(pw,pg);
00125 if((ri == 1) && (ci == 0))
00126 k = get_kgw(pw,pg);
00127 if((ri == 1) && (ci == 1))
00128 k = get_kgg(pw,pg);
00129
00130 d[0][0] = k;
00131 }
00132
00133
00134
00135
00136
00137
00138
00139
00140
00141 void con_awf2mat::matcond2d (matrix &d,long ri,long ci,long ipp)
00142 {
00143 double k;
00144 double pw,pg;
00145 k = 0.0;
00146
00147 pw = Tm->ip[ipp].av[0];
00148 pg = Tm->ip[ipp].av[1];
00149
00150 if((ri == 0) && (ci == 0))
00151 k = get_kww(pw,pg);
00152 if((ri == 0) && (ci == 1))
00153 k = get_kwg(pw,pg);
00154 if((ri == 1) && (ci == 0))
00155 k = get_kgw(pw,pg);
00156 if((ri == 1) && (ci == 1))
00157 k = get_kgg(pw,pg);
00158
00159 fillm(0.0,d);
00160
00161 d[0][0] = k; d[0][1] = 0.0;
00162 d[1][0] = 0.0; d[1][1] = k;
00163 }
00164
00165
00166
00167
00168
00169
00170
00171
00172
00173 void con_awf2mat::matcond3d (matrix &d,long ri,long ci,long ipp)
00174 {
00175 double k;
00176 double pw,pg;
00177 k = 0.0;
00178
00179 pw = Tm->ip[ipp].av[0];
00180 pg = Tm->ip[ipp].av[1];
00181
00182 if((ri == 0) && (ci == 0))
00183 k = get_kww(pw,pg);
00184 if((ri == 0) && (ci == 1))
00185 k = get_kwg(pw,pg);
00186 if((ri == 1) && (ci == 0))
00187 k = get_kgw(pw,pg);
00188 if((ri == 1) && (ci == 1))
00189 k = get_kgg(pw,pg);
00190
00191 fillm(0.0,d);
00192
00193 d[0][0]=k; d[0][1]=0.0; d[0][2]=0.0;
00194 d[1][0]=0.0; d[1][1]=k; d[1][2]=0.0;
00195 d[2][0]=0.0; d[2][1]=0.0; d[2][2]=k;
00196 }
00197
00198
00199
00200
00201
00202
00203
00204
00205
00206
00207 void con_awf2mat::matcap (double &c,long ri,long ci,long ipp)
00208 {
00209 double pw,pg;
00210 c = 0.0;
00211
00212 pw = Tm->ip[ipp].av[0];
00213 pg = Tm->ip[ipp].av[1];
00214
00215 if((ri == 0) && (ci == 0))
00216 c = get_capww(pw,pg);
00217 if((ri == 0) && (ci == 1))
00218 c = get_capwg(pw,pg);
00219 if((ri == 1) && (ci == 0))
00220 c = get_capgw(pw,pg);
00221 if((ri == 1) && (ci == 1))
00222 c = get_capgg(pw,pg);
00223 }
00224
00225
00226
00227
00228
00229
00230
00231
00232
00233
00234
00235 void con_awf2mat::rhs_volume (matrix &d,long ri,long ci,long ipp)
00236 {
00237 long m;
00238 m = d.m;
00239
00240 switch (m){
00241 case 1:{
00242 rhs1d1 (d,ri,ci,ipp);
00243 break;
00244 }
00245 case 2:{
00246 rhs2d1 (d,ri,ci,ipp);
00247 break;
00248 }
00249 case 3:{
00250 rhs3d1 (d,ri,ci,ipp);
00251 break;
00252 }
00253 default:{
00254 fprintf (stderr,"\n unknown number of components of stress tensor is required");
00255 fprintf (stderr,"\n in function (%s, line %d).\n",__FILE__,__LINE__);
00256 }
00257 }
00258 }
00259
00260
00261
00262
00263
00264
00265
00266
00267
00268
00269 void con_awf2mat::rhs1d1 (matrix &d,long ri,long ci,long ipp)
00270 {
00271 double f,pw,pg;
00272 double g;
00273
00274 pw = Tm->ip[ipp].av[0];
00275 pg = Tm->ip[ipp].av[1];
00276
00277 if(ri == 0){
00278 f = get_fw1(pw,pg);
00279 g = Tp->gr1;
00280
00281 fillm(0.0,d);
00282 d[0][0] = f*g;
00283 }
00284 if(ri == 1){
00285 f = get_fg1(pw,pg);
00286 g = Tp->gr1;
00287
00288 fillm(0.0,d);
00289 d[0][0] = f*g;
00290 }
00291 }
00292
00293
00294
00295
00296
00297
00298
00299
00300 void con_awf2mat::rhs2d1 (matrix &d,long ri,long ci,long ipp)
00301 {
00302 double f,pw,pg;
00303 double *g;
00304
00305 pw = Tm->ip[ipp].av[0];
00306 pg = Tm->ip[ipp].av[1];
00307
00308 if(ri == 0){
00309 g = new double [2];
00310 f = get_fw1(pw,pg);
00311
00312 g[0] = Tp->gr1;
00313 g[1] = Tp->gr2;
00314
00315 fillm(0.0,d);
00316 d[0][0] = f*g[0];
00317 d[1][0] = f*g[1];
00318 }
00319 if(ri == 1){
00320 g = new double [2];
00321 f = get_fg1(pw,pg);
00322
00323 g[0] = Tp->gr1;
00324 g[1] = Tp->gr2;
00325
00326 fillm(0.0,d);
00327 d[0][0] = f*g[0];
00328 d[1][0] = f*g[1];
00329 }
00330 }
00331
00332
00333
00334
00335
00336
00337
00338
00339 void con_awf2mat::rhs3d1 (matrix &d,long ri,long ci,long ipp)
00340 {
00341 double f,pw,pg;
00342 double *g;
00343
00344 pw = Tm->ip[ipp].av[0];
00345 pg = Tm->ip[ipp].av[1];
00346
00347 if(ri == 0){
00348 g = new double [3];
00349 f = get_fw1(pw,pg);
00350
00351 g[0] = Tp->gr1;
00352 g[1] = Tp->gr2;
00353 g[2] = Tp->gr3;
00354
00355 fillm(0.0,d);
00356 d[0][0] = f*g[0];
00357 d[1][0] = f*g[1];
00358 d[2][0] = f*g[2];
00359 }
00360 if(ri == 1){
00361 g = new double [2];
00362 f = get_fg1(pw,pg);
00363
00364 g[0] = Tp->gr1;
00365 g[1] = Tp->gr2;
00366 g[2] = Tp->gr3;
00367
00368 fillm(0.0,d);
00369 d[0][0] = f*g[0];
00370 d[1][0] = f*g[1];
00371 d[2][0] = f*g[2];
00372 }
00373 }
00374
00375
00376
00377
00378
00379
00380
00381
00382
00383
00384
00385 double con_awf2mat::get_sw(double pw, double pg)
00386 {
00387 double sw,pc;
00388
00389
00390 pc = -pw;
00391
00392 sw = 1.0 - (1.9722e-11)*pow(pc,2.4279);
00393
00394 return(sw);
00395 }
00396
00397
00398
00399
00400
00401
00402
00403
00404
00405
00406
00407
00408 double con_awf2mat::get_cs(double pw, double pg)
00409 {
00410 double cs,pc;
00411
00412
00413 pc = -pw;
00414
00415 cs = -2.4279*(1.9722e-11)*pow(pc,(2.4279-1.0));
00416
00417 return(cs);
00418 }
00419
00420
00421
00422
00423
00424
00425
00426
00427
00428
00429
00430
00431 double con_awf2mat::get_se(double sw)
00432 {
00433 double se;
00434
00435 se = (sw - 0.2)*(1.0 - 0.2);
00436
00437 return(se);
00438 }
00439
00440
00441
00442
00443
00444
00445
00446
00447
00448
00449
00450
00451
00452 double con_awf2mat::get_krl(double pw, double pg)
00453 {
00454 double sw,krl;
00455
00456 sw = get_sw(pw,pg);
00457 krl = 1.0 + 2207*pow((1.0-sw),1.0121);
00458
00459 return(krl);
00460 }
00461
00462
00463
00464
00465
00466
00467
00468
00469
00470
00471
00472
00473 double con_awf2mat::get_krg(double pw, double pg)
00474 {
00475 double krg,sw,se;
00476
00477 sw = get_sw(pw,pg);
00478 se = get_se(sw);
00479 krg = pow((1.0-se),2.0)*(1.0 - pow(se,(5.0/3.0)));
00480
00481 if((krg < 0.0001) || (sw >= 0.998))
00482 krg = 0.0001;
00483
00484 return(krg);
00485 }
00486
00487
00488
00489
00490
00491
00492
00493
00494
00495
00496
00497
00498 double con_awf2mat::get_kww(double pw, double pg)
00499 {
00500 double krw,kww;
00501
00502 krw = get_krl(pw,pg);
00503 kww = krw*kintr/muw0;
00504
00505 return(kww);
00506 }
00507
00508
00509
00510
00511
00512
00513
00514
00515
00516
00517
00518 double con_awf2mat::get_kwg(double pw, double pg)
00519 {
00520 double kwg;
00521
00522 kwg = 0.0;
00523
00524 return(kwg);
00525
00526 }
00527
00528
00529
00530
00531
00532
00533
00534
00535
00536
00537
00538 double con_awf2mat::get_kgw(double pw, double pg)
00539 {
00540 double kgw;
00541
00542 kgw = 0.0;
00543
00544 return(kgw);
00545
00546 }
00547
00548
00549
00550
00551
00552
00553
00554
00555
00556
00557
00558 double con_awf2mat::get_kgg(double pw, double pg)
00559 {
00560 double krg,kgg;
00561
00562 krg = get_krg(pw,pg);
00563 kgg = krg*kintr/mug;
00564
00565 return(kgg);
00566 }
00567
00568
00569
00570
00571
00572
00573
00574
00575
00576
00577
00578
00579 double con_awf2mat::get_capww(double pw, double pg)
00580 {
00581 double sw,cs,n,capww;
00582
00583 sw = get_sw(pw,pg);
00584 cs = get_cs(pw,pg);
00585 n = phi0;
00586
00587
00588 capww = (alpha-n)/ks*sw*(sw - pw*cs + pg*cs) + n*sw/kw - n*cs;
00589
00590 return(capww);
00591
00592 }
00593
00594
00595
00596
00597
00598
00599
00600
00601
00602
00603
00604 double con_awf2mat::get_capwg(double pw, double pg)
00605 {
00606 double sw,sg,cs,n,capwg;
00607
00608 sw = get_sw(pw,pg);
00609 sg = 1.0 - sw;
00610 cs = get_cs(pw,pg);
00611 n = phi0;
00612
00613
00614 capwg = (alpha-n)/ks*sw*(sg + pw*cs - pg*cs) + n*cs;
00615
00616 return(capwg);
00617
00618 }
00619
00620
00621
00622
00623
00624
00625
00626
00627
00628
00629
00630 double con_awf2mat::get_capgw(double pw, double pg)
00631 {
00632 double sw,sg,cs,n,capgw;
00633
00634 sw = get_sw(pw,pg);
00635 sg = 1.0 - sw;
00636 cs = get_cs(pw,pg);
00637 n = phi0;
00638
00639
00640 capgw = (alpha-n)/ks*sg*(sw + pg*cs - pw*cs) + n*cs;
00641
00642 return(capgw);
00643 }
00644
00645
00646
00647
00648
00649
00650
00651
00652
00653
00654
00655 double con_awf2mat::get_capgg(double pw, double pg)
00656 {
00657 double sw,sg,cs,n,capgg;
00658
00659 sw = get_sw(pw,pg);
00660 sg = 1.0 - sw;
00661 cs = get_cs(pw,pg);
00662 n = phi0;
00663
00664
00665 capgg = (alpha-n)/ks*sg*(sg - pg*cs + pw*cs) + n*sg/pg - n*cs;
00666
00667 return(capgg);
00668
00669 }
00670
00671
00672
00673
00674
00675
00676
00677
00678
00679
00680
00681 double con_awf2mat::get_fw1(double pw, double pg)
00682 {
00683 double fw1,krl;
00684
00685 krl = get_krl(pw,pg);
00686
00687 fw1 = krl*kintr/muw0*rhow;
00688
00689 return(fw1);
00690
00691 }
00692
00693
00694
00695
00696
00697
00698
00699
00700
00701
00702
00703 double con_awf2mat::get_fg1(double pw, double pg)
00704 {
00705 double fg1,krg,rhog;
00706
00707 krg = get_krg(pw,pg);
00708 rhog = 1.25;
00709
00710 fg1 = krg*kintr/mug*rhog;
00711
00712 return(fg1);
00713
00714 }
00715
00716
00717
00718
00719
00720
00721
00722
00723
00724
00725
00726
00727
00728
00729 double con_awf2mat::get_othervalue(long compother,double pw, double pg, long ipp)
00730 {
00731 double other;
00732 state_eq tt;
00733
00734 switch (compother){
00735 case 0:{
00736 other = pg - pw;
00737 break;
00738 }
00739 case 1:{
00740 other = pg;
00741 break;
00742 }
00743 case 2:{
00744 other = get_sw(pw,pg);
00745 break;
00746 }
00747 case 3:{
00748 other = pw;
00749 break;
00750 }
00751 case 4:{
00752 other = get_cs(pw,pg);
00753 break;
00754 }
00755 default:{
00756 fprintf (stderr,"\n\n unknown type of component is required in function (%s, line %d).\n",__FILE__,__LINE__);
00757 }
00758 }
00759 return (other);
00760
00761 }
00762
00763
00764
00765
00766
00767
00768
00769
00770
00771
00772 void con_awf2mat::print_othervalue_name(FILE *out,long compother)
00773 {
00774 switch (compother){
00775 case 0:{
00776 fprintf (out,"Capillary pressure (Pa)");
00777 break;
00778 }
00779 case 1:{
00780 fprintf (out,"Gas pressure (Pa) ");
00781 break;
00782 }
00783 case 2:{
00784 fprintf (out,"Degree of saturation () ");
00785 break;
00786 }
00787 case 3:{
00788 fprintf (out,"Pore water pressure (Pa) ");
00789 break;
00790 }
00791 case 4:{
00792 fprintf (out,"Moisture content (kg/kg) ");
00793 break;
00794 }
00795 default:{
00796 fprintf (stderr,"\n\n unknown type of component is required in function (%s, line %d).\n",__FILE__,__LINE__);
00797 }
00798 }
00799 }
00800
00801
00802
00803
00804
00805
00806
00807
00808
00809
00810
00811 void con_awf2mat::values_correction (vector &nv, long ipp)
00812 {
00813
00814 water_pressure_check(nv[0],nv[1],ipp);
00815 }
00816
00817
00818
00819
00820
00821
00822
00823
00824
00825
00826 void con_awf2mat::water_pressure_check(double &pw,double &pg,long ipp)
00827 {
00828 if (pw >= 0.0)
00829 pw = -1.0;
00830
00831 if (pg <= 0.0)
00832 pg = 1.0;
00833
00834
00835 Tm->ip[ipp].av[0] = pw;
00836 Tm->ip[ipp].av[1] = pg;
00837 }