📄 iprcal.cpp
字号:
return float(0.0003458*fTemp*fZ/fPress);
}
//===============================计算气体粘度计算==============================
//fPress——压力(兆帕)
//fTemp——温度(开氏度)
//frg——气体比重
//fZ——气体偏差因子
//返回值——气体粘度(mPa.s)
float Ug_Lee_Gonzalez_Eakin( float fPress, float fTemp, float frg, float fZ )
{
double fMg = 28.96*frg;
double fItem1, fItem2, fItem3, fItem4;
fItem1=(9.4+0.02*fMg)*pow(1.8*fTemp,1.5)/(209+19*fMg+1.8*fTemp);//UK
fItem2=3.5+986/(1.8*fTemp)+0.01*fMg; //UX
fItem3=2.4-0.2*fItem2; //UY
//fItem4=16.0185*0.001*2.699*frg*142.23*fPress/fZ/fTemp; //密度 Rg----天然气比重
fItem4=3.4844*frg*fPress/fZ/fTemp; //密度 Rg----天然气比重
return float(fItem1*exp(fItem2*pow(fItem4,fItem3))/10000);
}
//fPc——拟临界压力(兆帕)
//fTc——拟临界温度(开氏度)
//fN2——氮气摩尔分数
//fCO2——二氧化碳摩尔分数
//fH2S——硫化氢摩尔分数
float Ug_Standing(float fPress, float fTemp, float fPc, float fTc,
float frg, float fN2, float fCO2, float fH2S)
{
double fPr = fPress / fPc;
double fTr = fTemp / fTc;
double U1, dUN2, dUCO2, dUH2S;
double Ug;
U1 = ( 1.709*pow(10.0,-5) - 2.062*pow(10.0,-6)*frg)
* ( 1.8*(fTemp-273.15)+32 ) + 8.188*pow(10.0,-3)
- 6.15*pow(10.0,-3)*log10(frg);
dUN2 = fN2 * ( 8.48*log10(frg) + 9.59 ) * 0.001;
dUCO2 = fCO2 * ( 9.08*log10(frg) + 6.24 ) * 0.001;
dUH2S = fH2S * ( 8.49*log10(frg) + 3.73 ) * 0.001;
U1 = U1 + dUN2 + dUCO2 + dUH2S;
Ug = -2.4621182 + 2.97054714*fPr - 0.286264054*fPr*fPr + 0.00850420522*fPr*fPr*fPr
+ fTr*( 2.80860949 - 3.49803305*fPr + 0.36037302*fPr*fPr - 0.01044324*fPr*fPr*fPr )
+ fTr*fTr*( -0.793385684 + 1.39643306*fPr - 0.149144925*fPr*fPr + 0.00441015512*fPr*fPr*fPr )
+fTr*fTr*fTr*( 0.0839387178 - 0.186408848*fPr + 0.0203367881*fPr*fPr - 0.000609579263*fPr*fPr*fPr );
return float(U1*exp(Ug)/fTr);
}
float Ug_Calculate( float fTc, float fPc, float fTemp, float fPress, float fZ,
float frg, float fN2, float fCO2, float fH2S, int iUgMethod)
{
switch( iUgMethod )
{
case UG_LEE_GONZALEA_EAKIN:
return Ug_Lee_Gonzalez_Eakin(fPress,fTemp,frg,fZ);
case UG_STANDING:
return Ug_Standing(fPress,fTemp,fPc,fTc,frg,fN2,fCO2,fH2S);
}
return Ug_Standing(fPress,fTemp,fPc,fTc,frg,fN2,fCO2,fH2S);
}
//==========================设置拟压力计算的压力段数据=========================
//pfPress——保存设置的压力序列
//iNum——压力数据点数
void SetPMP_Pressure(float* pfPress, int iNum)
{
float fPressure=(float)-0.15, fStep;
for( int i=0; i<iNum; i++)
{
if( fPressure < 0.8 )
fStep = (float)0.2;
else if( fPressure < 14.5 )
fStep=(float)0.5;
else if( fPressure < 35 )
fStep = (float)1.0;
else fStep = (float)2.0;
fPressure += fStep;
pfPress[i] = fPressure;
}
}
//==========================单相气体拟压力计算=================================
//iNum——压力数据点数
//pfPress——压力序列(MPa)
//pfMp——计算出的拟压力序列(MPa^2/mPa.s)
//fPseudoPc——拟临界压力(兆帕)
//fPseudoTc——拟临界温度(开氏度)
//fResTemp——地层温度(开氏度)
//frg——气体比重
//fN2——氮气摩尔分数
//fCO2——二氧化碳摩尔分数
//fH2S——硫化氢摩尔分数
//iUgMethod——计算粘度的方法
//iZgMethod——计算Z因子的方法
void CalGasWellPseudoPressure( int iNum,
float *pfPress,
float *pfMp,
float fPseudoPc,
float fPseudoTc,
float fResTemp,
float frg,
float fN2,
float fCO2,
float fH2S,
int iUgMethod,
int iZgMethod)
{
int i;
double startFun,midFun,endFun,fSum;
float* pfZg=new float[iNum];
float* pfUg=new float[iNum];
for(i=0;i<iNum;i++)
{
pfZg[i] = Zg_Calculate( fPseudoTc,fPseudoPc,fResTemp,pfPress[i],iZgMethod );
pfUg[i] = Ug_Calculate( fPseudoTc,fPseudoPc,fResTemp,pfPress[i],
pfZg[i], frg, fN2, fCO2, fH2S, iUgMethod );
}
pfMp[0] = pfPress[0] * pfPress[0] / pfUg[0] / pfZg[0];
fSum = pfMp[0];
for( i=1; i<iNum; i++ )
{
midFun = 4.0 * ( pfPress[i-1] + pfPress[i] )
/ ( pfUg[i-1] + pfUg[i] ) / ( pfZg[i-1] + pfZg[i] );
startFun = 2.0 * pfPress[i-1] / pfUg[i-1] / pfZg[i-1];
endFun = 2.0 * pfPress[i] / pfUg[i] / pfZg[i];
fSum = fSum + ( startFun + 2 * midFun + endFun )
* ( pfPress[i] - pfPress[i-1] ) / 4;
pfMp[i] = (float)fSum;
}
delete []pfZg;
delete []pfUg;
}
//============================单点线性插值=====================================
float InsertValue_Linear( int num, //理论点数
float *TheoryP, //插值对象
float *TheoryMp,
float RealP ) //被插之值
{
int i;
float ans;
for( i=0; i<num-1; i++ )
{
if( TheoryP[i] <= RealP && TheoryP[i+1] > RealP )
break;
if( TheoryP[i] >= RealP && TheoryP[i+1] < RealP )
break;
}
if( i==num-1 )
{
if(TheoryP[0]<TheoryP[num-1])
{
if(RealP<=TheoryP[0]) i=0;
if(RealP>=TheoryP[num-1]) i=num-2;
}
else
{
if(RealP<=TheoryP[num-1]) i=num-2;
if(RealP>=TheoryP[0]) i=0;
}
}
ans=(RealP-TheoryP[i])*(TheoryMp[i+1]-TheoryMp[i])/(TheoryP[i+1]-TheoryP[i])+TheoryMp[i];
return ans;
}
//============================多点线性插值=====================================
void MultiInsertValue_Linear( int iNum, //理论点数
float *TheoryP, //插值对象
float *TheoryMp,
int DestNum, //被插值点数
float *pRealP, //被插的值
float *pRealMp ) //输出结果
{
for( int i=0; i<DestNum; i++ )
pRealMp[i] = InsertValue_Linear( iNum, TheoryP, TheoryMp, pRealP[i] );
}
//=============================================================================
//根据凝析油气PVT参数计算油相相对渗透率与气相相对渗透率之比
//pfKroDivKrg——保存Kro/Krg
//iGCPvtNum——凝析油气PVT数据点数
//pfV——气相摩尔分数,
//pfDo——油相摩尔密度
//pfDg——气相摩尔密度
//pfUo——油相粘度(mPa.s)
//pfUg——气相粘度(mPa.s)
int KroDivKrgByPvt( float* pfKroDivKrg, int iGCPvtNum, float* pfV,
float* pfDo, float* pfDg, float* pfUo, float* pfUg )
{
double fTmp1, fTmp2;
for( int i=0; i<iGCPvtNum; i++ )
{
if( fabs(pfV[i]-1) > 1e-15 )
{
fTmp1 = (1-pfV[i])*pfDg[i]*pfUo[i];
fTmp2 = pfV[i]*pfDo[i]*pfUg[i];
pfKroDivKrg[i] = float(fTmp1/fTmp2);
}
else
pfKroDivKrg[i] = 0;
}
return 0;
}
//=============================================================================
//根据凝析油气相对渗透率曲线计算油相相对渗透率与气相相对渗透率之比
//pfKroDivKrg——保存Kro/Krg
//iSatNum——相对渗透率数据点数
//pfKro——油相相对渗透率
//pfKrg——气相相对渗透率
int KroDivKrgBySat( float* pfKroDivKrg, int iSatNum, float* pfKro, float* pfKrg )
{
for( int i=0; i<iSatNum; i++ )
pfKroDivKrg[i] = pfKro[i]/pfKrg[i];
return 0;
}
//=============================================================================
//插值计算某压力对应的油相相对渗透率与气相相对渗透率之比Kro/Krg
//iSatNum——相对渗透率数据点数
//pfSg、pfKroDivKrgSat——气相饱和度Sg、对应的油相相对渗透率与气相相对渗透率之比Kro/Krg
//iGCPvtNum——凝析油气PVT数据点数
//pfP(兆帕)、pfKroDivKrgPVT——压力、对应的油相相对渗透率与气相相对渗透率之比Kro/Krg
//fPress(兆帕)——所求Kro/Krg对应的压力
//返回值为待求的气相饱和度
float GetSgByPress( int iSatNum, float* pfKroDivKrgSat, float* pfSg,
int iGCPvtNum, float* pfKroDivKrgPVT, float* pfP, float fPress )
{
float fSg;
float fKroDivKrg;
fKroDivKrg = InsertValue_Linear( iGCPvtNum, pfP, pfKroDivKrgPVT, fPress );
fSg = InsertValue_Linear( iSatNum, pfKroDivKrgSat, pfSg, fKroDivKrg );
return fSg;
}
//=============================================================================
//根据相对渗透率数据求取某一饱和度值对应的油相相对渗透率
//iSatNum——相对渗透率数据点数
//pfSg、pfKro——气相饱和度、对应的油相相对渗透率
//fSg——所求油相相对渗透率对应的气相饱和度
//返回值为待求的油相相对渗透率
float GetKro( int iSatNum, float* pfKro, float* pfSg, float fSg )
{
return InsertValue_Linear( iSatNum, pfSg, pfKro, fSg );
}
//=============================================================================
//根据相对渗透率数据求取某一饱和度值对应的气相相对渗透率
//iSatNum——相对渗透率数据点数
//pfSg、pfKro——气相饱和度、对应的气相相对渗透率
//fSg——所求气相相对渗透率对应的气相饱和度
//返回值为待求的气相相对渗透率
float GetKrg( int iSatNum, float* pfKrg, float* pfSg, float fSg )
{
return InsertValue_Linear( iSatNum, pfSg, pfKrg, fSg );
}
//=============================================================================
//计算两相拟压力
//iGCPvtNum——凝析油气PVT数据点数
//pfP(兆帕)——求解PVT参数的压力
//pfV(摩尔分数)——计算的气相摩尔分数
//pfDg(千克/立方米)——气相密度
//pfUg(毫帕.秒)——气相粘度
//pfDo(千克/立方米)——油相密度
//pfUo(毫帕.秒)——油相粘度
//iSatNum——相对渗透率数据点数
//pfSg——气相饱和度系列
//pfKro——油相相对渗透率
//pfKrg——气相相对渗透率
//iMpNum——计算的拟压力点数
//pfPress(兆帕)——求解拟压力对应的压力
//pfMp——求解的拟压力
int PseudoTwoPhase( int iGCPvtNum, float* pfP, float* pfV,
float* pfDg, float* pfUg, float* pfDo, float* pfUo,
int iSatNum, float* pfSg, float* pfKro, float* pfKrg,
int iMpNum, float* pfPress, float* pfMp)
{
float* pfKroDivKrgPvt = new float[iGCPvtNum];
float* pfKroDivKrgSat = new float[iSatNum];
KroDivKrgByPvt( pfKroDivKrgPvt, iGCPvtNum, pfV, pfDo, pfDg, pfUo, pfUg );
KroDivKrgBySat( pfKroDivKrgSat, iSatNum, pfKro, pfKrg );
float fPress, fV, fSg, fDg, fDo, fUg, fUo, fKro, fKrg;
double fStartFun, fMidFun, fEndFun, fSum;
fV = InsertValue_Linear( iGCPvtNum, pfP, pfV, pfPress[0] );
fDg = InsertValue_Linear( iGCPvtNum, pfP, pfDg, pfPress[0] );
fDo = InsertValue_Linear( iGCPvtNum, pfP, pfDo, pfPress[0] );
fUg = InsertValue_Linear( iGCPvtNum, pfP, pfUg, pfPress[0] );
fUo = InsertValue_Linear( iGCPvtNum, pfP, pfUo, pfPress[0] );
fSg = GetSgByPress( iSatNum, pfKroDivKrgSat, pfSg, iGCPvtNum, pfKroDivKrgPvt, pfP, pfPress[0] );
fKrg = GetKrg( iSatNum, pfKrg, pfSg, fSg );
fKro = GetKro( iSatNum, pfKro, pfSg, fSg );
if(fUo<1e-6)
fStartFun = fDg*fKrg/fUg;
else
fStartFun = fDo*fKro/fUo+fDg*fKrg/fUg;
pfMp[0] = float(0.5*fStartFun*pfPress[0]);
fSum = pfMp[0];
for( int i=1; i<iMpNum; i++ )
{
fPress = pfPress[i];
fV = InsertValue_Linear( iGCPvtNum, pfP, pfV, fPress );
fDg = InsertValue_Linear( iGCPvtNum, pfP, pfDg, fPress );
fDo = InsertValue_Linear( iGCPvtNum, pfP, pfDo, fPress );
fUg = InsertValue_Linear( iGCPvtNum, pfP, pfUg, fPress );
fUo = InsertValue_Linear( iGCPvtNum, pfP, pfUo, fPress );
fSg = GetSgByPress( iSatNum, pfKroDivKrgSat, pfSg, iGCPvtNum, pfKroDivKrgPvt, pfP, fPress );
fKrg = GetKrg( iSatNum, pfKrg, pfSg, fSg );
fKro = GetKro( iSatNum, pfKro, pfSg, fSg );
if(fUo<1e-6)
fEndFun = fDg*fKrg/fUg;
else
fEndFun = fDo*fKro/fUo+fDg*fKrg/fUg;
fPress = float(0.5*(pfPress[i-1]+pfPress[i]));
fV = InsertValue_Linear( iGCPvtNum, pfP, pfV, fPress );
fDg = InsertValue_Linear( iGCPvtNum, pfP, pfDg, fPress );
fDo = InsertValue_Linear( iGCPvtNum, pfP, pfDo, fPress );
fUg = InsertValue_Linear( iGCPvtNum, pfP, pfUg, fPress );
fUo = InsertValue_Linear( iGCPvtNum, pfP, pfUo, fPress );
fSg = GetSgByPress( iSatNum, pfKroDivKrgSat, pfSg, iGCPvtNum, pfKroDivKrgPvt, pfP, fPress );
fKrg = GetKrg( iSatNum, pfKrg, pfSg, fSg );
fKro = GetKro( iSatNum, pfKro, pfSg, fSg );
if(fUo<1e-6)
fMidFun = fDg*fKrg/fUg;
else
fMidFun = fDo*fKro/fUo+fDg*fKrg/fUg;
fSum = fSum + ( fStartFun+4*fMidFun+fEndFun ) * ( pfPress[i]-pfPress[i-1] ) / 6;
pfMp[i] = (float)fSum;
fStartFun = fEndFun;
}
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -