iprcal.cpp

这是本人两年前兼职为某个公司做的石油钻进设计软件

	return 0;
}






//=====================二项式压力方法计算压力与产量的关系======================
//fResPressure——地层压力(MPa)
//fA——二项式产能方程系数(MPa.d/万方)
//fB——二项式产能方程系数(MPa.d^2/万方^2)
//pfPwf——井底流压(MPa)
//pfRate——流量(万方)
//iNum——计算的压力、流量数据点数
int IPRCurve_AB_Pressure(float fResPressure, float fA, float fB, float* pfPwf, float* pfRate, int iNum)
{
	//if(fResPressure<0.2) return ERROR_RESPRESSURE;
	//if(fB<1e-10) return ERROR_B;
	//if(fA<1e-10) return ERROR_A;

	double  fQAof  = (-fA+sqrt(fA*fA+4*fB*fResPressure))/(2*fB);
	for(int i=0; i<iNum; i++)
	{
		pfRate[i] = float(fQAof/(iNum-1)*i);
		pfPwf[i]  = fResPressure - fA*pfRate[i] - fB*pfRate[i]*pfRate[i];
	}
	return ERROR_NON;
}


//====================二项式压力平方方法计算压力与产量的关系===================
//fResPressure——地层压力(MPa)
//fA——二项式产能方程系数(MPa^2.d/万方)
//fB——二项式产能方程系数(MPa^2.d^2/万方^2)
//pfPwf——井底流压(MPa)
//pfRate——流量(万方)
//iNum——计算的压力、流量数据点数
int IPRCurve_AB_SquaredPressure(float fResPressure, float fA, float fB, float* pfPwf, float* pfRate, int iNum)
{
	//if(fResPressure<0.2) return ERROR_RESPRESSURE;
	//if(fB<1e-10) return ERROR_B;
	//if(fA<1e-10) return ERROR_A;

	double  fTemp, fSqdPressure = fResPressure*fResPressure;
	double  fQAof  = (-fA+sqrt(fA*fA+4*fB*fSqdPressure))/(2*fB);
	for(int i=0; i<iNum; i++)
	{
		pfRate[i] = float(fQAof/(iNum-1)*i);
		fTemp     = fSqdPressure - fA*pfRate[i] - fB*pfRate[i]*pfRate[i];
		if(fTemp>0)
			pfPwf[i] = (float)sqrt(fTemp);
		else
			pfPwf[i] = 0;
	}
	return ERROR_NON;
}




//====================二项式拟压力方法计算压力与产量的关系=====================
//fResPressure——地层压力(MPa)
//fA——二项式产能方程系数(MPa^2/mPa.s . d/万方)
//fB——二项式产能方程系数(MPa^2/mPa.s . d^2/万方^2)
//pfPwf——井底流压(MPa)
//pfRate——流量(万方)
//iNum——计算的压力、流量数据点数
//pfPMp_Pressure——压力拟压力序列中的压力值
//pfPMp_Mp——压力拟压力序列中的拟压力值
//pfPMp_Num——压力拟压力序列数据点数
int IPRCurve_AB_PseudoPressure(float fResPressure, float fA, float fB, float* pfPwf, float* pfRate, int iNum, float* pfPMp_Pressure, float* pfPMp_Mp, int pfPMp_Num)
{
	//if(fResPressure<0.2) return ERROR_RESPRESSURE;
	//if(fB<1e-10) return ERROR_B;
	//if(fA<1e-10) return ERROR_A;

	double  fTemp, fMpressure = InsertValue_Linear( pfPMp_Num, pfPMp_Pressure, pfPMp_Mp, fResPressure );
	double  fQAof  = (-fA+sqrt(fA*fA+4*fB*fMpressure))/(2*fB);
	for(int i=0; i<iNum; i++)
	{
		pfRate[i] = float(fQAof/(iNum-1)*i);
		fTemp     = fMpressure - fA*pfRate[i] - fB*pfRate[i]*pfRate[i];
		pfPwf[i]  = InsertValue_Linear( pfPMp_Num, pfPMp_Mp, pfPMp_Pressure, float(fTemp) );
	}
	return ERROR_NON;
}






//=====================指数式压力方法计算压力与产量的关系======================
//fResPressure——地层压力(MPa)
//fC——指数式产能方程系数C
//fN——指数式产能方程系数n
//pfPwf——井底流压(MPa)
//pfRate——流量(万方)
//iNum——计算的压力、流量数据点数
int IPRCurve_CN_Pressure(float fResPressure, float fC, float fN, float* pfPwf, float* pfRate, int iNum)
{
	//if(fResPressure<0.2) return ERROR_RESPRESSURE;
	//if(fC<1e-10) return ERROR_C;
	//if(fN<1e-10) return ERROR_N;

	double  fQAof  = fC*pow(fResPressure,fN);
	for(int i=0; i<iNum; i++)
	{
		pfRate[i] = float(fQAof/(iNum-1)*i);
		pfPwf[i]  = float(fResPressure - pow(10, (log10(pfRate[i])-log10(fC))/fN));
	}
	return ERROR_NON;
}



//====================指数式压力平方方法计算压力与产量的关系===================
//fResPressure——地层压力(MPa)
//fC——指数式产能方程系数C
//fN——指数式产能方程系数n
//pfPwf——井底流压(MPa)
//pfRate——流量(万方)
//iNum——计算的压力、流量数据点数
int IPRCurve_CN_SquaredPressure(float fResPressure, float fC, float fN, float* pfPwf, float* pfRate, int iNum)
{
	//if(fResPressure<0.2) return ERROR_RESPRESSURE;
	//if(fC<1e-10) return ERROR_C;
	//if(fN<1e-10) return ERROR_N;

	double  fTemp, fSqdPressure = fResPressure*fResPressure;
	double  fQAof  = fC*pow(fSqdPressure,fN);
	for(int i=0; i<iNum; i++)
	{
		pfRate[i] = float(fQAof/(iNum-1)*i);
		fTemp     = fSqdPressure - pow(10, (log10(pfRate[i])-log10(fC))/fN);
		if(fTemp>0)
			pfPwf[i] = (float)sqrt(fTemp);
		else
			pfPwf[i] = 0;
	}
	return ERROR_NON;
}



//====================指数式拟压力方法计算压力与产量的关系=====================
//fResPressure——地层压力(MPa)
//fC——指数式产能方程系数C
//fN——指数式产能方程系数n
//pfPwf——井底流压(MPa)
//pfRate——流量(万方)
//iNum——计算的压力、流量数据点数
//pfPMp_Pressure——压力拟压力序列中的压力值
//pfPMp_Mp——压力拟压力序列中的拟压力值
//pfPMp_Num——压力拟压力序列数据点数
int IPRCurve_CN_PseudoPressure(float fResPressure, float fC, float fN, float* pfPwf, float* pfRate, int iNum, float* pfPMp_Pressure, float* pfPMp_Mp, int pfPMp_Num)
{
	//if(fResPressure<0.2) return ERROR_RESPRESSURE;
	//if(fC<1e-10) return ERROR_C;
	//if(fN<1e-10) return ERROR_N;

	double  fTemp, fMpressure = InsertValue_Linear( pfPMp_Num, pfPMp_Pressure, pfPMp_Mp, fResPressure );
	double  fQAof  = fC*pow(fMpressure,fN);
	for(int i=0; i<iNum; i++)
	{
		pfRate[i] = float(fQAof/(iNum-1)*i);
		fTemp     = fMpressure - pow(10, (log10(pfRate[i])-log10(fC))/fN);
		pfPwf[i]  = InsertValue_Linear( pfPMp_Num, pfPMp_Mp, pfPMp_Pressure, float(fTemp) );
	}
	return ERROR_NON;
}





//==================根据地层参数计算压力方法二项式系数A========================
//fU——粘度(mPa.s)
//fZ——气体偏差因子
//fT——气层温度(开氏度)
//fP——气层压力(MPa)
//fK——气层渗透率(mD)
//fh——气层厚度(m)
//fRe——气层供给半径(m)
//frw——井半径(m)
//fS——表皮系数
//计算压力、产量关系时，压力单位为MPa，产量单位为万方
float CalCoeff_A_Pressure(float fU, float fZ, float fT, float fP, float fK, float fh, float fRe, float frw, float fS)
{
	double fCoeff_A = 6.3667*fU*fZ*fT/fP/fK/fh;
	fCoeff_A *= (log(fRe/frw)+fS);
	return float(fCoeff_A);
}
//==================根据地层参数计算压力方法二项式系数B========================
//frg——气体比重
float CalCoeff_B_Pressure(float frg, float fZ, float fT, float fP, float fK, float fh, float frw)
{
	double fBeta    = 7.644*pow(10,10)/pow(fK,1.5);
	double fCoeff_B = 1.409*pow(10,-13)*fBeta*frg*fZ*fT/frw/fh/fh/fP;
	return float(fCoeff_B);
}





//==================根据地层参数计算压力平方方法二项式系数A====================
//fU——粘度(mPa.s)
//fZ——气体偏差因子
//fT——气层温度(开氏度)
//fK——气层渗透率(mD)
//fh——气层厚度(m)
//fRe——气层供给半径(m)
//frw——井半径(m)
//fS——表皮系数
//计算压力、产量关系时，压力单位为MPa，产量单位为万方
float CalCoeff_A_SquaredPressure(float fU, float fZ, float fT, float fK, float fh, float fRe, float frw, float fS)
{
	double fCoeff_A = 12.734*fU*fZ*fT/fK/fh;
	fCoeff_A *= (log(fRe/frw)+fS);
	return float(fCoeff_A);
}
//==================根据地层参数计算压力平方方法二项式系数B====================
//frg——气体比重
float CalCoeff_B_SquaredPressure(float frg, float fZ, float fT, float fK, float fh, float frw)
{
	double fBeta    = 7.644*pow(10,10)/pow(fK,1.5);
	double fCoeff_B = 2.825*pow(10,-13)*fBeta*frg*fZ*fT/frw/fh/fh;
	return float(fCoeff_B);
}





//==================根据地层参数计算拟压力方法二项式系数A======================
//fT——气层温度(开氏度)
//fK——气层渗透率(mD)
//fh——气层厚度(m)
//fRe——气层供给半径(m)
//frw——井半径(m)
//fS——表皮系数
//计算压力、产量关系时，压力单位为MPa，产量单位为万方
float CalCoeff_A_PseudoPressure(float fT, float fK, float fh, float fRe, float frw, float fS)
{
	double fCoeff_A = 12.734*fT/fK/fh;
	fCoeff_A *= (log(fRe/frw)+fS);
	return float(fCoeff_A);
}
//==================根据地层参数计算拟压力方法二项式系数B======================
//fU——粘度(mPa.s)
//frg——气体比重
float CalCoeff_B_PseudoPressure(float fU, float frg, float fT, float fK, float fh, float frw)
{
	double fBeta    = 7.644*pow(10,10)/pow(fK,1.5);
	double fCoeff_B = 2.825*pow(10,-13)*fBeta*frg*fT/fU/frw/fh/fh;
	return float(fCoeff_B);
}





//======================计算两相拟压力下地层中的压力分布=======================
//计算半径r(m)处的两相拟压力(MPa.kg/m^3/mPa.s)
//fr——半径(m)
//fRw——井有效半径(m)
//fRe——气层供给半径(m)
//fMpe——外边界处两相拟压力(MPa.kg/m^3/mPa.s)
//fMpwf——井底处两相拟压力(MPa.kg/m^3/mPa.s)
float CalTMpbyRadius(float fr, float fRw, float fRe, float fMpe, float fMpwf)
{
	double fTemp1 = log(fr/fRw) / log(fRe/fRw);
	double fTemp2 = fMpe - fMpwf;
	return float(fTemp1*fTemp2+fMpwf);
}




//=================根据地层参数计算两相拟压力方法二项式系数A===================
//fGamag——气体比重
//fK——气层渗透率(mD)
//fh——气层厚度(m)
//fRe——气层供给半径(m)
//frw——井半径(m)
//fS——表皮系数
float CalCoeff_A_TwoPseudoPressure(float fGamag, float fK, float fh, float fRe, float frw, float fS)
{
	double fCoeff_A = 1.842*10000*1.205*fGamag/fK/fh;
	fCoeff_A *= (log(fRe/frw)+fS);
	return float(fCoeff_A);
}

//========================计算两相拟压力方法二项式系数B========================
//计算两相拟压力方法二项式系数B的被积函数
//fr——半径(m)
//fPress——半径r处的压力(MPa)
//fK——气层渗透率(mD)
//iGCPvtNum——PVT数据点数
//pfP——计算PVT数据的压力(MPa)
//pflg——计算的气相质量分数
//pfUg——计算粘度(mPa.s)
//pfKroDivKrgPvt——根据PVT参数计算的Kro除Krg
//iSatNum——相对渗透率曲线数据点数
//pfSg——相对渗透率曲线中饱和度值
//pfKrg——相对渗透率曲线气相相对渗透率值
//pfKroDivKrgSat——根据相对渗透率曲线计算的Kro除Krg
float CalNonDarcyFun_TwoPseudoPressure(float  fr,      float fPress, float  fK,    int    iGCPvtNum, 
									   float* pfP,     float* pflg,  float* pfUg,  float* pfKroDivKrgPvt,
									   int    iSatNum, float* pfSg,  float* pfKrg, float* pfKroDivKrgSat)
{
	float  flg    = InsertValue_Linear( iGCPvtNum, pfP, pflg, fPress );
	float  fUg    = InsertValue_Linear( iGCPvtNum, pfP, pfUg, fPress );
	//if(fUg<0) 
	//fUg = 0.005;
	float  fSg    = GetSgByPress( iSatNum, pfKroDivKrgSat, pfSg, iGCPvtNum, pfKroDivKrgPvt, pfP, fPress );
	float  fKrg   = GetKrg( iSatNum, pfKrg, pfSg, fSg );
	double fBeta  = 7.644*pow(10,10)/pow(fK*fKrg,1.5);
	//double fBeta  = 7.644*pow(10,10)/pow(fK,1.5);
	double fTemp  = fBeta*fKrg*flg*flg/fUg/fr/fr;

	if( fTemp < 0 ) fTemp = 0;
	return float(fTemp);
}
//辛普森积分计算两相拟压力方法二项式系数B的积分部分
//fK——气层渗透率(mD)
//fa——井半径(m)
//fb——气层供给半径(m)
//fMpwf——井底处两相拟压力(MPa.kg/m^3/mPa.s)
//fMpe——外边界处两相拟压力(MPa.kg/m^3/mPa.s)
//iMpNum——两相拟压力数据点数
//pfMP_M——两相拟压力中的拟压力
//pfMP_P——两相拟压力中的压力
//iGCPvtNum——PVT数据点数
//pfP——计算PVT数据的压力(MPa)
//pflg——计算的气相质量分数
//pfUg——计算粘度(mPa.s)
//pfKroDivKrgPvt——根据PVT参数计算的Kro除Krg
//iSatNum——相对渗透率曲线数据点数
//pfSg——相对渗透率曲线中饱和度值
//pfKrg——相对渗透率曲线气相相对渗透率值
//pfKroDivKrgSat——根据相对渗透率曲线计算的Kro除Krg
float SimpsonIntegrate(float fK,      float fa,      float fb,      float fMpwf, float fMpe,
					   int iMpNum,    float* pfMP_M, float* pfMP_P, 
					   int iGCPvtNum, float* pfP,    float* pflg,   float* pfUg,  float* pfKroDivKrgPvt,
					   int iSatNum,   float* pfSg,   float* pfKrg,  float* pfKroDivKrgSat)
{
	int    n=1, iReNum=50;
	double fh, ft1, ft2, fx1=fa, fx2, fPx, fMPx, fSum=0;
	double fLogfa = log(fa), fLogfb = log(fb);

	fh  = (fLogfb-fLogfa)/iReNum;
	fPx = InsertValue_Linear( iMpNum, pfMP_M, pfMP_P, fMpwf );
	ft1 = CalNonDarcyFun_TwoPseudoPressure(fa, (float)fPx, fK, iGCPvtNum, pfP, pflg, pfUg, 
		                                   pfKroDivKrgPvt, iSatNum, pfSg, pfKrg, pfKroDivKrgSat);

	for( n=1; n<=iReNum; n++)
	{
		fx2  = exp(fLogfa + fh*n);
		fMPx = CalTMpbyRadius((float)fx2, fa, fb, fMpe, fMpwf);
		fPx  = InsertValue_Linear( iMpNum, pfMP_M, pfMP_P, (float)fMPx );
		ft2  = CalNonDarcyFun_TwoPseudoPressure((float)fx2, (float)fPx,  fK, iGCPvtNum, pfP, pflg, pfUg, 
		                                   pfKroDivKrgPvt, iSatNum, pfSg, pfKrg, pfKroDivKrgSat);
		fSum = fSum + (ft1+ft2)*(fx2-fx1)/2;