00001
00002
00003
00004
00005
00006
00007
00008 #include <stdio.h>
00009 #include <stdlib.h>
00010 #include <math.h>
00011 #include "constrel.h"
00012 #include "sejtkrmat.h"
00013 #include "globalt.h"
00014
00015 sejtkrmat::sejtkrmat()
00016 {
00017
00018 emod = 10.0e6;
00019
00020 nu = 0.4;
00021
00022 alpha = 0.5;
00023
00024 kz = 18.0e12;
00025
00026 ks = 3.6e6;
00027
00028 kk = 2.0e9;
00029
00030 phi0 = 0.5;
00031
00032 k = 1.13e-10;
00033
00034 rhok = 998.0;
00035
00036 g = 9.806;
00037
00038 lambda = 0.0;
00039
00040 c = 4.1667e-8;
00041 }
00042
00043 sejtkrmat::~sejtkrmat()
00044 {}
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057 void sejtkrmat::matcond (matrix &d,long ri,long ci,long ipp)
00058 {
00059 long n;
00060 n = d.n;
00061
00062 switch (n){
00063 case 1:{
00064 matcond1d (d,ri,ci,ipp);
00065 break;
00066 }
00067 case 2:{
00068 matcond2d (d,ri,ci,ipp);
00069 break;
00070 }
00071 case 3:{
00072 matcond3d (d,ri,ci,ipp);
00073 break;
00074 }
00075 default:{
00076 print_err("unknown number of components of conductivity tensor is required",__FILE__,__LINE__,__func__);
00077 }
00078 }
00079 }
00080
00081
00082
00083
00084
00085
00086
00087
00088
00089 void sejtkrmat::matcond1d (matrix &d,long ri,long ci,long ipp)
00090 {
00091 double kk;
00092 double pw;
00093
00094 pw = Tm->ip[ipp].av[0];
00095 kk = get_kww(pw);
00096
00097 fillm(0.0,d);
00098
00099 d[0][0] = kk;
00100 }
00101
00102
00103
00104
00105
00106
00107
00108
00109
00110 void sejtkrmat::matcond2d (matrix &d,long ri,long ci,long ipp)
00111 {
00112 double kk;
00113 double pw;
00114
00115 pw = Tm->ip[ipp].av[0];
00116 kk = get_kww(pw);
00117
00118 fillm(0.0,d);
00119
00120 d[0][0] = kk; d[0][1] = 0.0;
00121 d[1][0] = 0.0; d[1][1] = kk;
00122 }
00123
00124
00125
00126
00127
00128
00129
00130
00131
00132
00133 void sejtkrmat::matcond3d (matrix &d,long ri,long ci,long ipp)
00134 {
00135 double kk;
00136 double pw;
00137
00138 pw = Tm->ip[ipp].av[0];
00139 kk = get_kww(pw);
00140
00141 fillm(0.0,d);
00142
00143 d[0][0]=kk; d[0][1]=0.0; d[0][2]=0.0;
00144 d[1][0]=0.0; d[1][1]=kk; d[1][2]=0.0;
00145 d[2][0]=0.0; d[2][1]=0.0; d[2][2]=kk;
00146 }
00147
00148
00149
00150
00151
00152
00153
00154
00155
00156
00157 void sejtkrmat::matcap (double &cc,long ri,long ci,long ipp)
00158 {
00159 double pw;
00160
00161 pw = Tm->ip[ipp].av[0];
00162 cc = get_capww(pw);
00163 }
00164
00165
00166
00167
00168
00169
00170
00171
00172 void sejtkrmat::read(XFILE *in)
00173 {
00174
00175 alpha = 1.0 - ks/kz;
00176 }
00177
00178
00179
00180
00181
00182
00183
00184
00185
00186
00187
00188 void sejtkrmat::rhs_volume (matrix &d,long ri,long ci,long ipp)
00189 {
00190 long m;
00191 m = d.m;
00192
00193 switch (m){
00194 case 1:{
00195 rhs1d1 (d,ri,ci,ipp);
00196 break;
00197 }
00198 case 2:{
00199 rhs2d1 (d,ri,ci,ipp);
00200 break;
00201 }
00202 case 3:{
00203 rhs3d1 (d,ri,ci,ipp);
00204 break;
00205 }
00206 default:{
00207 fprintf (stderr,"\n unknown number of components of stress tensor is required");
00208 fprintf (stderr,"\n in function (%s, line %d).\n",__FILE__,__LINE__);
00209 }
00210 }
00211 }
00212
00213
00214
00215
00216
00217
00218
00219
00220
00221 void sejtkrmat::rhs1d1 (matrix &d,long ri,long ci,long ipp)
00222 {
00223 double f,pw;
00224 double g;
00225
00226 pw = Tm->ip[ipp].av[0];
00227
00228 f = get_fw1(pw);
00229 g = Tp->gr1;
00230
00231 fillm(0.0,d);
00232 d[0][0] = f*g;
00233 }
00234
00235
00236
00237
00238
00239
00240
00241
00242 void sejtkrmat::rhs2d1 (matrix &d,long ri,long ci,long ipp)
00243 {
00244 double f,pw;
00245 double *g;
00246
00247 pw = Tm->ip[ipp].av[0];
00248
00249 g = new double [2];
00250 f = get_fw1(pw);
00251
00252 g[0] = Tp->gr1;
00253 g[1] = Tp->gr2;
00254
00255 fillm(0.0,d);
00256 d[0][0] = f*g[0];
00257 d[1][0] = f*g[1];
00258 }
00259
00260
00261
00262
00263
00264
00265
00266
00267 void sejtkrmat::rhs3d1 (matrix &d,long ri,long ci,long ipp)
00268 {
00269 double f,pw;
00270 double *g;
00271
00272 pw = Tm->ip[ipp].av[0];
00273
00274 g = new double [3];
00275 f = get_fw1(pw);
00276
00277 g[0] = Tp->gr1;
00278 g[1] = Tp->gr2;
00279 g[2] = Tp->gr3;
00280
00281 fillm(0.0,d);
00282 d[0][0] = f*g[0];
00283 d[1][0] = f*g[1];
00284 d[2][0] = f*g[2];
00285 }
00286
00287
00288
00289
00290
00291
00292
00293
00294
00295
00296
00297 double sejtkrmat::get_kuw(double pw)
00298 {
00299 double kuw;
00300
00301 kuw = -alpha;
00302
00303 return(kuw);
00304
00305 }
00306
00307
00308
00309
00310
00311
00312
00313
00314
00315
00316 double sejtkrmat::get_kwu(double pw)
00317 {
00318 double kwu;
00319
00320 kwu = 0.0;
00321
00322 return(kwu);
00323
00324 }
00325
00326
00327
00328
00329
00330
00331
00332
00333
00334
00335 double sejtkrmat::get_kww(double pw)
00336 {
00337 double kww;
00338
00339 kww = c/ks;
00340
00341 return(kww);
00342
00343 }
00344
00345
00346
00347
00348
00349
00350
00351
00352
00353
00354 double sejtkrmat::get_capuw(double pw)
00355 {
00356 double capuw;
00357
00358 capuw = 0.0;
00359
00360 return(capuw);
00361
00362 }
00363
00364
00365
00366
00367
00368
00369
00370
00371
00372
00373 double sejtkrmat::get_capwu(double pw)
00374 {
00375 double capwu;
00376
00377 capwu = alpha;
00378
00379 return(capwu);
00380
00381 }
00382
00383
00384
00385
00386
00387
00388
00389
00390
00391
00392 double sejtkrmat::get_capww(double pw)
00393 {
00394 double capww;
00395
00396 lambda = phi0/kk + (1-phi0)/kz - (1 - alpha)/kz;
00397
00398 capww = lambda;
00399
00400 return(capww);
00401
00402 }
00403
00404
00405
00406
00407
00408
00409
00410
00411
00412
00413 double sejtkrmat::get_fw1(double pw)
00414 {
00415 double fw1;
00416
00417 fw1 = c/ks;
00418
00419 return(fw1);
00420
00421 }
00422
00423
00424
00425
00426
00427
00428
00429
00430
00431
00432
00433 double sejtkrmat::get_fu1(double pw)
00434 {
00435 double fu1;
00436
00437 fu1 = 0.0;
00438
00439 return(fu1);
00440
00441 }
00442
00443
00444
00445
00446
00447
00448
00449
00450
00451
00452
00453
00454
00455 double sejtkrmat::get_othervalue(long compother,double pw, long ipp)
00456 {
00457 double other;
00458 state_eq tt;
00459
00460 switch (compother){
00461 case 0:{
00462 other = -pw;
00463 break;
00464 }
00465 case 1:{
00466 other = 1.0;
00467 break;
00468 }
00469 case 2:{
00470 other = pw;
00471 break;
00472 }
00473 default:{
00474 fprintf (stderr,"\n\n unknown type of component is required in function (%s, line %d).\n",__FILE__,__LINE__);
00475 }
00476 }
00477 return (other);
00478
00479 }
00480
00481
00482
00483
00484
00485
00486
00487
00488
00489
00490 void sejtkrmat::print_othervalue_name(FILE *out,long compother)
00491 {
00492 switch (compother){
00493 case 0:{
00494 fprintf (out,"Capillary pressure (Pa)");
00495 break;
00496 }
00497 case 1:{
00498 fprintf (out,"Degree of saturation () ");
00499 break;
00500 }
00501 case 2:{
00502 fprintf (out,"Pore water pressure (Pa) ");
00503 break;
00504 }
00505 default:{
00506 fprintf (stderr,"\n\n unknown type of component is required in function (%s, line %d).\n",__FILE__,__LINE__);
00507 }
00508 }
00509 }
00510
00511
00512
00513
00514
00515
00516
00517
00518
00519
00520
00521 void sejtkrmat::values_correction (vector &nv, long ipp)
00522 {
00523
00524 water_pressure_check(nv[0],ipp);
00525 }
00526
00527
00528
00529
00530
00531
00532
00533
00534
00535
00536 void sejtkrmat::water_pressure_check(double &pw,long ipp)
00537 {
00538 if (pw >= 0.0)
00539 pw = -1.0;
00540
00541
00542 Tm->ip[ipp].av[0] = pw;
00543 }