00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022 #include <stdio.h>
00023 #include <stdlib.h>
00024 #include <math.h>
00025 #include "aliast.h"
00026 #include "multiphase.h"
00027 #include "constrel.h"
00028 #include "globalt.h"
00029
00030 multiph::multiph()
00031 {
00032 mw = 18.01528;
00033 ma = 28.9645;
00034 gasr = 8314.41;
00035
00036 scale_pc = Tp->scale[0];
00037 scale_pg = Tp->scale[1];
00038 scale_t = Tp->scale[2];
00039
00040
00041
00042 }
00043 multiph::~multiph()
00044 {}
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057 void multiph::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 fprintf (stderr,"\n unknown number of components of stress tensor is required");
00077 fprintf (stderr,"\n in function (%s, line %d).\n",__FILE__,__LINE__);
00078 }
00079 }
00080 }
00081
00082
00083
00084
00085
00086
00087
00088
00089
00090
00091 void multiph::matcond1d (matrix &d,long ri,long ci,long ipp)
00092 {
00093 double k;
00094 double pc,pg,t;
00095 k = 0.0;
00096
00097 pc = Tm->ip[ipp].av[0];
00098 pg = Tm->ip[ipp].av[1];
00099 t = Tm->ip[ipp].av[2];
00100
00101 if((ri == 0) && (ci == 0))
00102 k = get_kcc(pc,pg,t,ipp);
00103 if((ri == 0) && (ci == 1))
00104 k = get_kcg(pc,pg,t,ipp);
00105 if((ri == 0) && (ci == 2))
00106 k = get_kct(pc,pg,t,ipp);
00107
00108 if((ri == 1) && (ci == 0))
00109 k = get_kgc(pc,pg,t,ipp);
00110 if((ri == 1) && (ci == 1))
00111 k = get_kgg(pc,pg,t,ipp);
00112 if((ri == 1) && (ci == 2))
00113 k = get_kgt(pc,pg,t,ipp);
00114
00115 if((ri == 2) && (ci == 0))
00116 k = get_ktc(pc,pg,t,ipp);
00117 if((ri == 2) && (ci == 1))
00118 k = get_ktg(pc,pg,t,ipp);
00119 if((ri == 2) && (ci == 2))
00120 k = get_ktt1(pc,pg,t,ipp);
00121
00122 d[0][0] = k;
00123 }
00124
00125
00126
00127
00128
00129
00130
00131
00132
00133 void multiph::matcond2d (matrix &d,long ri,long ci,long ipp)
00134 {
00135 double k;
00136 double pc,pg,t;
00137 k = 0.0;
00138
00139 pc = Tm->ip[ipp].av[0];
00140 pg = Tm->ip[ipp].av[1];
00141 t = Tm->ip[ipp].av[2];
00142
00143 if((ri == 0) && (ci == 0))
00144 k = get_kcc(pc,pg,t,ipp);
00145 if((ri == 0) && (ci == 1))
00146 k = get_kcg(pc,pg,t,ipp);
00147 if((ri == 0) && (ci == 2))
00148 k = get_kct(pc,pg,t,ipp);
00149
00150 if((ri == 1) && (ci == 0))
00151 k = get_kgc(pc,pg,t,ipp);
00152 if((ri == 1) && (ci == 1))
00153 k = get_kgg(pc,pg,t,ipp);
00154 if((ri == 1) && (ci == 2))
00155 k = get_kgt(pc,pg,t,ipp);
00156
00157 if((ri == 2) && (ci == 0))
00158 k = get_ktc(pc,pg,t,ipp);
00159 if((ri == 2) && (ci == 1))
00160 k = get_ktg(pc,pg,t,ipp);
00161 if((ri == 2) && (ci == 2))
00162 k = get_ktt1(pc,pg,t,ipp);
00163
00164 fillm(0.0,d);
00165
00166 d[0][0] = k; d[0][1] = 0.0;
00167 d[1][0] = 0.0; d[1][1] = k;
00168
00169 fprintf(Outt,"\n i=%ld j=%ld kij=%e",ri,ci,k);
00170
00171 }
00172
00173
00174
00175
00176
00177
00178
00179
00180
00181
00182 void multiph::matcond3d (matrix &d,long ri,long ci,long ipp)
00183 {
00184 double k;
00185 double pc,pg,t;
00186 k = 0.0;
00187
00188 pc = Tm->ip[ipp].av[0];
00189 pg = Tm->ip[ipp].av[1];
00190 t = Tm->ip[ipp].av[2];
00191
00192 if((ri == 0) && (ci == 0))
00193 k = get_kcc(pc,pg,t,ipp);
00194 if((ri == 0) && (ci == 1))
00195 k = get_kcg(pc,pg,t,ipp);
00196 if((ri == 0) && (ci == 2))
00197 k = get_kct(pc,pg,t,ipp);
00198
00199 if((ri == 1) && (ci == 0))
00200 k = get_kgc(pc,pg,t,ipp);
00201 if((ri == 1) && (ci == 1))
00202 k = get_kgg(pc,pg,t,ipp);
00203 if((ri == 1) && (ci == 2))
00204 k = get_kgt(pc,pg,t,ipp);
00205
00206 if((ri == 2) && (ci == 0))
00207 k = get_ktc(pc,pg,t,ipp);
00208 if((ri == 2) && (ci == 1))
00209 k = get_ktg(pc,pg,t,ipp);
00210 if((ri == 2) && (ci == 2))
00211 k = get_ktt1(pc,pg,t,ipp);
00212
00213 fillm(0.0,d);
00214
00215 d[0][0]=k; d[0][1]=0.0; d[0][2]=0.0;
00216 d[1][0]=0.0; d[1][1]=k; d[1][2]=0.0;
00217 d[2][0]=0.0; d[2][1]=0.0; d[2][2]=k;
00218 }
00219
00220
00221
00222
00223
00224
00225
00226
00227
00228
00229 void multiph::matcap (double &c,long ri,long ci,long ipp)
00230 {
00231 double pc,pg,t;
00232 c = 0.0;
00233
00234 pc = Tm->ip[ipp].av[0];
00235 pg = Tm->ip[ipp].av[1];
00236 t = Tm->ip[ipp].av[2];
00237
00238 if((ri == 0) && (ci == 0))
00239 c = get_capcc(pc,pg,t,ipp);
00240 if((ri == 0) && (ci == 1))
00241 c = get_capcg(pc,pg,t,ipp);
00242 if((ri == 0) && (ci == 2))
00243 c = get_capct(pc,pg,t,ipp);
00244
00245 if((ri == 1) && (ci == 0))
00246 c = get_capgc(pc,pg,t,ipp);
00247 if((ri == 1) && (ci == 1))
00248 c = get_capgg(pc,pg,t,ipp);
00249 if((ri == 1) && (ci == 2))
00250 c = get_capgt(pc,pg,t,ipp);
00251
00252 if((ri == 2) && (ci == 0))
00253 c = get_captc(pc,pg,t,ipp);
00254 if((ri == 2) && (ci == 1))
00255 c = get_captg(pc,pg,t,ipp);
00256 if((ri == 2) && (ci == 2))
00257 c = get_captt(pc,pg,t,ipp);
00258 }
00259
00260
00261
00262
00263
00264
00265
00266
00267
00268
00269
00270
00271
00272 void multiph::matcond2 (matrix &d,long ri,long ci,long ipp)
00273 {
00274 long n;
00275 n = d.n;
00276
00277 switch (n){
00278 case 1:{
00279 matcond1d_2 (d,ri,ci,ipp);
00280 break;
00281 }
00282 case 2:{
00283 matcond2d_2 (d,ri,ci,ipp);
00284 break;
00285 }
00286 case 3:{
00287 matcond3d_2 (d,ri,ci,ipp);
00288 break;
00289 }
00290 default:{
00291 fprintf (stderr,"\n unknown number of components of stress tensor is required");
00292 fprintf (stderr,"\n in function (%s, line %d).\n",__FILE__,__LINE__);
00293 }
00294 }
00295 }
00296
00297
00298
00299
00300
00301
00302
00303
00304
00305
00306
00307
00308 void multiph::matcond1d_2 (matrix &d,long ri,long ci,long ipp)
00309 {
00310 double a,b,c,g;
00311 double pc,pg,t;
00312
00313 pc = Tm->ip[ipp].av[0];
00314 pg = Tm->ip[ipp].av[1];
00315 t = Tm->ip[ipp].av[2];
00316
00317 fillm(0.0,d);
00318
00319 if((ri == 2) && (ci == 0)){
00320 }
00321
00322 if((ri == 2) && (ci == 1)){
00323 }
00324
00325 if((ri == 2) && (ci == 2)){
00326 a = get_ktt2a(pc,pg,t,ipp);
00327 b = get_ktt2b(pc,pg,t,ipp);
00328 c = get_ktt2c(pc,pg,t,ipp);
00329 g = Tp->gr1;
00330
00331 d[0][0] = -1.0*(a*(Tm->ip[ipp].grad[1][0] - Tm->ip[ipp].grad[0][0] - b*g) + c*Tm->ip[ipp].grad[1][0]);
00332 }
00333 }
00334
00335
00336
00337
00338
00339
00340
00341
00342
00343
00344 void multiph::matcond2d_2 (matrix &d,long ri,long ci,long ipp)
00345 {
00346 double a,b,c,*g;
00347 double pc,pg,t;
00348
00349 pc = Tm->ip[ipp].av[0];
00350 pg = Tm->ip[ipp].av[1];
00351 t = Tm->ip[ipp].av[2];
00352
00353 fillm(0.0,d);
00354
00355 if((ri == 2) && (ci == 0)){
00356 }
00357 if((ri == 2) && (ci == 1)){
00358 }
00359
00360 if((ri == 2) && (ci == 2)){
00361 a = get_ktt2a(pc,pg,t,ipp);
00362 b = get_ktt2b(pc,pg,t,ipp);
00363 c = get_ktt2c(pc,pg,t,ipp);
00364
00365 g = new double [2];
00366 g[0] = Tp->gr1;
00367 g[1] = Tp->gr2;
00368
00369 d[0][0] = -1.0*(a*(Tm->ip[ipp].grad[1][0] - Tm->ip[ipp].grad[0][0] - b*g[0]) + c*Tm->ip[ipp].grad[1][0]);
00370 d[0][1] = -1.0*(a*(Tm->ip[ipp].grad[1][1] - Tm->ip[ipp].grad[0][1] - b*g[1]) + c*Tm->ip[ipp].grad[1][1]);
00371
00372 delete [] g;
00373 }
00374 }
00375
00376
00377
00378
00379
00380
00381
00382
00383
00384
00385
00386 void multiph::matcond3d_2 (matrix &d,long ri,long ci,long ipp)
00387 {
00388 double a,b,c,*g;
00389 double pc,pg,t;
00390
00391 pc = Tm->ip[ipp].av[0];
00392 pg = Tm->ip[ipp].av[1];
00393 t = Tm->ip[ipp].av[2];
00394
00395 fillm(0.0,d);
00396
00397 if((ri == 2) && (ci == 0)){
00398 }
00399 if((ri == 2) && (ci == 1)){
00400 }
00401
00402 if((ri == 2) && (ci == 2)){
00403 a = get_ktt2a(pc,pg,t,ipp);
00404 b = get_ktt2b(pc,pg,t,ipp);
00405 c = get_ktt2c(pc,pg,t,ipp);
00406
00407 g = new double [3];
00408 g[0] = Tp->gr1;
00409 g[1] = Tp->gr2;
00410 g[2] = Tp->gr3;
00411
00412 d[0][0] = -1.0*(a*(Tm->ip[ipp].grad[1][0] - Tm->ip[ipp].grad[0][0] - b*g[0]) + c*Tm->ip[ipp].grad[1][0]);
00413 d[0][1] = -1.0*(a*(Tm->ip[ipp].grad[1][1] - Tm->ip[ipp].grad[0][1] - b*g[1]) + c*Tm->ip[ipp].grad[1][1]);
00414 d[0][2] = -1.0*(a*(Tm->ip[ipp].grad[1][2] - Tm->ip[ipp].grad[0][2] - b*g[2]) + c*Tm->ip[ipp].grad[1][2]);
00415
00416 delete [] g;
00417 }
00418 }
00419
00420
00421
00422
00423
00424
00425
00426
00427
00428
00429
00430
00431
00432 void multiph::rhs_volume (matrix &d,long ri,long ci,long ipp)
00433 {
00434 long m;
00435 m = d.m;
00436
00437 switch (m){
00438 case 1:{
00439 rhs1d1 (d,ri,ci,ipp);
00440 break;
00441 }
00442 case 2:{
00443 rhs2d1 (d,ri,ci,ipp);
00444 break;
00445 }
00446 case 3:{
00447 rhs3d1 (d,ri,ci,ipp);
00448 break;
00449 }
00450 default:{
00451 fprintf (stderr,"\n unknown number of components of stress tensor is required");
00452 fprintf (stderr,"\n in function (%s, line %d).\n",__FILE__,__LINE__);
00453 }
00454 }
00455 }
00456
00457
00458
00459
00460
00461
00462
00463 void multiph::rhs1d1 (matrix &d,long ri,long ci,long ipp)
00464 {
00465 double f,pc,pg,t;
00466 double g;
00467
00468 pc = Tm->ip[ipp].av[0];
00469 pg = Tm->ip[ipp].av[1];
00470 t = Tm->ip[ipp].av[2];
00471
00472 if(ri == 0){
00473 f = get_fc1(pc,pg,t,ipp);
00474 g = Tp->gr1;
00475
00476 fillm(0.0,d);
00477 d[0][0] = f*g;
00478 }
00479 if(ri == 1){
00480 f = get_fg(pc,pg,t,ipp);
00481 g = Tp->gr1;
00482
00483 fillm(0.0,d);
00484 d[0][0] = f*g;
00485 }
00486 if(ri == 2){
00487 f = get_ft1(pc,pg,t,ipp);
00488 g = Tp->gr1;
00489
00490 fillm(0.0,d);
00491 d[0][0] = f*g;
00492 }
00493 }
00494
00495
00496
00497
00498
00499
00500
00501 void multiph::rhs2d1 (matrix &d,long ri,long ci,long ipp)
00502 {
00503 double f,pc,pg,t;
00504 double *g;
00505
00506 pc = Tm->ip[ipp].av[0];
00507 pg = Tm->ip[ipp].av[1];
00508 t = Tm->ip[ipp].av[2];
00509
00510 if(ri == 0){
00511 g = new double [2];
00512 f = get_fc1(pc,pg,t,ipp);
00513
00514 g[0] = Tp->gr1;
00515 g[1] = Tp->gr2;
00516
00517 fillm(0.0,d);
00518 d[0][0] = f*g[0];
00519 d[1][0] = f*g[1];
00520 }
00521 if(ri == 1){
00522 g = new double [2];
00523 f = get_fg(pc,pg,t,ipp);
00524
00525 g[0] = Tp->gr1;
00526 g[1] = Tp->gr2;
00527
00528 fillm(0.0,d);
00529 d[0][0] = f*g[0];
00530 d[1][0] = f*g[1];
00531 }
00532 if(ri == 2){
00533 g = new double [2];
00534 f = get_ft1(pc,pg,t,ipp);
00535
00536 g[0] = Tp->gr1;
00537 g[1] = Tp->gr2;
00538
00539 fillm(0.0,d);
00540 d[0][0] = f*g[0];
00541 d[1][0] = f*g[1];
00542 }
00543 }
00544
00545
00546
00547
00548
00549
00550
00551 void multiph::rhs3d1 (matrix &d,long ri,long ci,long ipp)
00552 {
00553 double f,pc,pg,t;
00554 double *g;
00555
00556 pc = Tm->ip[ipp].av[0];
00557 pg = Tm->ip[ipp].av[1];
00558 t = Tm->ip[ipp].av[2];
00559
00560 if(ri == 0){
00561 g = new double [3];
00562 f = get_fc1(pc,pg,t,ipp);
00563
00564 g[0] = Tp->gr1;
00565 g[1] = Tp->gr2;
00566 g[2] = Tp->gr3;
00567
00568 fillm(0.0,d);
00569 d[0][0] = f*g[0];
00570 d[1][0] = f*g[1];
00571 d[2][0] = f*g[2];
00572 }
00573 if(ri == 1){
00574 g = new double [2];
00575 f = get_fg(pc,pg,t,ipp);
00576
00577 g[0] = Tp->gr1;
00578 g[1] = Tp->gr2;
00579 g[2] = Tp->gr3;
00580
00581 fillm(0.0,d);
00582 d[0][0] = f*g[0];
00583 d[1][0] = f*g[1];
00584 d[2][0] = f*g[2];
00585 }
00586 if(ri == 2){
00587 g = new double [3];
00588 f = get_ft1(pc,pg,t,ipp);
00589
00590 g[0] = Tp->gr1;
00591 g[1] = Tp->gr2;
00592 g[2] = Tp->gr3;
00593
00594 fillm(0.0,d);
00595 d[0][0] = f*g[0];
00596 d[1][0] = f*g[1];
00597 d[2][0] = f*g[2];
00598 }
00599 }
00600
00601
00602
00603
00604
00605
00606
00607
00608
00609
00610
00611
00612
00613
00614
00615 double multiph::transmission_transcoeff(double trc,long ri,long ci,long nn,long bc,long ipp)
00616 {
00617 long k;
00618 double c,pc,pg,t;
00619 c = 0.0;
00620
00621 k=Gtt->give_dof(nn,0);
00622 if (k>0) {pc = Lsrst->lhs[k-1]+Lsrst->lhsi[k-1];}
00623 if (k==0) {pc = 0.0;}
00624 if (k<0) {pc = Tb->lc[0].pv[0-k-1].getval();}
00625 k=Gtt->give_dof(nn,1);
00626 if (k>0) {pg = Lsrst->lhs[k-1]+Lsrst->lhsi[k-1];}
00627 if (k==0) {pg = 0.0;}
00628 if (k<0) {pg = Tb->lc[0].pv[0-k-1].getval();}
00629 k=Gtt->give_dof(nn,2);
00630 if (k>0) {t = Lsrst->lhs[k-1]+Lsrst->lhsi[k-1];}
00631 if (k==0) {t = 0.0;}
00632 if (k<0) {t = Tb->lc[0].pv[0-k-1].getval();}
00633
00634 if((ri == 0) && (ci == 0))
00635 c = get_transmission_transcoeff_cc(pc,pg,t,bc,ipp);
00636 if((ri == 0) && (ci == 1))
00637 c = 0.0;
00638 if((ri == 0) && (ci == 2))
00639 c = 0.0;
00640
00641 if((ri == 1) && (ci == 0))
00642 c = 0.0;
00643 if((ri == 1) && (ci == 1))
00644 c = 0.0;
00645 if((ri == 1) && (ci == 2))
00646 c = 0.0;
00647
00648 if((ri == 2) && (ci == 0))
00649 c = 0.0;
00650 if((ri == 2) && (ci == 1))
00651 c = 0.0;
00652 if((ri == 2) && (ci == 2))
00653 c = get_transmission_transcoeff_tt(pc,pg,t,bc,ipp);
00654
00655 c = c*trc;
00656
00657 return (c);
00658 }
00659
00660
00661
00662
00663
00664
00665
00666
00667
00668
00669
00670
00671
00672 double multiph::transmission_transcoeff(double trc,long ri,long ci,long nn,long bc,long ipp, int flag)
00673 {
00674 long k;
00675 double c,pc,pg,t;
00676 c = 0.0;
00677
00678 k=Gtt->give_dof(nn,0);
00679 if (k>0) {pc = Lsrst->lhs[k-1]+Lsrst->lhsi[k-1];}
00680 if (k==0) {pc = 0.0;}
00681 if (k<0) {pc = Tb->lc[0].pv[0-k-1].getval();}
00682 k=Gtt->give_dof(nn,1);
00683 if (k>0) {pg = Lsrst->lhs[k-1]+Lsrst->lhsi[k-1];}
00684 if (k==0) {pg = 0.0;}
00685 if (k<0) {pg = Tb->lc[0].pv[0-k-1].getval();}
00686 k=Gtt->give_dof(nn,2);
00687 if (k>0) {t = Lsrst->lhs[k-1]+Lsrst->lhsi[k-1];}
00688 if (k==0) {t = 0.0;}
00689 if (k<0) {t = Tb->lc[0].pv[0-k-1].getval();}
00690
00691 if((ri == 0) && (ci == 0))
00692 c = get_transmission_transcoeff_cc(pc,pg,t,bc,ipp,flag);
00693 if((ri == 0) && (ci == 1))
00694 c = 0.0;
00695 if((ri == 0) && (ci == 2))
00696 c = 0.0;
00697
00698 if((ri == 1) && (ci == 0))
00699 c = 0.0;
00700 if((ri == 1) && (ci == 1))
00701 c = 0.0;
00702 if((ri == 1) && (ci == 2))
00703 c = 0.0;
00704
00705 if((ri == 2) && (ci == 0))
00706 c = 0.0;
00707 if((ri == 2) && (ci == 1))
00708 c = 0.0;
00709 if((ri == 2) && (ci == 2))
00710 c = get_transmission_transcoeff_tt(pc,pg,t,bc,ipp);
00711
00712 c = c*trc;
00713
00714 return (c);
00715 }
00716
00717
00718
00719
00720
00721
00722
00723
00724
00725
00726
00727
00728
00729
00730 double multiph::transmission_nodval(double nodval,double trc2,long ri,long ci,long nn,long bc,long ipp)
00731 {
00732 long k;
00733 double c,pc,pg,t;
00734 c = 0.0;
00735
00736 k=Gtt->give_dof(nn,0);
00737 if (k>0) {pc = Lsrst->lhs[k-1]+Lsrst->lhsi[k-1];}
00738 if (k==0) {pc = 0.0;}
00739 if (k<0) {pc = Tb->lc[0].pv[0-k-1].getval();}
00740 k=Gtt->give_dof(nn,1);
00741 if (k>0) {pg = Lsrst->lhs[k-1]+Lsrst->lhsi[k-1];}
00742 if (k==0) {pg = 0.0;}
00743 if (k<0) {pg = Tb->lc[0].pv[0-k-1].getval();}
00744 k=Gtt->give_dof(nn,2);
00745 if (k>0) {t = Lsrst->lhs[k-1]+Lsrst->lhsi[k-1];}
00746 if (k==0) {t = 0.0;}
00747 if (k<0) {t = Tb->lc[0].pv[0-k-1].getval();}
00748
00749 if((ri == 0) && (ci == 0))
00750 c = get_transmission_nodval_cc(nodval,pc,pg,t,bc,ipp);
00751 if((ri == 0) && (ci == 1))
00752 c = 0.0;
00753 if((ri == 0) && (ci == 2))
00754 c = 0.0;
00755
00756 if((ri == 1) && (ci == 0))
00757 c = 0.0;
00758 if((ri == 1) && (ci == 1))
00759 c = 0.0;
00760 if((ri == 1) && (ci == 2))
00761 c = 0.0;
00762
00763 if((ri == 2) && (ci == 0))
00764 c = 0.0;
00765 if((ri == 2) && (ci == 1))
00766 c = 0.0;
00767 if((ri == 2) && (ci == 2))
00768 c = get_transmission_nodval_tt(nodval,trc2,pc,pg,t,bc,ipp);
00769
00770 return (c);
00771 }
00772
00773
00774
00775
00776
00777
00778
00779
00780
00781
00782
00783
00784
00785
00786
00787 double multiph::transmission_flux(double nodval,double trc2,long ri,long ci,long nn,long bc,long ipp)
00788 {
00789 long k;
00790 double c,pc,pg,t;
00791 c = 0.0;
00792
00793 k=Gtt->give_dof(nn,0);
00794 if (k>0) {pc = Lsrst->lhs[k-1]+Lsrst->lhsi[k-1];}
00795 if (k==0) {pc = 0.0;}
00796 if (k<0) {pc = Tb->lc[0].pv[0-k-1].getval();}
00797 k=Gtt->give_dof(nn,1);
00798 if (k>0) {pg = Lsrst->lhs[k-1]+Lsrst->lhsi[k-1];}
00799 if (k==0) {pg = 0.0;}
00800 if (k<0) {pg = Tb->lc[0].pv[0-k-1].getval();}
00801 k=Gtt->give_dof(nn,2);
00802 if (k>0) {t = Lsrst->lhs[k-1]+Lsrst->lhsi[k-1];}
00803 if (k==0) {t = 0.0;}
00804 if (k<0) {t = Tb->lc[0].pv[0-k-1].getval();}
00805
00806 if((ri == 0) && (ci == 0))
00807 c = get_transmission_flux_cc(nodval,pc,pg,t,bc,ipp);
00808 if((ri == 0) && (ci == 1))
00809 c = 0.0;
00810 if((ri == 0) && (ci == 2))
00811 c = 0.0;
00812
00813 if((ri == 1) && (ci == 0))
00814 c = 0.0;
00815 if((ri == 1) && (ci == 1))
00816 c = 0.0;
00817 if((ri == 1) && (ci == 2))
00818 c = 0.0;
00819
00820 if((ri == 2) && (ci == 0))
00821 c = 0.0;
00822 if((ri == 2) && (ci == 1))
00823 c = 0.0;
00824 if((ri == 2) && (ci == 2))
00825 c = get_transmission_flux_tt(nodval,trc2,pc,pg,t,bc,ipp);
00826
00827 return (c);
00828 }
00829
00830
00831
00832
00833
00834
00835
00836
00837
00838
00839
00840
00841 void multiph::gaspress_check(double pc,double &pg,double t,long ipp)
00842 {
00843 double pgw;
00844 state_eq tt;
00845
00846
00847
00848 pgw = tt.get_pgw(pc,t);
00849 if (pgw > pg)
00850 pg = pgw;
00851
00852
00853
00854
00855 if(pc < 100.0)
00856 pg = 101325.0;
00857
00858
00859 if (pg<0.0){
00860 pg = 0.0;
00861 print_err("gas pressure was modified",__FILE__,__LINE__,__func__);
00862 }
00863
00864 }
00865
00866
00867
00868
00869
00870
00871
00872
00873
00874
00875
00876
00877
00878 void multiph::cappress_check(double &pc,double pg,double t,long ipp)
00879 {
00880
00881
00882
00883
00884
00885
00886
00887
00888
00889
00890
00891
00892
00893
00894
00895
00896
00897
00898
00899
00900
00901
00902
00903
00904
00905
00906
00907
00908 if (pc<0.0){
00909 pc=0.0;
00910 print_err("capillary pressure was modified",__FILE__,__LINE__,__func__);
00911 }
00912
00913 }
00914
00915
00916
00917
00918
00919
00920
00921
00922
00923
00924
00925
00926 void multiph::cappress_stop(double &pc,double pg,double t,long ipp)
00927 {
00928
00929
00930 }
00931
00932
00933
00934
00935
00936
00937
00938 void multiph::values_correction (vector &nv)
00939 {
00940
00941 cappress_check(nv[0],nv[1],nv[2],0);
00942
00943
00944 gaspress_check(nv[0],nv[1],nv[2],0);
00945 }
00946
00947
00948
00949
00950
00951
00952
00953
00954
00955
00956 double multiph::get_kcc(double pc,double pg,double t,long ipp)
00957 {
00958 double rhow,kintr,krw,muw,rhog,deff,dpgw_dpc,kcc,mg,ddbw,ds_dpc;
00959 state_eq tt;
00960
00961
00962 gaspress_check(pc,pg,t,ipp);
00963
00964 cappress_check(pc,pg,t,ipp);
00965
00966 rhow = tt.get_rhow(t);
00967 kintr = tt.get_kintr(pc,pg,t,ipp);
00968 krw = tt.get_krw(pc,t,ipp);
00969 muw = tt.get_muw(t);
00970 rhog = tt.get_rhog(pc,pg,t);
00971 deff = tt.get_deff(pc,pg,t,ipp);
00972 dpgw_dpc = tt.get_dpgw_dpc(pc,t);
00973 mg = tt.get_mg(pc,pg,t);
00974 ds_dpc = tt.get_ds_dpc(pc,t,ipp);
00975 ddbw = tt.get_ddbw(pc,pg,t,ipp);
00976
00977 kcc = rhow*kintr*krw/muw + ddbw*rhow;
00978 kcc = kcc - deff*dpgw_dpc/t/mg;
00979
00980 kcc = -1.0*kcc;
00981
00982 return(kcc);
00983 }
00984
00985
00986
00987
00988
00989
00990
00991
00992
00993
00994 double multiph::get_capcc(double pc,double pg,double t,long ipp)
00995 {
00996 double capcc,s,drhogw_dpc,rhogw,drhow_dpc,rhow,phi,ds_dpc;
00997 state_eq tt;
00998
00999
01000 gaspress_check(pc,pg,t,ipp);
01001
01002 cappress_check(pc,pg,t,ipp);
01003
01004 rhow = tt.get_rhow(t);
01005 phi = tt.get_phi(t,ipp);
01006 ds_dpc = tt.get_ds_dpc(pc,t,ipp);
01007 s = tt.get_s(pc,t,ipp);
01008 drhow_dpc = 0.0;
01009 drhogw_dpc = tt.get_drhogw_dpc(pc,t);
01010 rhogw = tt.get_rhogw(pc,t);
01011
01012 capcc = rhow*phi*ds_dpc + phi*s*drhow_dpc;
01013 capcc = capcc + (1.0 - s)*phi*drhogw_dpc - phi*rhogw*ds_dpc;
01014
01015 return(capcc);
01016 }
01017
01018
01019
01020
01021
01022
01023
01024
01025
01026
01027 double multiph::get_kcg(double pc,double pg,double t,long ipp)
01028 {
01029 double rhogw,kintr,krg,mug,deff,pgw,rhow,krw,muw,kcg,mg;
01030 state_eq tt;
01031
01032
01033 gaspress_check(pc,pg,t,ipp);
01034
01035 cappress_check(pc,pg,t,ipp);
01036
01037 rhogw = tt.get_rhogw(pc,t);
01038 kintr = tt.get_kintr(pc,pg,t,ipp);
01039 krg = tt.get_krg(pc,t,ipp);
01040 mug = tt.get_mug(pc,pg,t);
01041 deff = tt.get_deff(pc,pg,t,ipp);
01042 pgw = tt.get_pgw(pc,t);
01043 rhow = tt.get_rhow(t);
01044 krw = tt.get_krw(pc,t,ipp);
01045 muw = tt.get_muw(t);
01046 mg = tt.get_mg(pc,pg,t);
01047
01048 kcg = -1.0*rhow*kintr*krw/muw;
01049 kcg = -1.0*rhogw*kintr*krg/mug + deff*pgw/pg/t/mg;
01050
01051 kcg = -1.0*kcg;
01052
01053 return(kcg);
01054 }
01055
01056
01057
01058
01059
01060
01061
01062
01063
01064
01065 double multiph::get_capcg(double pc,double pg,double t,long ipp)
01066 {
01067 double capcg,s,phi,drhow_dpg;
01068 state_eq tt;
01069
01070
01071 gaspress_check(pc,pg,t,ipp);
01072
01073 cappress_check(pc,pg,t,ipp);
01074
01075 s = tt.get_s(pc,t,ipp);
01076 phi = tt.get_phi(t,ipp);
01077 drhow_dpg = 0.0;
01078
01079 capcg = phi*s*drhow_dpg;
01080
01081 return(capcg);
01082 }
01083
01084
01085
01086
01087
01088
01089
01090
01091
01092
01093 double multiph::get_kct(double pc,double pg,double t,long ipp)
01094 {
01095 double rhog,deff,dpgw_dt,kct,mg;
01096 state_eq tt;
01097
01098
01099 gaspress_check(pc,pg,t,ipp);
01100
01101 cappress_check(pc,pg,t,ipp);
01102
01103 rhog = tt.get_rhog(pc,pg,t);
01104 deff = tt.get_deff(pc,pg,t,ipp);
01105 dpgw_dt = tt.get_dpgw_dt(pc,t);
01106 mg = tt.get_mg(pc,pg,t);
01107
01108 kct = 0.0;
01109 kct = kct - deff*dpgw_dt/t/mg;
01110
01111 kct = -1.0*kct;
01112
01113 return(kct);
01114 }
01115
01116
01117
01118
01119
01120
01121
01122
01123
01124
01125
01126 double multiph::get_capct(double pc,double pg,double t,long ipp)
01127 {
01128 double capct,betasw,drhow_dt,dphi_dt,dehydw_dt,phi,s,drhogw_dt,ds_dt,rhow,rhogw;
01129 state_eq tt;
01130
01131
01132 gaspress_check(pc,pg,t,ipp);
01133
01134 cappress_check(pc,pg,t,ipp);
01135
01136 phi = tt.get_phi(t,ipp);
01137 s = tt.get_s(pc,t,ipp);
01138 drhogw_dt = tt.get_drhogw_dt(pc,t);
01139 ds_dt = tt.get_ds_dt(pc,t,ipp);
01140 rhow = tt.get_rhow(t);
01141 rhogw = tt.get_rhogw(pc,t);
01142 betasw = 0.0;
01143 drhow_dt = tt.get_drhow_dt(pc,t);
01144 dphi_dt = tt.get_dphi_dt(pc,pg,t,ipp);
01145 dehydw_dt = tt.get_dehydw_dt(pc,pg,t,ipp);
01146
01147 capct = betasw*rhow + rhow*phi*ds_dt + phi*s*drhow_dt + dphi_dt*s*rhow - dehydw_dt;
01148 capct = capct + (1-s)*phi*drhogw_dt - phi*rhogw*ds_dt + dphi_dt*(1.0 -s)*rhogw;
01149
01150 return(capct);
01151 }
01152
01153
01154
01155
01156
01157
01158
01159
01160
01161
01162 double multiph::get_kgg(double pc,double pg,double t,long ipp)
01163 {
01164 double rhoga,kintr,krg,mug,rhog,deff,kgg,mg,pgw;
01165 state_eq tt;
01166
01167
01168 gaspress_check(pc,pg,t,ipp);
01169
01170 cappress_check(pc,pg,t,ipp);
01171
01172 rhoga = tt.get_rhoga(pc,pg,t);
01173 kintr = tt.get_kintr(pc,pg,t,ipp);
01174 krg = tt.get_krg(pc,t,ipp);
01175 mug = tt.get_mug(pc,pg,t);
01176 rhog = tt.get_rhog(pc,pg,t);
01177 deff = tt.get_deff(pc,pg,t,ipp);
01178 mg = tt.get_mg(pc,pg,t);
01179 pgw = tt.get_pgw(pc,t);
01180
01181 kgg = -1.0*rhoga*kintr*krg/mug - deff*pgw/pg/t/mg;
01182
01183 kgg = -1.0*kgg;
01184
01185 return(kgg);
01186 }
01187
01188
01189
01190
01191
01192
01193
01194
01195
01196
01197 double multiph::get_capgg(double pc,double pg,double t,long ipp)
01198 {
01199 double capgg,phi,s,drhoga_dpg;
01200 state_eq tt;
01201
01202
01203 gaspress_check(pc,pg,t,ipp);
01204
01205 cappress_check(pc,pg,t,ipp);
01206
01207 phi = tt.get_phi(t,ipp);
01208 s = tt.get_s(pc,t,ipp);
01209 drhoga_dpg = tt.get_drhoga_dpg(pc,pg,t);
01210
01211 capgg = phi*(1-s)*drhoga_dpg;
01212
01213 return(capgg);
01214 }
01215
01216
01217
01218
01219
01220
01221
01222
01223
01224
01225 double multiph::get_kgc(double pc,double pg,double t,long ipp)
01226 {
01227 double rhog,deff,dpgw_dpc,kgc,mg;
01228 state_eq tt;
01229
01230
01231 gaspress_check(pc,pg,t,ipp);
01232
01233 cappress_check(pc,pg,t,ipp);
01234
01235 rhog = tt.get_rhog(pc,pg,t);
01236 deff = tt.get_deff(pc,pg,t,ipp);
01237 dpgw_dpc = tt.get_dpgw_dpc(pc,t);
01238 mg = tt.get_mg(pc,pg,t);
01239
01240 kgc = deff*dpgw_dpc/t/mg;
01241
01242 kgc = -1.0*kgc;
01243
01244 return(kgc);
01245 }
01246
01247
01248
01249
01250
01251
01252
01253
01254
01255
01256
01257 double multiph::get_capgc(double pc,double pg,double t,long ipp)
01258 {
01259 double capgc,phi,s,drhoga_dpc,rhoga,ds_dpc;
01260 state_eq tt;
01261
01262
01263 gaspress_check(pc,pg,t,ipp);
01264
01265 cappress_check(pc,pg,t,ipp);
01266
01267 phi = tt.get_phi(t,ipp);
01268 s = tt.get_s(pc,t,ipp);
01269 drhoga_dpc = tt.get_drhoga_dpc(pc,pg,t);
01270 rhoga = tt.get_rhoga(pc,pg,t);
01271 ds_dpc = tt.get_ds_dpc(pc,t,ipp);
01272
01273 capgc = phi*(1-s)*drhoga_dpc - phi*rhoga*ds_dpc;
01274
01275 return(capgc);
01276 }
01277
01278
01279
01280
01281
01282
01283
01284
01285
01286
01287 double multiph::get_kgt(double pc,double pg,double t,long ipp)
01288 {
01289 double rhog,deff,dpgw_dt,kgt,mg;
01290 state_eq tt;
01291
01292
01293 gaspress_check(pc,pg,t,ipp);
01294
01295 cappress_check(pc,pg,t,ipp);
01296
01297 rhog = tt.get_rhog(pc,pg,t);
01298 deff = tt.get_deff(pc,pg,t,ipp);
01299 dpgw_dt = tt.get_dpgw_dt(pc,t);
01300 mg = tt.get_mg(pc,pg,t);
01301
01302 kgt = deff*dpgw_dt/t/mg;
01303
01304 kgt = -1.0*kgt;
01305
01306 return(kgt);
01307 }
01308
01309
01310
01311
01312
01313
01314
01315
01316
01317
01318
01319 double multiph::get_capgt(double pc,double pg,double t,long ipp)
01320 {
01321 double capgt,phi,s,drhoga_dt,rhoga,ds_dt,dphi_dt;
01322 state_eq tt;
01323
01324
01325 gaspress_check(pc,pg,t,ipp);
01326
01327 cappress_check(pc,pg,t,ipp);
01328
01329 phi = tt.get_phi(t,ipp);
01330 s = tt.get_s(pc,t,ipp);
01331 drhoga_dt = tt.get_drhoga_dt(pc,pg,t);
01332 rhoga = tt.get_rhoga(pc,pg,t);
01333 ds_dt = tt.get_ds_dt(pc,t,ipp);
01334 dphi_dt = tt.get_dphi_dt(pc,pg,t,ipp);
01335
01336 capgt = phi*(1.0-s)*drhoga_dt - phi*rhoga*ds_dt + dphi_dt*(1.0 - s)*rhoga;
01337
01338 return(capgt);
01339 }
01340
01341
01342
01343
01344
01345
01346
01347
01348
01349
01350 double multiph::get_ktt1(double pc,double pg,double t,long ipp)
01351 {
01352 double lambdaeff,ktt1;
01353 state_eq tt;
01354
01355
01356 gaspress_check(pc,pg,t,ipp);
01357
01358 cappress_check(pc,pg,t,ipp);
01359
01360 lambdaeff = tt.get_lambdaeff(pc,pg,t,ipp);
01361
01362 ktt1 = -1.0*lambdaeff;
01363
01364 ktt1 = -1.0*ktt1;
01365
01366 return(ktt1);
01367 }
01368
01369
01370
01371
01372
01373
01374
01375
01376
01377
01378
01379 double multiph::get_ktt2a(double pc,double pg,double t,long ipp)
01380 {
01381 double cpw,rhow,kintr,krw,muw,ktt2a;
01382 state_eq tt;
01383
01384
01385 gaspress_check(pc,pg,t,ipp);
01386
01387 cappress_check(pc,pg,t,ipp);
01388
01389 cpw = tt.get_cpw();
01390 rhow = tt.get_rhow(t);
01391 kintr = tt.get_kintr(pc,pg,t,ipp);
01392 krw = tt.get_krw(pc,t,ipp);
01393 muw = tt.get_muw(t);
01394
01395 ktt2a = cpw*rhow*kintr*krw/muw;
01396
01397 return(ktt2a);
01398 }
01399
01400
01401
01402
01403
01404
01405
01406
01407
01408
01409 double multiph::get_ktt2b(double pc,double pg,double t,long ipp)
01410 {
01411 double rhow,ktt2b;
01412 state_eq tt;
01413
01414
01415 gaspress_check(pc,pg,t,ipp);
01416
01417 cappress_check(pc,pg,t,ipp);
01418
01419 rhow = tt.get_rhow(t);
01420
01421 ktt2b = rhow;
01422
01423 return(ktt2b);
01424 }
01425
01426
01427
01428
01429
01430
01431
01432
01433
01434
01435 double multiph::get_ktt2c(double pc,double pg,double t,long ipp)
01436 {
01437 double rhocpg,kintr,krg,mug,ktt2c;
01438 state_eq tt;
01439
01440
01441 gaspress_check(pc,pg,t,ipp);
01442
01443 cappress_check(pc,pg,t,ipp);
01444
01445 rhocpg = tt.get_rhocpg(pc,pg,t);
01446 kintr = tt.get_kintr(pc,pg,t,ipp);
01447 krg = tt.get_krg(pc,t,ipp);
01448 mug = tt.get_mug(pc,pg,t);
01449
01450 ktt2c = rhocpg*kintr*krg/mug;
01451
01452 return(ktt2c);
01453 }
01454
01455
01456
01457
01458
01459
01460
01461
01462
01463
01464 double multiph::get_ktt2d(double pc,double pg,double t,long ipp)
01465 {
01466 double rhog,ktt2d;
01467 state_eq tt;
01468
01469
01470 gaspress_check(pc,pg,t,ipp);
01471
01472 cappress_check(pc,pg,t,ipp);
01473
01474 rhog = tt.get_rhog(pc,pg,t);
01475
01476 ktt2d = rhog;
01477
01478 return(ktt2d);
01479 }
01480
01481
01482
01483
01484
01485
01486
01487
01488
01489
01490
01491 double multiph::get_captt(double pc,double pg,double t,long ipp)
01492 {
01493 double captt,s,drhow_dt,dphi_dt,dehydw_dt,fste,hydren,rhocp,dhvap,rhow,phi,ds_dt,betasw;
01494 state_eq tt;
01495
01496
01497 gaspress_check(pc,pg,t,ipp);
01498
01499 cappress_check(pc,pg,t,ipp);
01500
01501 rhocp = tt.get_rhocp(pc,pg,t,ipp);
01502 dhvap = tt.get_dhvap(t);
01503 rhow = tt.get_rhow(t);
01504 phi = tt.get_phi(t,ipp);
01505 ds_dt = tt.get_ds_dt(pc,t,ipp);
01506
01507 betasw = 0.0;
01508 s = tt.get_s(pc,t,ipp);
01509 drhow_dt = tt.get_drhow_dt(pc,t);
01510 dphi_dt = tt.get_dphi_dt(pc,pg,t,ipp);
01511 dehydw_dt = tt.get_dehydw_dt(pc,pg,t,ipp);
01512 fste = tt.get_fste(pc,pg,t,ipp);
01513 hydren = tt.get_hydren(pc,pg,t,ipp);
01514
01515 captt = rhocp - dhvap*(betasw*rhow + rhow*phi*ds_dt + phi*s*drhow_dt + dphi_dt*s*rhow - dehydw_dt) + dehydw_dt*hydren/fste;
01516
01517 return(captt);
01518 }
01519
01520
01521
01522
01523
01524
01525
01526
01527
01528
01529 double multiph::get_ktc(double pc,double pg,double t,long ipp)
01530 {
01531 double dhvap,ddbw,rhow,kintr,krw,muw,ktc,ds_dpc;
01532 state_eq tt;
01533
01534
01535 gaspress_check(pc,pg,t,ipp);
01536
01537 cappress_check(pc,pg,t,ipp);
01538
01539 dhvap = tt.get_dhvap(t);
01540 rhow = tt.get_rhow(t);
01541 kintr = tt.get_kintr(pc,pg,t,ipp);
01542 krw = tt.get_krw(pc,t,ipp);
01543 muw = tt.get_muw(t);
01544 ddbw = tt.get_ddbw(pc,pg,t,ipp);
01545 ds_dpc = tt.get_ds_dpc(pc,t,ipp);
01546
01547 ktc = dhvap*(rhow*kintr*krw/muw + ddbw*rhow);
01548
01549 return(ktc);
01550 }
01551
01552
01553
01554
01555
01556
01557
01558
01559
01560
01561 double multiph::get_captc(double pc,double pg,double t,long ipp)
01562 {
01563 double captc,s,drhow_dpc,dhvap,rhow,phi,ds_dpc;
01564 state_eq tt;
01565
01566
01567 gaspress_check(pc,pg,t,ipp);
01568
01569 cappress_check(pc,pg,t,ipp);
01570
01571 dhvap = tt.get_dhvap(t);
01572 rhow = tt.get_rhow(t);
01573 phi = tt.get_phi(t,ipp);
01574 ds_dpc = tt.get_ds_dpc(pc,t,ipp);
01575 s = tt.get_s(pc,t,ipp);
01576 drhow_dpc = 0.0;
01577
01578 captc = -1.0*dhvap*(rhow*phi*ds_dpc + phi*s*drhow_dpc);
01579
01580 return(captc);
01581 }
01582
01583
01584
01585
01586
01587
01588
01589
01590
01591
01592 double multiph::get_ktg(double pc,double pg,double t,long ipp)
01593 {
01594 double dhvap,rhow,kintr,krw,muw,ktg;
01595 state_eq tt;
01596
01597
01598 gaspress_check(pc,pg,t,ipp);
01599
01600 cappress_check(pc,pg,t,ipp);
01601
01602 dhvap = tt.get_dhvap(t);
01603 rhow = tt.get_rhow(t);
01604 kintr = tt.get_kintr(pc,pg,t,ipp);
01605 krw = tt.get_krw(pc,t,ipp);
01606 muw = tt.get_muw(t);
01607
01608 ktg = -1.0*(dhvap*rhow*kintr*krw/muw);
01609
01610 return(ktg);
01611 }
01612
01613
01614
01615
01616
01617
01618
01619
01620
01621
01622 double multiph::get_captg(double pc,double pg,double t,long ipp)
01623 {
01624 double captg,s,phi,drhow_dpg,dhvap;
01625 state_eq tt;
01626
01627
01628 gaspress_check(pc,pg,t,ipp);
01629
01630 cappress_check(pc,pg,t,ipp);
01631
01632 s = tt.get_s(pc,t,ipp);
01633 phi = tt.get_phi(t,ipp);
01634 drhow_dpg = 0.0;
01635 dhvap = tt.get_dhvap(t);
01636
01637 captg = -1.0*dhvap*phi*s*drhow_dpg;
01638
01639 return(captg);
01640 }
01641
01642
01643
01644
01645
01646
01647
01648
01649
01650 double multiph::get_fc1(double pc,double pg,double t,long ipp)
01651 {
01652 double fc1;
01653 double rhow,rhogw,rhog,kintr,krw,muw,krg,mug;
01654 state_eq tt;
01655
01656
01657 gaspress_check(pc,pg,t,ipp);
01658
01659 cappress_check(pc,pg,t,ipp);
01660
01661 rhow = tt.get_rhow(t);
01662 rhogw = tt.get_rhogw(pc,t);
01663 rhog = tt.get_rhog(pc,pg,t);
01664 kintr = tt.get_kintr(pc,pg,t,ipp);
01665 krw = tt.get_krw(pc,t,ipp);
01666 muw = tt.get_muw(t);
01667 krg = tt.get_krg(pc,t,ipp);
01668 mug = tt.get_mug(pc,pg,t);
01669
01670 fc1 = rhogw*kintr*krg*rhog/mug + rhow*kintr*krw*rhow/muw;
01671
01672 return(fc1);
01673 }
01674
01675
01676
01677
01678
01679
01680
01681
01682
01683
01684 double multiph::get_fg(double pc,double pg,double t,long ipp)
01685 {
01686 double fg;
01687 double rhoga,rhog,kintr,krg,mug;
01688 state_eq tt;
01689
01690
01691 gaspress_check(pc,pg,t,ipp);
01692
01693 cappress_check(pc,pg,t,ipp);
01694
01695 rhoga = tt.get_rhoga(pc,pg,t);
01696 rhog = tt.get_rhog(pc,pg,t);
01697 kintr = tt.get_kintr(pc,pg,t,ipp);
01698 krg = tt.get_krg(pc,t,ipp);
01699 mug = tt.get_mug(pc,pg,t);
01700
01701 fg = rhoga*kintr*krg*rhog/mug;
01702
01703 return(fg);
01704 }
01705
01706
01707
01708
01709
01710
01711
01712
01713
01714
01715 double multiph::get_ft1(double pc,double pg,double t,long ipp)
01716 {
01717 double ft1;
01718 double dhvap,rhow,kintr,krw,muw;
01719 state_eq tt;
01720
01721
01722 gaspress_check(pc,pg,t,ipp);
01723
01724 cappress_check(pc,pg,t,ipp);
01725
01726 dhvap = tt.get_dhvap(t);
01727 rhow = tt.get_rhow(t);
01728 kintr = tt.get_kintr(pc,pg,t,ipp);
01729 krw = tt.get_krw(pc,t,ipp);
01730 muw = tt.get_muw(t);
01731
01732 ft1 = -dhvap*rhow*kintr*krw*rhow/muw;
01733
01734 return(ft1);
01735 }
01736
01737
01738
01739
01740
01741
01742
01743
01744
01745
01746
01747
01748
01749
01750 double multiph::get_transmission_transcoeff_cc(double pc,double pg,double t,long bc,long ipp)
01751 {
01752 double fc3,pgws,rhow;
01753 state_eq tt;
01754
01755
01756 gaspress_check(pc,pg,t,ipp);
01757
01758 cappress_check(pc,pg,t,ipp);
01759
01760 switch (bc){
01761 case 30:{
01762
01763
01764 pgws = tt.get_pgws(t);
01765 rhow = tt.get_rhow(t);
01766
01767 fc3=-1.0*pgws*mw*mw/gasr/gasr/t/t/rhow;
01768
01769 break;
01770 }
01771 case 31:{
01772
01773 pgws = tt.get_pgws(t);
01774 rhow = tt.get_rhow(t);
01775
01776 fc3=-1.0*pgws*mw*mw/gasr/gasr/t/t/rhow;
01777
01778 break;
01779 }
01780 case 32:{
01781
01782 pgws = tt.get_pgws(t);
01783 rhow = tt.get_rhow(t);
01784
01785 fc3=-1.0*pgws*mw*mw/gasr/gasr/t/t/rhow;
01786
01787 break;
01788 }
01789 case 33:{
01790 fc3=0.0;
01791 break;
01792 }
01793 case 34:{
01794 fc3=0.0;
01795 break;
01796 }
01797 case 35:{
01798 fc3=0.0;
01799 break;
01800 }
01801 case 36:{
01802 fc3=0.0;
01803 break;
01804 }
01805 case 90:{
01806 fc3 = 1.0;
01807 break;
01808 }
01809 default:{
01810 fprintf (stderr,"\n\n No real boundary condition is prescribed (%s, line %d).\n",__FILE__,__LINE__);
01811 exit(0);
01812 }
01813 }
01814 return(fc3);
01815 }
01816
01817
01818
01819
01820
01821
01822
01823
01824
01825
01826
01827
01828 double multiph::get_transmission_transcoeff_cc(double pc,double pg,double t,long bc,long ipp,int flag)
01829 {
01830 double fc3,pgws,rhow;
01831 state_eq tt;
01832
01833
01834 gaspress_check(pc,pg,t,ipp);
01835
01836 cappress_check(pc,pg,t,ipp);
01837
01838 switch (bc){
01839 case 30:{
01840
01841
01842
01843
01844
01845 fc3 = 1.0;
01846
01847
01848 break;
01849 }
01850 case 31:{
01851
01852 pgws = tt.get_pgws(t);
01853 rhow = tt.get_rhow(t);
01854
01855 fc3=-1.0*pgws*mw*mw/gasr/gasr/t/t/rhow;
01856
01857 break;
01858 }
01859 case 32:{
01860
01861 pgws = tt.get_pgws(t);
01862 rhow = tt.get_rhow(t);
01863
01864 fc3=-1.0*pgws*mw*mw/gasr/gasr/t/t/rhow;
01865
01866 break;
01867 }
01868 case 33:{
01869 if(flag == 1)
01870 fc3=0.0;
01871 else{
01872 pgws = tt.get_pgws(t);
01873 fc3=1.0;
01874
01875 }
01876 break;
01877 }
01878 case 34:{
01879 fc3=0.0;
01880 break;
01881 }
01882 case 35:{
01883 fc3=0.0;
01884 break;
01885 }
01886 case 36:{
01887 fc3=0.0;
01888 break;
01889 }
01890 case 90:{
01891 fc3 = 1.0;
01892 break;
01893 }
01894 default:{
01895 fprintf (stderr,"\n\n No real boundary condition is prescribed (%s, line %d).\n",__FILE__,__LINE__);
01896 exit(0);
01897 }
01898 }
01899 return(fc3);
01900 }
01901
01902
01903
01904
01905
01906
01907
01908
01909
01910
01911
01912 double multiph::get_transmission_nodval_cc(double bv,double pc,double pg,double t,long bc,long ipp)
01913 {
01914 double new_nodval,fc3,pgws,rhow,rhogw_surf,alpha,temp;
01915 state_eq tt;
01916
01917
01918 gaspress_check(pc,pg,t,ipp);
01919
01920 cappress_check(pc,pg,t,ipp);
01921
01922 switch (bc){
01923 case 30:{
01924
01925
01926
01927
01928
01929
01930
01931 fc3 = tt.get_pcrh(bv,t);
01932
01933
01934
01935
01936
01937 break;
01938 }
01939 case 31:{
01940
01941
01942 pgws = tt.get_pgws(t);
01943 rhow = tt.get_rhow(t);
01944 alpha = 0.30;
01945
01946 fc3 = bv - pgws*mw/gasr/t;
01947 fc3 = fc3 - mw*mw*mw*pgws/t/t/t/gasr/gasr/gasr/rhow/rhow*(exp(-alpha*pc*mw/rhow/t/gasr))*pc*pc/2.0;
01948
01949 break;
01950 }
01951 case 32:{
01952
01953
01954 pgws = tt.get_pgws(t);
01955 rhow = tt.get_rhow(t);
01956 alpha = 0.30;
01957
01958 fc3 = bv*mw/gasr/t;
01959 fc3 = fc3 - pgws*mw/gasr/t;
01960 fc3 = fc3 - mw*mw*mw*pgws/t/t/t/gasr/gasr/gasr/rhow/rhow*(exp(-alpha*pc*mw/rhow/t/gasr))*pc*pc/2.0;
01961
01962 break;
01963 }
01964 case 33:{
01965
01966 pgws = tt.get_pgws(t);
01967 rhogw_surf = tt.get_rhogw(pc,t);
01968
01969 fc3 = bv*pgws*mw/gasr/t;
01970
01971 fc3 = fc3 - rhogw_surf;
01972
01973 break;
01974 }
01975 case 34:{
01976
01977 rhogw_surf = tt.get_rhogw(pc,t);
01978
01979 fc3 = bv - rhogw_surf;
01980 break;
01981 }
01982
01983 case 35:{
01984
01985 rhogw_surf = tt.get_rhogw(pc,t);
01986
01987 fc3 = bv*mw/gasr/t;
01988 fc3 = fc3 - rhogw_surf;
01989 break;
01990 }
01991 case 36:{
01992
01993 pgws = tt.get_pgws(t);
01994 rhogw_surf = tt.get_rhogw(pc,t);
01995
01996 fc3 = bv*pgws*mw/gasr/t;
01997
01998 fc3 = fc3 - rhogw_surf;
01999 break;
02000 }
02001 case 90:{
02002 temp = iso_fire(Tp->time,293.0,0.0);
02003 rhogw_surf = tt.get_rhogw(pc,t);
02004
02005 fc3 = rhogw_surf - bv/temp/gasr*mw;
02006 break;
02007 }
02008 default:{
02009 fprintf (stderr,"\n\n No real boundary condition is prescribed (%s, line %d).\n",__FILE__,__LINE__);
02010 exit(0);
02011 }
02012 }
02013 new_nodval = fc3;
02014 return(new_nodval);
02015 }
02016
02017
02018
02019
02020
02021
02022
02023
02024
02025
02026
02027
02028
02029 double multiph::get_transmission_flux_cc(double bv,double pc,double pg,double t,long bc,long ipp)
02030 {
02031 double flux,fc3,pgws,rhogw_surf,temp;
02032 state_eq tt;
02033
02034
02035 gaspress_check(pc,pg,t,ipp);
02036
02037 cappress_check(pc,pg,t,ipp);
02038
02039 switch (bc){
02040 case 30:{
02041
02042
02043
02044
02045
02046
02047
02048
02049 fc3 = tt.get_pcrh(bv,t);
02050 fc3 = fc3 - pc;
02051
02052
02053 break;
02054 }
02055 case 31:{
02056
02057 rhogw_surf = tt.get_rhogw(pc,t);
02058
02059 fc3 = bv - rhogw_surf;
02060 break;
02061 }
02062
02063 case 32:{
02064
02065 rhogw_surf = tt.get_rhogw(pc,t);
02066
02067 fc3 = bv*mw/gasr/t;
02068 fc3 = fc3 - rhogw_surf;
02069 break;
02070 }
02071 case 33:{
02072
02073 pgws = tt.get_pgws(t);
02074 rhogw_surf = tt.get_rhogw(pc,t);
02075
02076 fc3 = bv*pgws*mw/gasr/t;
02077
02078 fc3 = fc3 - rhogw_surf;
02079
02080 break;
02081 }
02082 case 34:{
02083
02084 rhogw_surf = tt.get_rhogw(pc,t);
02085
02086 fc3 = bv - rhogw_surf;
02087 break;
02088 }
02089
02090 case 35:{
02091
02092 rhogw_surf = tt.get_rhogw(pc,t);
02093
02094 fc3 = bv*mw/gasr/t;
02095 fc3 = fc3 - rhogw_surf;
02096 break;
02097 }
02098 case 36:{
02099
02100 pgws = tt.get_pgws(t);
02101 rhogw_surf = tt.get_rhogw(pc,t);
02102
02103 fc3 = bv*pgws*mw/gasr/t;
02104
02105 fc3 = fc3 - rhogw_surf;
02106 break;
02107 }
02108 case 90:{
02109 temp = iso_fire(Tp->time,293.0,0.0);
02110 pgws = tt.get_pgws(t);
02111 rhogw_surf = tt.get_rhogw(pc,t);
02112
02113 fc3 = bv/temp/gasr*mw;
02114 fc3 = rhogw_surf - fc3;
02115 break;
02116 }
02117 default:{
02118 fprintf (stderr,"\n\n No real boundary condition is prescribed (%s, line %d).\n",__FILE__,__LINE__);
02119 exit(0);
02120 }
02121 }
02122 flux = fc3;
02123 return(flux);
02124 }
02125
02126
02127
02128
02129
02130
02131
02132
02133
02134
02135
02136 double multiph::get_transmission_transcoeff_tt(double pc,double pg,double t,long bc,long ipp)
02137 {
02138 double ft3;
02139
02140
02141 gaspress_check(pc,pg,t,ipp);
02142
02143 cappress_check(pc,pg,t,ipp);
02144
02145 switch (bc){
02146 case 30:{
02147 ft3=1.0;
02148 break;
02149 }
02150 case 31:{
02151 ft3=0.0;
02152 break;
02153 }
02154 case 32:{
02155 ft3=1.0;
02156 break;
02157 }
02158 case 90:{
02159 ft3=1.0;
02160 break;
02161 }
02162 default:{
02163 fprintf (stderr,"\n\n No real boundary condition is prescribed (%s, line %d).\n",__FILE__,__LINE__);
02164 exit(0);
02165 }
02166 }
02167
02168 return(ft3);
02169 }
02170
02171
02172
02173
02174
02175
02176
02177
02178
02179
02180
02181
02182
02183
02184
02185 double multiph::get_transmission_nodval_tt(double bv,double trr,double pc,double pg,double t,long bc,long ipp)
02186 {
02187 double ft3,temp;
02188
02189
02190 gaspress_check(pc,pg,t,ipp);
02191
02192 cappress_check(pc,pg,t,ipp);
02193
02194 switch (bc){
02195 case 30:{
02196 ft3 = bv;
02197 break;
02198 }
02199 case 31:{
02200 ft3 = (bv - t);
02201 ft3 = -1.0*ft3;
02202 break;
02203 }
02204 case 32:{
02205 ft3 = heat_rate(0.0,Tp->time);
02206 break;
02207 }
02208 case 90:{
02209 temp = iso_fire(Tp->time,293.0,0.0);
02210 ft3 = (temp - t) + trr*(temp*temp*temp*temp - t*t*t*t);
02211 ft3 = -1.0*ft3;
02212 break;
02213 }
02214 default:{
02215 fprintf (stderr,"\n\n No real boundary condition is prescribed (%s, line %d).\n",__FILE__,__LINE__);
02216 exit(0);
02217 }
02218 }
02219 return(ft3);
02220 }
02221
02222
02223
02224
02225
02226
02227
02228
02229
02230
02231
02232
02233
02234 double multiph::get_transmission_flux_tt(double bv,double trr,double pc,double pg,double t,long bc,long ipp)
02235 {
02236 double ft3,temp;
02237
02238
02239 gaspress_check(pc,pg,t,ipp);
02240
02241 cappress_check(pc,pg,t,ipp);
02242
02243 switch (bc){
02244 case 30:{
02245 ft3 = (bv - t);
02246 break;
02247 }
02248 case 31:{
02249 ft3 = (bv - t);
02250 break;
02251 }
02252 case 32:{
02253 ft3 = (heat_rate(0.0,Tp->time) - t);
02254 break;
02255 }
02256 case 90:{
02257 temp = iso_fire(Tp->time,293.0,0.0);
02258 ft3 = (temp - t) + trr*(temp*temp*temp*temp - t*t*t*t);
02259 break;
02260 }
02261 default:{
02262 fprintf (stderr,"\n\n No real boundary condition is prescribed (%s, line %d).\n",__FILE__,__LINE__);
02263 exit(0);
02264 }
02265 }
02266
02267 return(ft3);
02268 }
02269
02270
02271
02272
02273
02274
02275
02276 double multiph::heat_rate(double stime,double time)
02277 {
02278 double tinf;
02279
02280 tinf = 298.15 + time/60.0;
02281
02282 if (tinf >= 573.15)
02283 tinf = 573.15;
02284
02285 return(tinf);
02286 }
02287
02288
02289
02290
02291
02292
02293
02294 double multiph::iso_fire (double time, double t0, double tfirestart)
02295 {
02296 double tinf;
02297
02298
02299 tinf=((t0-273.15)+345.0*log10(8.0*(time-tfirestart)/60.0+1.0))+273.15;
02300
02301 return(tinf);
02302 }
02303
02304
02305
02306
02307
02308
02309
02310
02311
02312
02313
02314 double multiph::get_othervalue(long compother,long ipp,double *r)
02315 {
02316 double other;
02317 state_eq tt;
02318
02319 switch (compother){
02320 case 0:{
02321 other = r[0];
02322 break;
02323 }
02324 case 1:{
02325 other = r[1];
02326 break;
02327 }
02328 case 2:{
02329 other = r[2];
02330 break;
02331 }
02332 case 3:{
02333 other = tt.get_rh(r[0],r[2]);
02334 break;
02335 }
02336 case 4:{
02337 other = tt.get_s(r[0],r[2],ipp);
02338 break;
02339 }
02340 case 5:{
02341 other = tt.get_pgw(r[0],r[2]);
02342 break;
02343 }
02344 case 6:{
02345 other = tt.get_pw(r[0],r[1],r[2]);
02346 break;
02347 }
02348 case 7:{
02349 other = tt.get_w(r[0],r[1],r[2],ipp);
02350 break;
02351 }
02352 case 8:{
02353 break;
02354 }
02355 default:{
02356 fprintf (stderr,"\n\n unknown type of component is required in function (%s, line %d).\n",__FILE__,__LINE__);
02357 }
02358 }
02359 return (other);
02360
02361 }
02362
02363
02364
02365
02366
02367
02368 void multiph::print_othervalue_name(FILE *out,long compother)
02369 {
02370 switch (compother){
02371 case 0:{
02372 fprintf (out,"Capilary pressure (Pa) ");
02373 break;
02374 }
02375 case 1:{
02376 fprintf (out,"Gas pressure (Pa) ");
02377 break;
02378 }
02379 case 2:{
02380 fprintf (out,"Temperature (K) ");
02381 break;
02382 }
02383 case 3:{
02384 fprintf (out,"Relative humidity () ");
02385 break;
02386 }
02387 case 4:{
02388 fprintf (out,"Degree of saturation () ");
02389 break;
02390 }
02391 case 5:{
02392 fprintf (out,"Pore water vapor pressure (Pa)");
02393 break;
02394 }
02395 case 6:{
02396 fprintf (out,"Pore water pressure (Pa) ");
02397 break;
02398 }
02399 case 7:{
02400 fprintf (out,"Moisture content (kg/kg) ");
02401 break;
02402 }
02403 case 8:{
02404 fprintf (out,"temp ");
02405 break;
02406 }
02407 default:{
02408 fprintf (stderr,"\n\n unknown type of component is required in function (%s, line %d).\n",__FILE__,__LINE__);
02409 }
02410 }
02411 }