00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021 #include <stdio.h>
00022 #include <stdlib.h>
00023 #include <math.h>
00024 #include "multiphasec.h"
00025 #include "constrel.h"
00026 #include "constrelcl.h"
00027 #include "constrelcu.h"
00028 #include "globalt.h"
00029 #include "globalc.h"
00030
00031 multiphc::multiphc()
00032 {
00033 scale_pc = Tp->scale[0];
00034 scale_pg = Tp->scale[1];
00035 scale_t = Tp->scale[2];
00036 }
00037 multiphc::~multiphc()
00038 {}
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055 void multiphc::matcond_u (matrix &d,long ri,long ci,long ipp)
00056 {
00057 long m;
00058
00059 m = d.m;
00060
00061 switch (m){
00062 case 1:{
00063 if (ci == 2){
00064 matrix s(d.m,d.m),sd(d.m,d.n);
00065
00066 matcond1d_u (sd,ri,ci,ipp);
00067
00068 Cmu->matstiff (s,ipp);
00069 mxm(s,sd,d);
00070
00071 destrm (s); destrm (sd);
00072 }
00073 else
00074 matcond1d_u (d,ri,ci,ipp);
00075 break;
00076 }
00077 case 3:{
00078 if (ci == 2){
00079 matrix s(d.m,d.m),sd(d.m,d.n);
00080
00081 matcond2d_u (sd,ri,ci,ipp);
00082
00083 Cmu->matstiff (s,ipp);
00084 mxm(s,sd,d);
00085
00086 destrm (s); destrm (sd);
00087 }
00088 else
00089 matcond2d_u (d,ri,ci,ipp);
00090 break;
00091 }
00092 case 6:{
00093 if (ci == 2){
00094 matrix s(d.m,d.m),sd(d.m,d.n);
00095
00096 matcond3d_u (sd,ri,ci,ipp);
00097
00098 Cmu->matstiff (s,ipp);
00099 mxm(s,sd,d);
00100
00101 destrm (s); destrm (sd);
00102 }
00103 else
00104 matcond3d_u (d,ri,ci,ipp);
00105 break;
00106 }
00107 default:{
00108 fprintf (stderr,"\n unknown number of components of stress tensor is required");
00109 fprintf (stderr,"\n in function (%s, line %d).\n",__FILE__,__LINE__);
00110 }
00111 }
00112 }
00113
00114
00115
00116
00117
00118
00119
00120
00121
00122
00123
00124
00125 void multiphc::matcap_u (matrix &d,long ri,long ci,long ipp)
00126 {
00127 long m;
00128
00129 m = d.m;
00130
00131 switch (m){
00132 case 1:{
00133 matcap1d_u (d,ri,ci,ipp);
00134 break;
00135 }
00136 case 3:{
00137 matcap2d_u (d,ri,ci,ipp);
00138 break;
00139 }
00140 case 6:{
00141 matcap3d_u (d,ri,ci,ipp);
00142 break;
00143 }
00144 default:{
00145 fprintf (stderr,"\n unknown number of components of stress tensor is required");
00146 fprintf (stderr,"\n in function (%s, line %d).\n",__FILE__,__LINE__);
00147 }
00148 }
00149 }
00150
00151
00152
00153
00154
00155
00156
00157
00158
00159
00160
00161
00162
00163 void multiphc::matcond1d_u (matrix &d,long ri,long ci,long ipp)
00164 {
00165 double k;
00166 double pc,pg,t;
00167 k = 0.0;
00168
00169 pc = Cmu->ip[ipp].av[0];
00170 pg = Cmu->ip[ipp].av[1];
00171 t = Cmu->ip[ipp].av[2];
00172
00173 if((ri == 0) && (ci == 0))
00174 k = get_kuc(pc,pg,t,ipp);
00175 if((ri == 0) && (ci == 1))
00176 k = get_kug(pc,pg,t,ipp);
00177 if((ri == 0) && (ci == 2))
00178 k = get_kut(pc,pg,t,ipp);
00179
00180 d[0][0] = k;
00181 }
00182
00183
00184
00185
00186
00187
00188
00189
00190
00191
00192
00193
00194 void multiphc::matcond2d_u (matrix &d,long ri,long ci,long ipp)
00195 {
00196 double k;
00197 double pc,pg,t;
00198 k = 0.0;
00199
00200 pc = Cmu->ip[ipp].av[0];
00201 pg = Cmu->ip[ipp].av[1];
00202 t = Cmu->ip[ipp].av[2];
00203
00204 if((ri == 0) && (ci == 0))
00205 k = get_kuc(pc,pg,t,ipp);
00206 if((ri == 0) && (ci == 1))
00207 k = get_kug(pc,pg,t,ipp);
00208 if((ri == 0) && (ci == 2))
00209 k = get_kut(pc,pg,t,ipp);
00210
00211 fillm(0.0,d);
00212
00213 d[0][0] = k;
00214 d[1][0] = k;
00215 d[2][0] = 0.0;
00216 }
00217
00218
00219
00220
00221
00222
00223
00224
00225
00226
00227
00228
00229 void multiphc::matcond3d_u (matrix &d,long ri,long ci,long ipp)
00230 {
00231 double k;
00232 double pc,pg,t;
00233 k = 0.0;
00234
00235 pc = Cmu->ip[ipp].av[0];
00236 pg = Cmu->ip[ipp].av[1];
00237 t = Cmu->ip[ipp].av[2];
00238
00239 if((ri == 0) && (ci == 0))
00240 k = get_kuc(pc,pg,t,ipp);
00241 if((ri == 0) && (ci == 1))
00242 k = get_kug(pc,pg,t,ipp);
00243 if((ri == 0) && (ci == 2))
00244 k = get_kut(pc,pg,t,ipp);
00245
00246 fillm(0.0,d);
00247
00248 d[0][0]=k;
00249 d[1][0]=k;
00250 d[2][0]=k;
00251 d[3][0]=0.0;
00252 d[4][0]=0.0;
00253 d[5][0]=0.0;
00254 }
00255
00256
00257
00258
00259
00260
00261
00262
00263
00264
00265
00266
00267
00268 void multiphc::matcap1d_u (matrix &d,long ri,long ci,long ipp)
00269 {
00270 double c;
00271 double pc,pg,t;
00272 c = 0.0;
00273
00274 pc = Cmu->ip[ipp].av[0];
00275 pg = Cmu->ip[ipp].av[1];
00276 t = Cmu->ip[ipp].av[2];
00277
00278 if((ri == 0) && (ci == 0))
00279 c = get_capuc(pc,pg,t,ipp);
00280 if((ri == 0) && (ci == 1))
00281 c = get_capug(pc,pg,t,ipp);
00282 if((ri == 0) && (ci == 2))
00283 c = get_caput(pc,pg,t,ipp);
00284
00285 fillm(0.0,d);
00286
00287 d[0][0] = c;
00288 }
00289
00290
00291
00292
00293
00294
00295
00296
00297
00298
00299
00300
00301
00302 void multiphc::matcap2d_u (matrix &d,long ri,long ci,long ipp)
00303 {
00304 double c;
00305 double pc,pg,t;
00306 c = 0.0;
00307
00308 pc = Cmu->ip[ipp].av[0];
00309 pg = Cmu->ip[ipp].av[1];
00310 t = Cmu->ip[ipp].av[2];
00311
00312 if((ri == 0) && (ci == 0))
00313 c = get_capuc(pc,pg,t,ipp);
00314 if((ri == 0) && (ci == 1))
00315 c = get_capug(pc,pg,t,ipp);
00316 if((ri == 0) && (ci == 2))
00317 c = get_caput(pc,pg,t,ipp);
00318
00319 fillm(0.0,d);
00320
00321 d[0][0] = c;
00322 d[1][0] = c;
00323 d[2][0] = 0.0;
00324 }
00325
00326
00327
00328
00329
00330
00331
00332
00333
00334
00335
00336
00337
00338 void multiphc::matcap3d_u (matrix &d,long ri,long ci,long ipp)
00339 {
00340 double c;
00341 double pc,pg,t;
00342 c = 0.0;
00343
00344 pc = Cmu->ip[ipp].av[0];
00345 pg = Cmu->ip[ipp].av[1];
00346 t = Cmu->ip[ipp].av[2];
00347
00348 if((ri == 0) && (ci == 0))
00349 c = get_capuc(pc,pg,t,ipp);
00350 if((ri == 0) && (ci == 1))
00351 c = get_capug(pc,pg,t,ipp);
00352 if((ri == 0) && (ci == 2))
00353 c = get_caput(pc,pg,t,ipp);
00354
00355 fillm(0.0,d);
00356
00357 d[0][0] = c;
00358 d[1][0] = c;
00359 d[2][0] = c;
00360 d[3][0] = 0.0;
00361 d[4][0] = 0.0;
00362 d[5][0] = 0.0;
00363 }
00364
00365
00366
00367
00368
00369
00370
00371
00372
00373
00374
00375
00376
00377
00378
00379
00380
00381 void multiphc::matcond_l (matrix &d,long ri,long ci,long ipp)
00382 {
00383 long m;
00384
00385 m = d.n;
00386
00387 switch (m){
00388 case 1:{
00389 matcond1d_l (d,ri,ci,ipp);
00390 break;
00391 }
00392 case 3:{
00393 matcond2d_l (d,ri,ci,ipp);
00394 break;
00395 }
00396 case 6:{
00397 matcond3d_l (d,ri,ci,ipp);
00398 break;
00399 }
00400 default:{
00401 fprintf (stderr,"\n unknown number of components of stress tensor is required");
00402 fprintf (stderr,"\n in function (%s, line %d).\n",__FILE__,__LINE__);
00403 }
00404 }
00405 }
00406
00407
00408
00409
00410
00411
00412
00413
00414
00415
00416
00417
00418 void multiphc::matcap_l (matrix &d,long ri,long ci,long ipp)
00419 {
00420 long m;
00421
00422 m = d.n;
00423
00424 switch (m){
00425 case 1:{
00426 matcap1d_l (d,ri,ci,ipp);
00427 break;
00428 }
00429 case 3:{
00430 matcap2d_l (d,ri,ci,ipp);
00431 break;
00432 }
00433 case 6:{
00434 matcap3d_l (d,ri,ci,ipp);
00435 break;
00436 }
00437 default:{
00438 fprintf (stderr,"\n unknown number of components of stress tensor is required");
00439 fprintf (stderr,"\n in function (%s, line %d).\n",__FILE__,__LINE__);
00440 }
00441 }
00442 }
00443
00444
00445
00446
00447
00448
00449
00450
00451
00452
00453
00454
00455
00456 void multiphc::matcond1d_l (matrix &d,long ri,long ci,long ipp)
00457 {
00458 double k;
00459 double pc,pg,t;
00460 k = 0.0;
00461
00462 pc = Cml->ip[ipp].av[0];
00463 pg = Cml->ip[ipp].av[1];
00464 t = Cml->ip[ipp].av[2];
00465
00466 if((ri == 0) && (ci == 0))
00467 k = get_kcu(pc,pg,t,ipp);
00468 if((ri == 1) && (ci == 0))
00469 k = get_kgu(pc,pg,t,ipp);
00470 if((ri == 2) && (ci == 0))
00471 k = get_ktu(pc,pg,t,ipp);
00472
00473 d[0][0] = k;
00474 }
00475
00476
00477
00478
00479
00480
00481
00482
00483
00484
00485
00486
00487 void multiphc::matcond2d_l (matrix &d,long ri,long ci,long ipp)
00488 {
00489 double k;
00490 double pc,pg,t;
00491 k = 0.0;
00492
00493 pc = Cml->ip[ipp].av[0];
00494 pg = Cml->ip[ipp].av[1];
00495 t = Cml->ip[ipp].av[2];
00496
00497 if((ri == 0) && (ci == 0))
00498 k = get_kcu(pc,pg,t,ipp);
00499 if((ri == 1) && (ci == 0))
00500 k = get_kgu(pc,pg,t,ipp);
00501 if((ri == 2) && (ci == 0))
00502 k = get_ktu(pc,pg,t,ipp);
00503
00504 fillm(0.0,d);
00505
00506
00507 d[0][0] = k;
00508 d[0][1] = k;
00509 d[0][2] = 0.0;
00510 }
00511
00512
00513
00514
00515
00516
00517
00518
00519
00520
00521
00522
00523 void multiphc::matcond3d_l (matrix &d,long ri,long ci,long ipp)
00524 {
00525 double k;
00526 double pc,pg,t;
00527 k = 0.0;
00528
00529 pc = Cml->ip[ipp].av[0];
00530 pg = Cml->ip[ipp].av[1];
00531 t = Cml->ip[ipp].av[2];
00532
00533 if((ri == 0) && (ci == 0))
00534 k = get_kcu(pc,pg,t,ipp);
00535 if((ri == 1) && (ci == 0))
00536 k = get_kgu(pc,pg,t,ipp);
00537 if((ri == 2) && (ci == 0))
00538 k = get_ktu(pc,pg,t,ipp);
00539
00540 fillm(0.0,d);
00541
00542 d[0][0]=k;
00543 d[0][1]=k;
00544 d[0][2]=k;
00545 d[0][3]=0.0;
00546 d[0][4]=0.0;
00547 d[0][5]=0.0;
00548 }
00549
00550
00551
00552
00553
00554
00555
00556
00557
00558
00559
00560
00561
00562 void multiphc::matcap1d_l (matrix &d,long ri,long ci,long ipp)
00563 {
00564 double c;
00565 double pc,pg,t;
00566 c = 0.0;
00567
00568 pc = Cml->ip[ipp].av[0];
00569 pg = Cml->ip[ipp].av[1];
00570 t = Cml->ip[ipp].av[2];
00571
00572 if((ri == 0) && (ci == 0))
00573 c = get_capcu(pc,pg,t,ipp);
00574 if((ri == 1) && (ci == 0))
00575 c = get_capgu(pc,pg,t,ipp);
00576 if((ri == 2) && (ci == 0))
00577 c = get_captu(pc,pg,t,ipp);
00578
00579 fillm(0.0,d);
00580
00581 d[0][0] = c;
00582 }
00583
00584
00585
00586
00587
00588
00589
00590
00591
00592
00593
00594
00595
00596 void multiphc::matcap2d_l (matrix &d,long ri,long ci,long ipp)
00597 {
00598 double c;
00599 double pc,pg,t;
00600 c = 0.0;
00601
00602 pc = Cml->ip[ipp].av[0];
00603 pg = Cml->ip[ipp].av[1];
00604 t = Cml->ip[ipp].av[2];
00605
00606 if((ri == 0) && (ci == 0))
00607 c = get_capcu(pc,pg,t,ipp);
00608 if((ri == 1) && (ci == 0))
00609 c = get_capgu(pc,pg,t,ipp);
00610 if((ri == 2) && (ci == 0))
00611 c = get_captu(pc,pg,t,ipp);
00612
00613 fillm(0.0,d);
00614
00615 d[0][0] = c;
00616 d[0][1] = c;
00617 d[0][2] = 0.0;
00618 }
00619
00620
00621
00622
00623
00624
00625
00626
00627
00628
00629
00630
00631
00632 void multiphc::matcap3d_l (matrix &d,long ri,long ci,long ipp)
00633 {
00634 double c;
00635 double pc,pg,t;
00636 c = 0.0;
00637
00638 pc = Cml->ip[ipp].av[0];
00639 pg = Cml->ip[ipp].av[1];
00640 t = Cml->ip[ipp].av[2];
00641
00642 if((ri == 0) && (ci == 0))
00643 c = get_capcu(pc,pg,t,ipp);
00644 if((ri == 1) && (ci == 0))
00645 c = get_capgu(pc,pg,t,ipp);
00646 if((ri == 2) && (ci == 0))
00647 c = get_captu(pc,pg,t,ipp);
00648
00649 fillm(0.0,d);
00650
00651 d[0][0] = c;
00652 d[0][1] = c;
00653 d[0][2] = c;
00654 d[0][3] = 0.0;
00655 d[0][4] = 0.0;
00656 d[0][5] = 0.0;
00657 }
00658
00659
00660
00661
00662
00663
00664
00665
00666
00667
00668
00669
00670
00671
00672 void multiphc::gaspress_check(double pc,double &pg,double t,long ipp)
00673 {
00674 double pgw;
00675 state_eq tt;
00676
00677
00678
00679 pgw = tt.get_pgw(pc,t);
00680 if (pgw > pg)
00681 pg = pgw;
00682
00683
00684
00685 if(pc < 100.0)
00686 pg = 101325.0;
00687
00688
00689 Tm->ip[ipp].av[1] = pg;
00690 }
00691
00692
00693
00694
00695
00696
00697
00698
00699
00700
00701
00702
00703 void multiphc::cappress_check(double &pc,double pg,double t,long ipp)
00704 {
00705
00706
00707
00708
00709
00710
00711
00712
00713
00714
00715
00716
00717
00718
00719
00720
00721
00722
00723
00724
00725
00726
00727
00728
00729
00730
00731 }
00732
00733
00734
00735
00736
00737
00738
00739
00740
00741
00742
00743
00744 void multiphc::cappress_stop(double &pc,double pg,double t,long ipp)
00745 {
00746
00747
00748 }
00749
00750
00751 void multiphc::values_correction (vector &nv)
00752 {
00753
00754 cappress_check(nv[0],nv[1],nv[2],0);
00755
00756
00757 gaspress_check(nv[0],nv[1],nv[2],0);
00758 }
00759
00760
00761
00762
00763
00764
00765
00766
00767
00768
00769 double multiphc::get_kcu(double pc,double pg,double t,long ipp)
00770 {
00771 double rhow,kintr,krw,muw,kcu;
00772 state_eq tt;
00773 state_eqcl cl;
00774
00775
00776
00777
00778
00779
00780 rhow = tt.get_rhow(t);
00781 kintr = cl.get_kintr(pc,pg,t,ipp);
00782 krw = cl.get_krw(pc,t,ipp);
00783 muw = tt.get_muw(t);
00784
00785 kcu = -1.0*(rhow*kintr*krw/muw*rhow);
00786
00787 return(kcu);
00788 }
00789
00790
00791
00792
00793
00794
00795
00796
00797
00798
00799
00800 double multiphc::get_capcu(double pc,double pg,double t,long ipp)
00801 {
00802 double capcu,s,rhogw,rhow,alpha;
00803 state_eq tt;
00804 state_eqcl cl;
00805
00806 s = cl.get_s(pc,pg,t,ipp);
00807 rhogw = tt.get_rhogw(pc,t);
00808 rhow = tt.get_rhow(t);
00809 alpha = cl.get_alpha(pc,pg,t,ipp);
00810
00811 capcu = ((1.0-s)*rhogw + s*rhow)*alpha;
00812
00813 return(capcu);
00814 }
00815
00816
00817
00818
00819
00820
00821
00822
00823
00824
00825
00826 double multiphc::get_kgu(double pc,double pg,double t,long ipp)
00827 {
00828 double kgu;
00829
00830 kgu = 0.0;
00831
00832 return(kgu);
00833
00834 }
00835
00836
00837
00838
00839
00840
00841
00842
00843 double multiphc::get_capgu(double pc,double pg,double t,long ipp)
00844 {
00845 double capgu,s,rhoga,alpha;
00846 state_eq tt;
00847 state_eqcl cl;
00848
00849 s = cl.get_s(pc,pg,t,ipp);
00850 rhoga = tt.get_rhoga(pc,pg,t);
00851 alpha = cl.get_alpha(pc,pg,t,ipp);
00852
00853 capgu = (1.0-s)*rhoga*alpha;
00854
00855 return(capgu);
00856 }
00857
00858
00859
00860
00861
00862
00863
00864
00865
00866
00867 double multiphc::get_ktu(double pc,double pg,double t,long ipp)
00868 {
00869 double dhvap,rhow,kintr,krw,muw,ktu;
00870 state_eq tt;
00871 state_eqcl cl;
00872
00873
00874
00875
00876
00877
00878 dhvap = tt.get_dhvap(t);
00879 rhow = tt.get_rhow(t);
00880 kintr = cl.get_kintr(pc,pg,t,ipp);
00881 krw = cl.get_krw(pc,t,ipp);
00882 muw = tt.get_muw(t);
00883
00884 ktu = dhvap*(rhow*kintr*krw/muw*rhow);
00885
00886 return(ktu);
00887 }
00888
00889
00890
00891
00892
00893
00894
00895
00896 double multiphc::get_captu(double pc,double pg,double t,long ipp)
00897 {
00898 double captu,dhvap,s,rhow,alpha;
00899 state_eq tt;
00900 state_eqcl cl;
00901
00902 dhvap = tt.get_dhvap(t);
00903 s = cl.get_s(pc,pg,t,ipp);
00904 rhow = tt.get_rhow(t);
00905 alpha = cl.get_alpha(pc,pg,t,ipp);
00906
00907 captu = -1.0*dhvap*s*rhow*alpha;
00908
00909 return(captu);
00910 }
00911
00912
00913
00914
00915
00916
00917
00918
00919
00920
00921
00922
00923 double multiphc::get_kuc(double pc,double pg,double t,long ipp)
00924 {
00925 double alpha,s,kuc;
00926 state_eqcu cu;
00927
00928 alpha = cu.get_alpha(pc,pg,t,ipp);
00929 s = cu.get_s(pc,pg,t,ipp);
00930
00931 kuc = alpha*s;
00932
00933 return(kuc);
00934
00935 }
00936
00937
00938
00939
00940
00941
00942
00943
00944 double multiphc::get_capuc(double pc,double pg,double t,long ipp)
00945 {
00946 double capuc;
00947
00948 capuc = 0.0;
00949
00950 return(capuc);
00951 }
00952
00953
00954
00955
00956
00957
00958
00959
00960
00961
00962
00963
00964 double multiphc::get_kug(double pc,double pg,double t,long ipp)
00965 {
00966 double alpha,kug;
00967 state_eqcu cu;
00968
00969 alpha = cu.get_alpha(pc,pg,t,ipp);
00970
00971 kug = -1.0*alpha;
00972
00973 return(kug);
00974
00975 }
00976
00977
00978
00979
00980
00981
00982
00983
00984
00985 double multiphc::get_capug(double pc,double pg,double t,long ipp)
00986 {
00987 double capug;
00988
00989 capug = 0.0;
00990
00991 return(capug);
00992 }
00993
00994
00995
00996
00997
00998
00999
01000
01001
01002
01003
01004
01005
01006 double multiphc::get_kut(double pc,double pg,double t,long ipp)
01007 {
01008 double betas,kut;
01009 state_eqcu cu;
01010
01011 betas = cu.get_betas(pc,pg,t,ipp);
01012
01013 kut = -1.0*betas/3.0;
01014
01015 return(kut);
01016 }
01017
01018
01019
01020
01021
01022
01023
01024
01025
01026 double multiphc::get_caput(double pc,double pg,double t,long ipp)
01027 {
01028 double caput;
01029
01030 caput = 0.0;
01031
01032 return(caput);
01033 }
01034
01035
01036
01037
01038
01039
01040
01041
01042
01043
01044
01045
01046
01047
01048
01049 void multiphc::rhs_u1 (matrix &d,long ri,long ci,long ipp)
01050 {
01051 long m;
01052
01053 m = d.m;
01054
01055 switch (m){
01056 case 1:{
01057 if (ci == 2){
01058 matrix s(d.m,d.m),sd(d.m,d.n);
01059
01060 rhs1d_u1 (sd,ri,ci,ipp);
01061
01062 Cmu->matstiff (s,ipp);
01063 mxm(s,sd,d);
01064
01065 destrm (s); destrm (sd);
01066 }
01067 else
01068 rhs1d_u1 (d,ri,ci,ipp);
01069 break;
01070 }
01071 case 3:{
01072 if (ci == 2){
01073 matrix s(d.m,d.m),sd(d.m,d.n);
01074
01075 rhs2d_u1 (sd,ri,ci,ipp);
01076
01077 Cmu->matstiff (s,ipp);
01078 mxm(s,sd,d);
01079
01080 destrm (s); destrm (sd);
01081 }
01082 else
01083 rhs2d_u1 (d,ri,ci,ipp);
01084 break;
01085 }
01086 case 6:{
01087 if (ci == 2){
01088 matrix s(d.m,d.m),sd(d.m,d.n);
01089
01090 rhs3d_u1 (sd,ri,ci,ipp);
01091
01092 Cmu->matstiff (s,ipp);
01093 mxm(s,sd,d);
01094
01095 destrm (s); destrm (sd);
01096 }
01097 else
01098 rhs3d_u1 (d,ri,ci,ipp);
01099 break;
01100 }
01101 default:{
01102 fprintf (stderr,"\n unknown number of components of stress tensor is required");
01103 fprintf (stderr,"\n in function (%s, line %d).\n",__FILE__,__LINE__);
01104 }
01105 }
01106 }
01107
01108
01109
01110
01111
01112
01113
01114
01115
01116
01117
01118
01119
01120 void multiphc::rhs1d_u1 (matrix &d,long ri,long ci,long ipp)
01121 {
01122 double k;
01123 double pc,pg,t;
01124 k = 0.0;
01125
01126 pc = Cmu->ip[ipp].av[0];
01127 pg = Cmu->ip[ipp].av[1];
01128 t = Cmu->ip[ipp].av[2];
01129
01130 if((ri == 0) && (ci == 0))
01131 k = get_fuc1(pc,pg,t,ipp);
01132 if((ri == 0) && (ci == 1))
01133 k = get_fug1(pc,pg,t,ipp);
01134 if((ri == 0) && (ci == 2))
01135 k = get_fut1(pc,pg,t,ipp);
01136
01137 d[0][0] = k;
01138 }
01139
01140
01141
01142
01143
01144
01145
01146
01147
01148
01149
01150
01151 void multiphc::rhs2d_u1 (matrix &d,long ri,long ci,long ipp)
01152 {
01153 double k;
01154 double pc,pg,t;
01155 k = 0.0;
01156
01157 pc = Cmu->ip[ipp].av[0];
01158 pg = Cmu->ip[ipp].av[1];
01159 t = Cmu->ip[ipp].av[2];
01160
01161 if((ri == 0) && (ci == 0))
01162 k = get_fuc1(pc,pg,t,ipp);
01163 if((ri == 0) && (ci == 1))
01164 k = get_fug1(pc,pg,t,ipp);
01165 if((ri == 0) && (ci == 2))
01166 k = get_fut1(pc,pg,t,ipp);
01167
01168 fillm(0.0,d);
01169
01170 d[0][0] = k;
01171 d[1][0] = k;
01172 d[2][0] = 0.0;
01173 }
01174
01175
01176
01177
01178
01179
01180
01181
01182
01183
01184
01185
01186 void multiphc::rhs3d_u1 (matrix &d,long ri,long ci,long ipp)
01187 {
01188 double k;
01189 double pc,pg,t;
01190 k = 0.0;
01191
01192 pc = Cmu->ip[ipp].av[0];
01193 pg = Cmu->ip[ipp].av[1];
01194 t = Cmu->ip[ipp].av[2];
01195
01196 if((ri == 0) && (ci == 0))
01197 k = get_fuc1(pc,pg,t,ipp);
01198 if((ri == 0) && (ci == 1))
01199 k = get_fug1(pc,pg,t,ipp);
01200 if((ri == 0) && (ci == 2))
01201 k = get_fut1(pc,pg,t,ipp);
01202
01203 fillm(0.0,d);
01204
01205 d[0][0]=k;
01206 d[1][0]=k;
01207 d[2][0]=k;
01208 d[3][0]=0.0;
01209 d[4][0]=0.0;
01210 d[5][0]=0.0;
01211 }
01212
01213
01214
01215
01216
01217
01218
01219
01220
01221
01222
01223
01224
01225
01226
01227 double multiphc::get_fuc1(double pc,double pg,double t,long ipp)
01228 {
01229 double alpha,s,fuc1;
01230 state_eqcu cu;
01231
01232 alpha = cu.get_alpha(pc,pg,t,ipp);
01233 s = cu.get_s(pc,pg,t,ipp);
01234
01235 fuc1 = alpha*s;
01236
01237 return(fuc1);
01238 }
01239
01240
01241
01242
01243
01244
01245
01246
01247
01248
01249
01250
01251 double multiphc::get_fug1(double pc,double pg,double t,long ipp)
01252 {
01253 double alpha,fug1;
01254 state_eqcu cu;
01255
01256 alpha = cu.get_alpha(pc,pg,t,ipp);
01257
01258 fug1 = -1.0*alpha;
01259
01260 return(fug1);
01261 }
01262
01263
01264
01265
01266
01267
01268
01269
01270
01271
01272
01273
01274
01275 double multiphc::get_fut1(double pc,double pg,double t,long ipp)
01276 {
01277 double fut1,betas;
01278 state_eqcu cu;
01279
01280 betas = cu.get_betas(pc,pg,t,ipp);
01281
01282 fut1 = -1.0*betas/3.0;
01283
01284 return(fut1);
01285 }
01286
01287
01288
01289
01290
01291
01292
01293
01294
01295
01296
01297
01298
01299
01300
01301
01302 void multiphc::rhs_u2 (matrix &d,long ri,long ci,long ipp)
01303 {
01304 long m;
01305
01306 m = d.m;
01307
01308 switch (m){
01309 case 1:{
01310 rhs1d_u2 (d,ri,ci,ipp);
01311 break;
01312 }
01313 case 3:{
01314 rhs2d_u2 (d,ri,ci,ipp);
01315 break;
01316 }
01317 case 6:{
01318 rhs3d_u2 (d,ri,ci,ipp);
01319 break;
01320 }
01321 default:{
01322 fprintf (stderr,"\n unknown number of components of stress tensor is required");
01323 fprintf (stderr,"\n in function (%s, line %d).\n",__FILE__,__LINE__);
01324 }
01325 }
01326 }
01327
01328
01329
01330
01331
01332
01333
01334
01335
01336
01337
01338
01339
01340 void multiphc::rhs1d_u2 (matrix &d,long ri,long ci,long ipp)
01341 {
01342 double k;
01343 double pc,pg,t;
01344 double g;
01345 k = 0.0;
01346
01347 pc = Cmu->ip[ipp].av[0];
01348 pg = Cmu->ip[ipp].av[1];
01349 t = Cmu->ip[ipp].av[2];
01350
01351 k = get_fu2(pc,pg,t,ipp);
01352
01353
01354 g = 0.0;
01355
01356 fillm(0.0,d);
01357 d[0][0] = k*g;
01358 }
01359
01360
01361
01362
01363
01364
01365
01366
01367
01368
01369
01370
01371 void multiphc::rhs2d_u2 (matrix &d,long ri,long ci,long ipp)
01372 {
01373 double k;
01374 double pc,pg,t;
01375 double *g;
01376 k = 0.0;
01377
01378 pc = Cmu->ip[ipp].av[0];
01379 pg = Cmu->ip[ipp].av[1];
01380 t = Cmu->ip[ipp].av[2];
01381
01382 k = get_fu2(pc,pg,t,ipp);
01383
01384 g = new double [2];
01385
01386 g[0] = 0.0;
01387 g[1] = 0.0;
01388
01389 fillm(0.0,d);
01390 d[0][0] = k*g[0];
01391 d[1][0] = k*g[1];
01392 }
01393
01394
01395
01396
01397
01398
01399
01400
01401
01402
01403
01404
01405 void multiphc::rhs3d_u2 (matrix &d,long ri,long ci,long ipp)
01406 {
01407 double k;
01408 double pc,pg,t;
01409 double *g;
01410 k = 0.0;
01411
01412 pc = Cmu->ip[ipp].av[0];
01413 pg = Cmu->ip[ipp].av[1];
01414 t = Cmu->ip[ipp].av[2];
01415
01416 k = get_fu2(pc,pg,t,ipp);
01417
01418 g = new double [3];
01419
01420 g[0] = 0.0;
01421 g[1] = 0.0;
01422 g[2] = 0.0;
01423
01424 fillm(0.0,d);
01425 d[0][0] = k*g[0];
01426 d[1][0] = k*g[1];
01427 d[2][0] = k*g[2];
01428 }
01429
01430
01431
01432
01433
01434
01435
01436
01437
01438
01439
01440
01441 double multiphc::get_fu2(double pc,double pg,double t,long ipp)
01442 {
01443 double fu2,phi,rhos,s,rhow,rhog;
01444 state_eq tt;
01445 state_eqcu cu;
01446
01447 phi = tt.get_phi(t,ipp);
01448 rhos = cu.get_rhos(pc,pg,t,ipp);
01449 s = cu.get_s(pc,pg,t,ipp);
01450 rhow = tt.get_rhow(t);
01451 rhog = tt.get_rhog(pc,pg,t);
01452
01453 fu2 = -1.0*((1.0 - phi)*rhos + phi*s*rhow + phi*(1.0 - s)*rhog);
01454
01455 fu2 = -1.0*fu2;
01456
01457 return(fu2);
01458 }
01459
01460