skin.cpp

石油行业油井损害软件,有详细的油井损害模型,对石油建模有重要帮助

#include "stdafx.h"
#include "储层损害矿场评价软件.h"

#include "SkinDoc.h"
#include "math.h"

static void fei(double, double&);

void CMyDoc::deskin()
{
	double dk, sqkdkv, thew;
	double delta,s2b3,s2b2,s2b1,x1p2,x2,x1,x,sdp,sg,sp;
	double swb,sv,hdd,rpd,drw,b,a1,a2,b1,b2,c1,c2,alph,rdp,dpk,a,sh;
	double sumt,tan1;
	int i;
	if(md!=0) //如果非达西渗流表皮系数由用户输入
		d=(1.e-9)*belt*rhls*hk/(542.86721*h*vis*rw);
	snd=d*qc; //计算非达西渗流表皮系数
	if(s>0.0) dk=0.5*hk;
	if(s<0.0) dk=2.*hk; 
	sqkdkv=sqrt(hk/vk);
	spt=(h/hp-1.)*(log(h*sqkdkv/rw)-2.0); //计算部分打开地层的表皮系数
	tan1=sin(wsl)/cos(wsl)/sqkdkv;
	thew=atan(tan1)*180./3.1415926; //计算井斜角
	sws=-pow((thew/56.),1.865)*log10(0.01*h*sqkdkv/rw);
	sws=sws-pow((thew/41.),2.06); //计算井斜表皮系数
	sca=0.5*log(31.62/ca);//计算地层边界形状表皮系数
	scp=(h/hp)*(1./rkfpc-1.)*log(rb/rw);//计算相变表皮系数
	if(npf == 0) // 裸眼井
		spt=0.0;
	else //射孔井
	{
		if(lphi == 0 || lphi == 360)
		{
			a1=-2.091;
			a2=0.0453;
			b1=5.1313;
			b2=1.8672;
			c1=0.16;
			c2=2.675;
			alph=0.25;
		}
		else
		{
			if(lphi == 180) 
			{
				a1=-2.025;
				a2=0.0943;
				b1=3.0373;
				b2=1.8115;
				c1=0.026;
				c2=4.532;
				alph=0.5;
			}
			else
			{
				if(lphi == 120) 
				{
					a1=-2.018;
					a2=0.0634;
					b1=1.6136;
					b2=1.777;
					c1=6.6e-3;
					c2=5.32;
					alph=0.648;
				}
				else
				{
					if(lphi == 90) 
					{
						a1=-1.905;
						a2=0.1038;
						b1=1.5674;
						b2=1.6953;
						c1=1.9e-3;
						c2=6.155;
						alph=0.726;
					}
					else
					{
						if(lphi == 60)
						{
							a1=-1.898;
							a2=0.1023;
							b1=1.3654;
							b2=1.649;
							c1=3.0e-4;
							c2=7.509;
							alph=0.813;
						}
						else
						{
							if(lphi == 45) 
							{
								a1=-1.788;
								a2=0.2398;
								b1=1.1915;
								b2=1.6392;
								c1=4.6e-5;
								c2=8.791;
								alph=0.86;
							}
						}
					}
				}
			}
		}
		rdp=0.0125+rp; //油井压实带半径
		dpk=0.175*dk;  //压实带渗透率
		a=a1*log10(rp*(1.+1./sqkdkv)/h)+a2;
		b=b1*rp*(1.+1./sqkdkv)/h+b2;
		if(lphi == 0)
			drw=0.25*flp;
		else
			drw=alph*(rw+flp); //有效的井筒半径
		rpd=rp*(1.+1./sqkdkv)/h;
		hdd=h*sqkdkv/(flp*(n));
		sh=log(rw/drw);
		sv=pow(10.,a)*pow(rpd,b)*pow(hdd,(b-1.));
		swb=c1*exp(c2*rw/(flp+rw));
		sp=sh+sv+swb;
		sg=2.*hk*h*flp/(gk*rp*rp*(n));//射孔充填线性流拟表皮系数
		sdp=2.0*hp*(hk/dpk-hk/dk)*log(rdp/rp)/(flp*(n));//射孔压实带拟表皮系数
		spf=sp+sg+sdp;//射孔表皮系数
	}

    if(n1b == 0)
		s1b=0.;
	else
	{
		x=phi*vis*ct*dwtob1*dwtob1/3.6/hk;
		fei(x,s1b);
		if(n1b == 1)
			s1b=0.5*s1b;
		else if(n1b == 2)
			s1b=-0.5*s1b;
		else
			s1b=0.5*(hk-hkout)/(hk+hkout)*s1b;
	}
    if(n2b == 0) 
		s2b=0.;
	else
	{
		x1=phi*vis*ct*dwtob1*dwtob1/3.6/hk;
		x2=phi*vis*ct*dwtob2*dwtob2/3.6/hk;
		x1p2=x1+x2;
		fei(x1,s2b1);
		fei(x2,s2b2);
		fei(x1p2,s2b3);
		if(n2b == 1) 
			s2b=0.5*(s2b1+s2b2+s2b3);
		else  if(n2b == 2) 
			s2b=-0.5*(s2b1+s2b2+s2b3);
		else
			s2b=0.5*(hk-hkout)/(hk+hkout)*(s2b1+s2b2+
			      (hk-3.*hkout)/(hk+3.*hkout)*s2b3);
	}
    if(nan == 0) //如果各向同性
	san=0.;
	else
	{
		delta=pow((hky/hkx),0.25);
		san=-log(0.5*(delta+1./delta)); //各向异性表皮系数
	}
    if(nvq == 0)  //如果不考虑产量变化
		svq=0.;  //变产量表皮系数为0
	else
	{
		sumt=0.0;
		for(i=1; i<=nq; i++)
			sumt=sumt+(q[i]-q[i-1])/q[nq]*log(1.+1./(tp-t[i-1]));
		svq=0.5*(log((tp+1.)/tp)-sumt); //变产量表皮系数
	}
	sym=(86.4*3.1416*h*hk*alf)*pow((pi-pwf),2)/(qc*b*vis);  ////压力敏感表皮系数;
    sd=s-spt-sws-spf-snd-scp-sca-s1b-s2b-san-svq-sym; //地层损害表皮系数
}

void CMyDoc::qvsds()
{
	double ds0;
	ds0=0.025*sd;
	for(int i=0; i<40; i++)		// 40是任意给定的
	{
		ds[i]=(i+1)*ds0;
		dqr[i]=0.87*hnslop*ds[i]/(pi-pwf-0.87*hnslop*ds[i]);
	}
	dpr0=0.87*hnslop*sd/(pi-pwf-0.87*hnslop*sd);
}
 
void CMyDoc::calrs()
{
	double a0,alnkkw,x10,x1,x20,x2,f,df,fddf,ve1,ve2;
	double sdn,rkdkw;
	sdn=sd*hp/h;

       if(sdn>0.&&sdn<=1.)   rkdkw=5.;
       if(sdn>1.&&sdn<=3.)   rkdkw=5.5;
      if(sdn>3.&&sdn<=5.)   rkdkw=6.5 ;
       if(sdn>5.&&sdn<=8.)   rkdkw=8. ;
       if(sdn>8.&&sdn<=10.)  rkdkw=10.;
       if(sdn>10.&&sdn<=15.) rkdkw=15. ;
       if(sdn>15.&&sdn<=25.) rkdkw=18.;
       if(sdn>25.) rkdkw=20.;

a10:   a0=12.5;
a50:   a0=a0-0.5;
       rs=a0*rw;
       alnkkw=log(rkdkw);
       if(rs <= rw) return;
a100:   x10=rw/(rs-rw);
       x1=x10*alnkkw;
       x20=rs/(rs-rw);
       x2=x20*alnkkw;
       if(rs <= rw) goto a50;
       fei(x1,ve1);
       fei(x2,ve2);
       f=ve1-ve2-(log(rs/rw)+sdn)*(pow(rkdkw,(-x20)));
       df=(1.+x1*(log(rs/rw)+sdn))*(pow(rkdkw,(-x20)));
       df=(pow(rkdkw,(-x10))-df)/(rs-rw);
		fddf=f/df;
       rs=rs-fddf;
       if(fabs(f) > 0.001 || fabs(fddf) > 0.001) goto a100;
       fls=rs-rw;
	if(fls>=0.01&&fls<=1.5) goto a120;
	if(fls>1.5) rkdkw=rkdkw+0.5;
	if(fls<0.01) rkdkw=rkdkw-0.5;
	goto a10;
a120:  ; 
}

void fei(double x, double& ei)
{
	double d1,d2,d3,d4,d5,f1,f2,f3,g1,g2,g3,ei1,ei2;

	d1=0.99999193;
	d2=-0.2499105;
	d3=0.05519968;
	d4=-0.00976004;
	d5=0.00107857;
	f1=21.0996530827;
	f2=25.6329561486;
	f3=9.5733223454;
	g1=8.6347608925;
	g2=18.059016973;
	g3=8.573328401;
    if(x <= 1.) 
       ei=-(log(x)+0.57721566)+x*(d1+x*(d2+x*(d3+x*(d4+d5*x))));
	else
	{
		ei1=x*(g1+x*(g2+x*(g3+x)))+0.2677737343;
		ei2=x*(f1+x*(f2+x*(f3+x)))+3.9584969228;
		ei=exp(-x)/x*ei1/ei2;
	}
}