📄 boastview.cpp
字号:
rw1=(rhoscw + rsw1*rhoscg)/bw[i-1][j][k];
rg1=rhoscg/bg[i-1][j][k];
fact=-d288*(el[i-1][j][k]-el[i][j][k]);
gow1=(ro1+ro)*fact;
gww1=(rw1+rw)*fact + pcow-pcow1;
ggw1=(rg1+rg)*fact + pcgo1-pcgo;
p11=p1-pp;
ho1=p11+gow1;
hw1=p11+gww1;
hg1=p11+ggw1;
if(ho1>=0.)kro1=interp(sat,krot,msat,so1s,kro1);
if(ho1<0.) kro1=interp(sat,krot,msat,sso,kro1);
if(hw1>=0.)krw1=interp(sat,krwt,msat,sw1s,krw1);
if(hw1<0.) krw1=interp(sat,krwt,msat,ssw,krw1);
if(hg1>=0.)krg1=interp(sat,krgt,msat,sg1s,krg1);
if(hg1<0.) krg1=interp(sat,krgt,msat,ssg,krg1);
mo1=4.0*kro1/((bo[i-1][j][k]+bo[i][j][k]) * (muo1+muo));
mw1=4.0*krw1/((bw[i-1][j][k]+bw[i][j][k]) * (muw1+muw));
mg1=4.0*krg1/((bg[i-1][j][k]+bg[i][j][k]) * (mug1+mug));
M115:
aow=tx[i][j][k]*mo1;
aww=tx[i][j][k]*mw1;
agw=tx[i][j][k]*mg1;
if(i==ii)goto M125;
p2=p[i+1][j][k];
bpt=pbot[i+1][j][k];
rso2=intpvt(bpt,rslope,pot,rsot,mpot,p2,rso2);
muo2=intpvt(bpt,vslope,pot,muot,mpot,p2,muo2);
rsw2=interp(pwt,rswt,mpwt,p2,rsw2);
muw2=interp(pwt,muwt,mpwt,p2,muw2);
mug2=interp(pgt,mugt,mpgt,p2,mug2);
so2=so[i+1][j][k];
sw2=sw[i+1][j][k];
sg2=sg[i+1][j][k];
pcow2=interp(sat,pcowt,msat,sw2,pcow2);
pcgo2=interp(sat,pcgot,msat,sg2,pcgo2);
ro2=(rhosco + rso2*rhoscg)/bo[i+1][j][k];
rw2=(rhoscw + rsw2*rhoscg)/bw[i+1][j][k];
rg2=rhoscg/bg[i+1][j][k];
fact=-d288*(el[i+1][j][k]-el[i][j][k]);
gow2=(ro2+ro)*fact;
gww2=(rw2+rw)*fact + pcow-pcow2;
ggw2=(rg2+rg)*fact + pcgo2-pcgo;
p22=p2-pp;
ho2=p22+gow2;
hw2=p22+gww2;
hg2=p22+ggw2;
if(ho2>=0.)kro2=interp(sat,krot,msat,so2,kro2);
if(ho2<0.) kro2=interp(sat,krot,msat,sso,kro2);
if(hw2>=0.)krw2=interp(sat,krwt,msat,sw2,krw2);
if(hw2<0.) krw2=interp(sat,krwt,msat,ssw,krw2);
if(hg2>=0.)krg2=interp(sat,krgt,msat,sg2,krg2);
if(hg2<0.) krg2=interp(sat,krgt,msat,ssg,krg2);
mo2=4.0*kro2/((bo[i+1][j][k]+bo[i][j][k]) * (muo2+muo));
mw2=4.0*krw2/((bw[i+1][j][k]+bw[i][j][k]) * (muw2+muw));
mg2=4.0*krg2/((bg[i+1][j][k]+bg[i][j][k]) * (mug2+mug));
M125:
aoe=tx[i+1][j][k]*mo2;
awe=tx[i+1][j][k]*mw2;
age=tx[i+1][j][k]*mg2;
if(j==1)goto M135;
p3=p[i][j-1][k];
bpt=pbot[i][j-1][k];
rso3=intpvt(bpt,rslope,pot,rsot,mpot,p3,rso3);
muo3=intpvt(bpt,vslope,pot,muot,mpot,p3,muo3);
rsw3=interp(pwt,rswt,mpwt,p3,rsw3);
muw3=interp(pwt,muwt,mpwt,p3,muw3);
mug3=interp(pgt,mugt,mpgt,p3,mug3);
so3=so[i][j-1][k];
sw3=sw[i][j-1][k];
sg3=sg[i][j-1][k];
pcow3=interp(sat,pcowt,msat,sw3,pcow3);
pcgo3=interp(sat,pcgot,msat,sg3,pcgo3);
ro3=(rhosco + rso3*rhoscg)/bo[i][j-1][k];
rw3=(rhoscw + rsw3*rhoscg)/bw[i][j-1][k];
rg3=rhoscg/bg[i][j-1][k];
fact=-d288*(el[i][j-1][k]-el[i][j][k]);
gow3=(ro3+ro)*fact;
gww3=(rw3+rw)*fact + pcow-pcow3;
ggw3=(rg3+rg)*fact + pcgo3-pcgo;
p33=p3-pp;
ho3=p33+gow3;
hw3=p33+gww3;
hg3=p33+ggw3;
if(ho3>=0.)kro3=interp(sat,krot,msat,so3,kro3);
if(ho3<0.) kro3=interp(sat,krot,msat,sso,kro3);
if(hw3>=0.)krw3=interp(sat,krwt,msat,sw3,krw3);
if(hw3<0.) krw3=interp(sat,krwt,msat,ssw,krw3);
if(hg3>=0.)krg3=interp(sat,krgt,msat,sg3,krg3);
if(hg3<0.) krg3=interp(sat,krgt,msat,ssg,krg3);
mo3=4.0*kro3/((bo[i][j-1][k]+bo[i][j][k]) * (muo3+muo));
mw3=4.0*krw3/((bw[i][j-1][k]+bw[i][j][k]) * (muw3+muw));
mg3=4.0*krg3/((bg[i][j-1][k]+bg[i][j][k]) * (mug3+mug));
M135:
aos=ty[i][j][k]*mo3;
aws=ty[i][j][k]*mw3;
ags=ty[i][j][k]*mg3;
if(j==jj)goto M140;
p4=p[i][j+1][k];
bpt=pbot[i][j+1][k];
rso4=intpvt(bpt,rslope,pot,rsot,mpot,p4,rso4);
muo4=intpvt(bpt,vslope,pot,muot,mpot,p4,muo4);
rsw4=interp(pwt,rswt,mpwt,p4,rsw4);
muw4=interp(pwt,muwt,mpwt,p4,muw4);
mug4=interp(pgt,mugt,mpgt,p4,mug4);
so4=so[i][j+1][k];
sw4=sw[i][j+1][k];
sg4=sg[i][j+1][k];
pcow4=interp(sat,pcowt,msat,sw4,pcow4);
pcgo4=interp(sat,pcgot,msat,sg4,pcgo4);
ro4=(rhosco + rso4*rhoscg)/bo[i][j+1][k];
rw4=(rhoscw + rsw4*rhoscg)/bw[i][j+1][k];
rg4=rhoscg/bg[i][j+1][k];
fact=-d288*(el[i][j+1][k]-el[i][j][k]);
gow4=(ro4+ro)*fact;
gww4=(rw4+rw)*fact + pcow-pcow4;
ggw4=(rg4+rg)*fact + pcgo4-pcgo;
p44=p4-pp;
ho4=p44+gow4;
hw4=p44+gww4;
hg4=p44+ggw4;
if(ho4>=0.)kro4=interp(sat,krot,msat,so4,kro4);
if(ho4<0.) kro4=interp(sat,krot,msat,sso,kro4);
if(hw4>=0.)krw4=interp(sat,krwt,msat,sw4,krw4);
if(hw4<0.) krw4=interp(sat,krwt,msat,ssw,krw4);
if(hg4>=0.)krg4=interp(sat,krgt,msat,sg4,krg4);
if(hg4<0.) krg4=interp(sat,krgt,msat,ssg,krg4);
mo4=4.0*kro4/((bo[i][j+1][k]+bo[i][j][k]) * (muo4+muo));
mw4=4.0*krw4/((bw[i][j+1][k]+bw[i][j][k]) * (muw4+muw));
mg4=4.0*krg4/((bg[i][j+1][k]+bg[i][j][k]) * (mug4+mug));
M140:
aon=ty[i][j+1][k]*mo4;
awn=ty[i][j+1][k]*mw4;
agn=ty[i][j+1][k]*mg4;
if(k==1)goto M145;
p5=p[i][j][k-1];
bpt=pbot[i][j][k-1];
rso5=intpvt(bpt,rslope,pot,rsot,mpot,p5,rso5);
muo5=intpvt(bpt,vslope,pot,muot,mpot,p5,muo5);
rsw5=interp(pwt,rswt,mpwt,p5,rsw5);
muw5=interp(pwt,muwt,mpwt,p5,muw5);
mug5=interp(pgt,mugt,mpgt,p5,mug5);
so5=so[i][j][k-1];
sw5=sw[i][j][k-1];
sg5=sg[i][j][k-1];
pcow5=interp(sat,pcowt,msat,sw5,pcow5);
pcgo5=interp(sat,pcgot,msat,sg5,pcgo5);
ro5=(rhosco + rso5*rhoscg)/bo[i][j][k-1];
rw5=(rhoscw + rsw5*rhoscg)/bw[i][j][k-1];
rg5=rhoscg/bg[i][j][k-1];
fact=-d288*(el[i][j][k-1]-el[i][j][k]);
gow5=(ro5+ro)*fact;
gww5=(rw5+rw)*fact + pcow-pcow5;
ggw5=(rg5+rg)*fact + pcgo5-pcgo;
p55=p5-pp;
ho5=p55+gow5;
hw5=p55+gww5;
hg5=p55+ggw5;
if(ho5>=0.)kro5=interp(sat,krot,msat,so5,kro5);
if(ho5<0.) kro5=interp(sat,krot,msat,sso,kro5);
if(hw5>=0.)krw5=interp(sat,krwt,msat,sw5,krw5);
if(hw5<0.) krw5=interp(sat,krwt,msat,ssw,krw5);
if(hg5>=0.)krg5=interp(sat,krgt,msat,sg5,krg5);
if(hg5<0.) krg5=interp(sat,krgt,msat,ssg,krg5);
mo5=4.0*kro5/((bo[i][j][k-1]+bo[i][j][k]) * (muo5+muo));
mw5=4.0*krw5/((bw[i][j][k-1]+bw[i][j][k]) * (muw5+muw));
mg5=4.0*krg5/((bg[i][j][k-1]+bg[i][j][k]) * (mug5+mug));
M145:
aot=tz[i][j][k]*mo5;
awt=tz[i][j][k]*mw5;
agt=tz[i][j][k]*mg5;
if(k==kk)goto M150;
p6=p[i][j][k+1];
bpt=pbot[i][j][k+1];
rso6=intpvt(bpt,rslope,pot,rsot,mpot,p6,rso6);
muo6=intpvt(bpt,vslope,pot,muot,mpot,p6,muo6);
rsw6=interp(pwt,rswt,mpwt,p6,rsw6);
muw6=interp(pwt,muwt,mpwt,p6,muw6);
mug6=interp(pgt,mugt,mpgt,p6,mug6);
so6=so[i][j][k+1];
sw6=sw[i][j][k+1];
sg6=sg[i][j][k+1];
pcow6=interp(sat,pcowt,msat,sw6,pcow6);
pcgo6=interp(sat,pcgot,msat,sg6,pcgo6);
ro6=(rhosco + rso6*rhoscg)/bo[i][j][k+1];
rw6=(rhoscw + rsw6*rhoscg)/bw[i][j][k+1];
rg6=rhoscg/bg[i][j][k+1];
fact=-d288*(el[i][j][k+1]-el[i][j][k]);
gow6=(ro6+ro)*fact;
gww6=(rw6+rw)*fact + pcow-pcow6;
ggw6=(rg6+rg)*fact + pcgo6-pcgo;
p66=p6-pp;
ho6=p66+gow6;
hw6=p66+gww6;
hg6=p66+ggw6;
if(ho6>=0.)kro6=interp(sat,krot,msat,so6,kro6);
if(ho6<0.) kro6=interp(sat,krot,msat,sso,kro6);
if(hw6>=0.)krw6=interp(sat,krwt,msat,sw6,krw6);
if(hw6<0.) krw6=interp(sat,krwt,msat,ssw,krw6);
if(hg6>=0.)krg6=interp(sat,krgt,msat,sg6,krg6);
if(hg6<0.) krg6=interp(sat,krgt,msat,ssg,krg6);
mo6=4.0*kro6/((bo[i][j][k+1]+bo[i][j][k]) * (muo6+muo));
mw6=4.0*krw6/((bw[i][j][k+1]+bw[i][j][k]) * (muw6+muw));
mg6=4.0*krg6/((bg[i][j][k+1]+bg[i][j][k]) * (mug6+mug));
M150:
aob=tz[i][j][k+1]*mo6;
awb=tz[i][j][k+1]*mw6;
agb=tz[i][j][k+1]*mg6;
rso1a=0.5*(rso1+rso);
rso2a=0.5*(rso2+rso);
rso3a=0.5*(rso3+rso);
rso4a=0.5*(rso4+rso);
rso5a=0.5*(rso5+rso);
rso6a=0.5*(rso6+rso);
rsw1a=0.5*(rsw1+rsw);
rsw2a=0.5*(rsw2+rsw);
rsw3a=0.5*(rsw3+rsw);
rsw4a=0.5*(rsw4+rsw);
rsw5a=0.5*(rsw5+rsw);
rsw6a=0.5*(rsw6+rsw);
ao1=aow*gow1;
ao2=aoe*gow2;
ao3=aos*gow3;
ao4=aon*gow4;
ao5=aot*gow5;
ao6=aob*gow6;
aw1=aww*gww1;
aw2=awe*gww2;
aw3=aws*gww3;
aw4=awn*gww4;
aw5=awt*gww5;
aw6=awb*gww6;
gowt[i][j][k]= ao1 + ao2 + ao3 + ao4 + ao5 + ao6;
gwwt[i][j][k]= aw1 + aw2 + aw3 + aw4 + aw5 + aw6;
ggwt[i][j][k]=agw*ggw1+age*ggw2+ags*ggw3+agn*ggw4+agt*ggw5+agb*ggw6+
rso1a*ao1+rso2a*ao2+rso3a*ao3+rso4a*ao4+rso5a*ao5+rso6a*ao6+
rsw1a*aw1+rsw2a*aw2+rsw3a*aw3+rsw4a*aw4+rsw5a*aw5+rsw6a*aw6;
qowg[i][j][k]=(bo[i][j][k]-bg[i][j][k]*rso)*(-gowt[i][j][k]+qo[i][j][k])+
(bw[i][j][k]-bg[i][j][k]*rsw)*(-gwwt[i][j][k]+qw[i][j][k])+
bg[i][j][k]*(-ggwt[i][j][k]+qg[i][j][k]);
aw[i][j][k]=(bo[i][j][k] + 0.5*bg[i][j][k]*(rso1-rso)) * aow +
(bw[i][j][k] + 0.5*bg[i][j][k]*(rsw1-rsw)) * aww +
bg[i][j][k]*agw;
ae[i][j][k]=(bo[i][j][k] + 0.5*bg[i][j][k]*(rso2-rso)) * aoe +
(bw[i][j][k] + 0.5*bg[i][j][k]*(rsw2-rsw)) * awe +
bg[i][j][k]*age;
as[i][j][k]=(bo[i][j][k] + 0.5*bg[i][j][k]*(rso3-rso)) * aos +
(bw[i][j][k] + 0.5*bg[i][j][k]*(rsw3-rsw)) * aws +
bg[i][j][k]*ags;
an[i][j][k]=(bo[i][j][k] + 0.5*bg[i][j][k]*(rso4-rso)) * aon +
(bw[i][j][k] + 0.5*bg[i][j][k]*(rsw4-rsw)) * awn +
bg[i][j][k]*agn;
at[i][j][k]=(bo[i][j][k] + 0.5*bg[i][j][k]*(rso5-rso)) * aot +
(bw[i][j][k] + 0.5*bg[i][j][k]*(rsw5-rsw)) * awt +
bg[i][j][k]*agt;
ab[i][j][k]=(bo[i][j][k] + 0.5*bg[i][j][k]*(rso6-rso)) * aob +
(bw[i][j][k] + 0.5*bg[i][j][k]*(rsw6-rsw)) * awb +
bg[i][j][k]*agb;
ow[i][j][k]=aow;
oe[i][j][k]=aoe;
os[i][j][k]=aos;
on[i][j][k]=aon;
ot[i][j][k]=aot;
ob[i][j][k]=aob;
ww[i][j][k]=aww;
we[i][j][k]=awe;
ws[i][j][k]=aws;
wn[i][j][k]=awn;
wt[i][j][k]=awt;
wb[i][j][k]=awb;
}
/**** calculate main diagonal and rhs vector***********/
for(k=1;k<=kk;k++)
for(j=1;j<=jj;j++)
for(i=1;i<=ii;i++)
{
sum[i][j][k]=aw[i][j][k]+ae[i][j][k]+as[i][j][k]+an[i][j][k]+
at[i][j][k]+ab[i][j][k];
gam[i][j][k]=vp[i][j][k]*ct[i][j][k]*div1;
e[i][j][k]=-sum[i][j][k] - gam[i][j][k];
b[i][j][k]= qowg[i][j][k] - gam[i][j][k]*p[i][j][k];
}
return;
}
/*----------- print inittial data---------*/
void CBOASTView::initprtp(void)
{
int i,j,k;
fprintf(fw,"THE-OIL-RESERVOIR-PRESSURE-DISTRIBUTION\n");
for(k=1;k<=kk;k++)
{
fprintf(fw," k=%d\n",k);
for(j=1;j<=jj;j++) {
for(i=1;i<=ii;i++)
{ fprintf(fw,"%8.2f",p[i][j][k]*0.0068948);
if(i%10==0) fprintf(fw,"\n");
}
fprintf(fw,"\n");
}
}
fprintf(fw,"THE-OIL-SATURATION-DISTRIBUTION\n");
for(k=1;k<=kk;k++)
{
fprintf(fw," k=%d\n",k);
for(j=1;j<=jj;j++) {
for(i=1;i<=ii;i++)
{ fprintf(fw,"%8.2f",so[i][j][k]);
if(i%10==0) fprintf(fw,"\n");
}
fprintf(fw,"\n");
}
}
fprintf(fw,"THE-WATER-SATURATION-DISTRIBUTION\n");
for(k=1;k<=kk;k++)
{
fprintf(fw," k=%d\n",k);
for(j=1;j<=jj;j++) {
for(i=1;i<=ii;i++)
{ fprintf(fw,"%8.2f",sw[i][j][k]);
if(i%10==0) fprintf(fw,"\n");
}
fprintf(fw,"\n");
}
}
fprintf(fw,"THE-GAS-SATURATION-DISTRIBUTION\n");
for(k=1;k<=kk;k++)
{
fprintf(fw," k=%d\n",k);
for(j=1;j<=jj;j++) {
for(i=1;i<=ii;i++)
{ fprintf(fw,"%8.2f",sg[i][j][k]);
if(i%10==0) fprintf(fw,"\n");
}
fprintf(fw,"\n");
}
}
fprintf(fw,"THE-RESER-BUBBLE0POINT-PRESSURE-DISTRIBUTION\n");
for(k=1;k<=kk;k++)
{
fprintf(fw," k=%d\n",k);
for(j=1;j<=jj;j++) {
for(i=1;i<=ii;i++)
{ fprintf(fw,"%8.2f",pbot[i][j][k]*0.0068948);
if(i%10==0) fprintf(fw,"\n");
}
fprintf(fw,"\n");
}
}
return;
}
/*----------solution method ---------------------*/
void CBOASTView::lsor()
{
float um[lll],azl[lll],bzl[lll],czl[lll],dzl[lll],uzl[lll];
/* float div=delt/delt0; */
int niter=0;
float dmax=1.0;
float rho1=0.0;
float theta=0.0;
float omega=(float)1.7;
float tol=(float)0.10;
float tol1=0.00;
int miter=250;
int nx=ii,ny=jj,nz=kk;
int jm,jp,km,kp,i,j,k;
float tw,dmax0,theta0,gslsor,dm,arg;
float delta,om;
M11:
tw=1.0-omega;
dmax0=dmax;
theta0=theta;
if(niter>=miter) {
fprintf(fw,"conversion not reached . %5d %10.4 %12.5\n",niter,tol,dmax);
return; }
niter=niter+1;
dmax=0.0;
for(k=1;k<=nz;k++)
for(j=1;j<=ny;j++)
{
for(i=1;i<=nx;i++)
{
um[i]=p[i][j][k];
azl[i]=aw[i][j][k];
bzl[i]=e[i][j][k];
czl[i]=ae[i][j][k];
dzl[i]=b[i][j][k];
if(ny==1) goto M14;
jm=j-1;
jp=j+1;
if(j==1) jm=1;
if(j==ny) jp=ny;
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -