📄 boastview.cpp
字号:
boder=interp(pot,bopt,mpot,pp,boder);
rsoder=interp(pot,rsopt,mpot,pp,rsoder);
bwder=interp(pwt,bwpt,mpwt,pp,bwder);
rswder=interp(pwt,rswpt,mpwt,pp,rswder);
bgder=interp(pgt,bgpt,mpgt,pp,bgder);
if(pp>pbot[i][j][k])boder=bslope;
if(pp>pbot[i][j][k])rsoder=rslope;
co=-(boder-bg[i][j][k]*rsoder)/bo[i][j][k];
cw=-(bwder-bg[i][j][k]*rswder)/bw[i][j][k];
cg=-bgder/bg[i][j][k];
cr=interp(pgt,crt,mpgt,pp,cr);
ct[i][j][k]=cr+co*so[i][j][k]+cw*sw[i][j][k]+cg*sg[i][j][k];
}
/* auto. time step control calc. of pressure and sat. maxima. */
ppm=0.0;
som=0.0;
swm=0.0;
sgm=0.0;
for(k=1;k<=kk;k++)
for(j=1;j<=jj;j++)
for(i=1;i<=ii;i++)
{
dpo=p[i][j][k]-pn[i][j][k];
dso=so[i][j][k]-son[i][j][k];
dsw=sw[i][j][k]-swn[i][j][k];
dsg=sg[i][j][k]-sgn[i][j][k];
if(fabs(dpo)>=fabs(ppm)) ppm=dpo;
if(fabs(dso)>=fabs(som)) som=dso;
if(fabs(dsw)>=fabs(swm)) swm=dsw;
if(fabs(dsg)>=fabs(sgm)) sgm=dsg;
}
dpmc=fabs(ppm);
dsmc=fabs(som);
if(dsmc<fabs(sgm))dsmc=fabs(sgm);
if(dsmc<fabs(swm))dsmc=fabs(swm);
/**** repeat time step *****/
if(dsmc<dsmax && dpmc<dpmax) goto M402;
if(delt<=dtmin || fact2>=1.0)goto M402;
itflag=itflag+1;
delt=delt*fact2;
if(delt<dtmin) delt=dtmin;
ft=eti+delt;
if(ft>ftmax)delt=ftmax-eti;
/**** reset variables. ****/
for(k=1;k<=kk;k++)
for(j=1;j<=jj;j++)
for(i=1;i<=ii;i++)
{
p[i][j][k]=pn[i][j][k];
so[i][j][k]=son[i][j][k];
sw[i][j][k]=swn[i][j][k];
sg[i][j][k]=sgn[i][j][k];
}
goto Mark1060;
M402:
/**** undersaturated grid block saturation caculation. *****/
for(k=1;k<=kk;k++)
for(j=1;j<=jj;j++)
for(i=1;i<=ii;i++)
{
if(p[i][j][k]>pn[i][j][k]) continue;
if(p[i][j][k]<pbot[i][j][k]) continue;
ip=i+1;
im=i-1;
jp=j+1;
jm=j-1;
kp=k+1;
km=k-1;
if(ip>ii) goto M412;
if(sgn[ip][j][k]>0.0001) continue;
M412: if(im<1) goto M414;
if(sgn[im][j][k]>0.0001) continue;
M414: if(jp>jj) goto M416;
if(sgn[i][jp][k]>0.0001) continue;
M416: if(jm<1) goto M418;
if(sgn[i][jm][k]>0.0001) continue;
M418: if(kp>kk) goto M420;
if(sgn[i][j][kp]>0.0001) continue;
M420: if(km<1) goto M422;
if(sgn[i][j][km]>0.0001) continue;
M422: sg[i][j][k]=0.0;
so[i][j][k]=1.0-sw[i][j][k];
}
/**** repressurization algorithm ******/
if(ireprs==1) goto M5022;
for(k=1;k<=kk;k++)
for(j=1;j<=jj;j++)
for(i=1;i<=ii;i++)
{
if(sg[i][j][k]<=0.0001) continue;
pp=p[i][j][k];
if(p[i][j][k]>pbot[i][j][k]) pp=pbot[i][j][k];
bbo=interp(pot,bot,mpot,pp,bbo);
rso=interp(pot,rsot,mpot,pp,rso);
bbg=interp(pgt,bgt,mpgt,pp,bbg);
if(so[i][j][k]<=0.01) continue;
rsonew=rso+sg[i][j][k]*bbo/(so[i][j][k]*bbg);
pbonew=interp(rsot,pot,mpot,rsonew,pbonew);
pbot[i][j][k]=pbonew;
}
M5022:
/***** update old fluid volumes for material balance. *****/
stboi=stbo;
stbwi=stbw;
mcfgi=mcfgt;
/**** update new fluid volumes.*********/
stbo=scfo*d5615;
stbw=scfw*d5615;
mcfg=scfg*0.001;
mcfg1=scfg1*0.001;
mcfgt=mcfg+mcfg1;
/**** well report ( all rate & pressure data applicable this step )********/
flag=0.0;
M590:
if(iwlrep==1)
fprintf(fw,"\n\n-----------------------WELL REPORT TIME=%8.1f DAYS FROM BEGINING----------------------------\n",ft);
ij=0;
tor=0.0;
tgr=0.0;
twr=0.0;
toc=0.0;
tgc=0.0;
twc=0.0;
for(j=1;j<=nvqn;j++)
{
gor=0.0;
wor=0.0;
iq1=iqn1[j];
iq2=iqn2[j];
iq3=iqn3[j];
ij=ij+1;
lay=iq3+(layer[j]-1);
for(k=iq3;k<=lay;k++)
{
qoo=qo[iq1][iq2][k]*d5615;
qww=qw[iq1][iq2][k]*d5615;
qgg=qg[iq1][iq2][k]*0.001;
cumo[j][k]=cumo[j][k]+qoo*delt*0.001;
cumw[j][k]=cumw[j][k]+qww*delt*0.001;
cumg[j][k]=cumg[j][k]+qgg*delt*0.001;
if(iwlrep==0) continue;
if(ij==1 && k==iq3)
fprintf(fw," wellname i j k calcpwfc realpwf pid qo qg qw gor wor cumo cumg cumw\n");
if(qoo!=0.)
{ gor=qgg*1000./qoo;
wor=qww/qoo;
}
fprintf(fw,"%10s%3d%3d%3d %5.2f %6.2f %7.3f %5.2f %6.2f %7.2f %6.2f %6.2f %6.2f %7.2f %8.2f\n",
wellid[j],iqn1[j],iqn2[j],k,pwfc[j][k]*0.0068948,pwf[j][k]*0.0068948,pid[j][k]*0.3048,
qoo*rhosco*0.01602*158.988/1000.,qgg*0.02832,qww*rhoscw*0.01602*158.988/1000.,
gor*0.1781,wor,cumo[j][k]*158.988*rhosco*0.01602/1000.,cumg[j][k]*0.02832,
cumw[j][k]*158.988*rhoscw*0.01602/1000.);
tor=tor+qoo;
tgr=tgr+qgg;
twr=twr+qww;
toc=toc+cumo[j][k];
tgc=tgc+cumg[j][k];
twc=twc+cumw[j][k];
}
}
if(iwlrep!=0){
fprintf(fw,"------------------------------------------------------------------------------------------------------\n");
fprintf(fw," TOTAL: %9.2f%8.1f%8.1f %8.0f%8.0f%8.0f\n",
tor*rhosco*0.01602*158.988/1000.,tgr*0.02832,twr*rhoscw*158.988*0.01602/1000.,
toc*158.988*rhosco*0.01602/1000.,tgc*0.02832,twc*rhoscw*158.988*0.01602/1000.);
}
/**** caculate material balance errors & average reservoir pressure****/
if(flag!=1.0)
{delt0=delt;
eti=eti+delt; }
matbal(stbo,stboi,stbw,stbwi,mcfg,mcfgi,delt0,resvol,n,d5615,mcfg1,mcfgt);
/****************** start add three array ***************/
Tarry[n]=ft;
Parry[n]=pavg;
Rarry[n]=wor/(1+wor);
/****************** end add ***********************/
if(flag==1.0) goto M2056;
if(wor>wormax) { fprintf(fw,"MAXIMUM WOR HAS BEEN EXCEEDED--PROG IS BEING TERMINATED\n");goto M1001;}
if(gor>gormax) { fprintf(fw,"MAXIMUM GOR HAS BEEN EXCEEDED--PROG IS BEING TERMINATED\n");goto M1001;}
if(pavg<pamin) { fprintf(fw,"MINIUMM PRE HAS BEEN EXCEEDED--PROG IS BEING TERMINATED\n");goto M1001;}
if(pavg>pamax) { fprintf(fw,"MAXIMUM PRE HAS BEEN EXCEEDED--PROG IS BEING TERMINATED\n");goto M1001;}
/**** summary report.***********/
ppj[pjnum]=pavg*0.0068948;
fwj[pjnum]=wor/(1+wor);
timej[pjnum]=ft;
pjnum=pjnum+1;
sprintf(ssn,"Ellape time is :%8.2f Day",ft);
pDC->TextOut(410,33,ssn);
if(isumry==0)goto M500;
M2056:
nlp=n+1;
prtps(nlp,ipmap,isomap,iswmap,isgmap,ipbmap,delt0,eti);
if(flag==1.0) {goto M1006;}
M500:
/**** update arrays.**************/
for(k=1;k<=kk;k++)
for(j=1;j<=jj;j++)
for(i=1;i<=ii;i++)
{
qo[i][j][k]=0.0;
qw[i][j][k]=0.0;
qg[i][j][k]=0.0;
pn[i][j][k]=p[i][j][k];
son[i][j][k]=so[i][j][k];
swn[i][j][k]=sw[i][j][k];
sgn[i][j][k]=sg[i][j][k];
}
continue;
M1001:
flag=1.0;
iwlrep=1;
goto M590;
} /*------- n cycle is end ----*/
M1006:
AfxMessageBox("计算结束");
fclose(fr);
fclose(fw);
fclose(oko);
fclose(okw);
fclose(okp);
ReleaseDC(pDC);
return;
}
void CBOASTView::initvary()
{
int i,j,k;
for(k=0;k<nnn;k++)
for(j=0;j<mmm;j++)
for(i=0;i<lll;i++)
{ aw[i][j][k]=0.;
as[i][j][k]=0.;
at[i][j][k]=0.;
ae[i][j][k]=0.;
an[i][j][k]=0.;
ab[i][j][k]=0.;
}
for(k=0;k<nnn;k++)
for(j=0;j<mmm;j++)
for(i=0;i<lll+1;i++)
{ tx[i][j][k]=0.;
ow[i][j][k]=0.;
oe[i][j][k]=0.;
ww[i][j][k]=0.;
we[i][j][k]=0.;
}
for(k=0;k<nnn;k++)
for(j=0;j<mmm+1;j++)
for(i=0;i<lll;i++)
{ ty[i][j][k]=0.;
os[i][j][k]=0.;
on[i][j][k]=0.;
ws[i][j][k]=0.;
wn[i][j][k]=0.;
}
for(k=0;k<nnn+1;k++)
for(j=0;j<mmm;j++)
for(i=0;i<lll;i++)
{ tz[i][j][k]=0.;
ot[i][j][k]=0.;
ob[i][j][k]=0.;
wt[i][j][k]=0.;
wb[i][j][k]=0.;
}
for(j=0;j<nnn;j++)
for(i=0;i<maxwell;i++)
{
cumo[i][j]=0.;
cumw[i][j]=0.;
cumg[i][j]=0.;
}
return;
}
/*-----------to read dz dx dy depth ---------------*/
void CBOASTView::grid1()
{
int i,j,k;
char words[maxword];
float length=(float)0.3048,rath;
float del,sum[lll][mmm];
fscanf(fr,"%s",words);
fprintf(fw,"%s\n",words);
fscanf(fr,"%d%d%d%f",&ii,&jj,&kk,&rath);
fprintf(fw,"%5d%5d%5d%5.2f\n",ii,jj,kk,rath);
fscanf(fr,"%s",words);
fprintf(fw,"%s\n",words);
for(i=1;i<=ii;i++) {
fscanf(fr,"%f",&dx[i][1][1]);
fprintf(fw,"%8.2f",dx[i][1][1]);
if(i%10==0) fprintf(fw,"\n");
dx[i][1][1]=dx[i][1][1]/length;
}
fprintf(fw,"\n");
for(k=1;k<=kk;k++)
for(j=1;j<=jj;j++)
for(i=1;i<=ii;i++)
dx[i][j][k]=dx[i][1][1];
for(j=1;j<=jj;j++) {
fscanf(fr,"%f",&dy[1][j][1]);
fprintf(fw,"%8.2f",dy[1][j][1]);
if(j%10==0) fprintf(fw,"\n");
dy[1][j][1]=dy[1][j][1]/length;
}
fprintf(fw,"\n");
for(k=1;k<=kk;k++)
for(j=1;j<=jj;j++)
for(i=1;i<=ii;i++)
dy[i][j][k]=dy[1][j][1];
fscanf(fr,"%s",words);
fprintf(fw,"%s\n",words);
for(k=1;k<=kk;k++)
for(j=1;j<=jj;j++) {
for(i=1;i<=ii;i++)
{ fscanf(fr,"%f",&dz[i][j][k]);
fprintf(fw,"%8.2f",dz[i][j][k]);
if(i%10==0) fprintf(fw,"\n");
dz[i][j][k]=dz[i][j][k]*rath/length;
}
fprintf(fw,"\n"); }
fscanf(fr,"%s",words);
fprintf(fw,"%s\n",words);
for(j=1;j<=jj;j++) {
for(i=1;i<=ii;i++)
{ fscanf(fr,"%f",&varel[i][j]);
fprintf(fw,"%8.2f",varel[i][j]);
if(i%10==0) fprintf(fw,"\n");
varel[i][j]=varel[i][j]/length;
}
fprintf(fw,"\n"); }
for(j=1;j<=jj;j++)
for(i=1;i<=ii;i++)
sum[i][j]=0.;
for(k=1;k<=kk;k++)
for(j=1;j<=jj;j++)
for(i=1;i<=ii;i++)
{ del=sum[i][j]+dz[i][j][k]*0.50;
el[i][j][k]=varel[i][j]+del;
sum[i][j]=dz[i][j][k]+sum[i][j];
}
return;
}
void CBOASTView::parm1() /*--------------to read vp and kx,ky,kz-----------*/
{
int i,j,k;
char words[maxword];
fscanf(fr,"%s",words);
fprintf(fw,"%s\n",words);
for(k=1;k<=kk;k++)
for(j=1;j<=jj;j++) {
for(i=1;i<=ii;i++)
{ fscanf(fr,"%f",&vp[i][j][k]);
phi[i][j][k]=vp[i][j][k];
fprintf(fw,"%8.2f",vp[i][j][k]);
if(i%10==0) fprintf(fw,"\n");
}
fprintf(fw,"\n"); }
for(k=1;k<=kk;k++)
for(j=1;j<=jj;j++) {
for(i=1;i<=ii;i++)
{ fscanf(fr,"%f",&kx[i][j][k]);
fprintf(fw,"%8.2f",kx[i][j][k]);
if(i%10==0) fprintf(fw,"\n");
}
fprintf(fw,"\n"); }
for(k=1;k<=kk;k++)
for(j=1;j<=jj;j++) {
for(i=1;i<=ii;i++)
{ fscanf(fr,"%f",&ky[i][j][k]);
fprintf(fw,"%8.2f",ky[i][j][k]);
if(i%10==0) fprintf(fw,"\n");
}
fprintf(fw,"\n"); }
for(k=1;k<=kk;k++)
{ fscanf(fr,"%f",&kz[1][1][k]);
fprintf(fw,"%8.4f\n",kz[1][1][k]);
}
for(k=1;k<=kk;k++)
for(j=1;j<=jj;j++)
for(i=1;i<=ii;i++)
kz[i][j][k]=kz[1][1][k];
return;
}
void CBOASTView::tran1() /*-caculate interblock transmissibilities--*/
{
int i,j,k;
float facx,facy,facz;
char ssn[100];
for(k=1;k<=kk;k++)
for(j=1;j<=jj;j++)
for(i=1;i<=ii;i++)
{
facx=1.0;
facy=1.0;
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -