📄 flash.cpp
字号:
#include "stdafx.h"
#include "CVenus.h"
#include "resource.h"
#include "math.h"
#include <afxwin.h>
#include "fstream.h"
#include <strstrea.h>
#include <iostream.h>
#include "Flash.h"
#include "IPRCal.h"
#ifdef _DEBUG
#define new DEBUG_NEW
#undef THIS_FILE
static char THIS_FILE[] = __FILE__;
#endif
//=============================================================================
//=============================闪蒸计算构造函数================================
EOSFlashClass::EOSFlashClass( )
{
DS = 1;
m_iNum = 18;
m_pfTci = new float[m_iNum];
m_pfPci = new float[m_iNum];
m_pfWi = new float[m_iNum];
m_pfMi = new float[m_iNum];
m_pfVci = new float[m_iNum];
m_pfZci = new float[m_iNum];
m_pfKij = (float **)calloc( m_iNum, sizeof( float* ) );
for( int i=0; i<m_iNum; i++ )
m_pfKij[i] = (float*)calloc( m_iNum, sizeof( float ) );
m_pffiL = new float[m_iNum];
m_pffiV = new float[m_iNum];
m_pffi = new float[m_iNum+1];
m_pfKi = new float[m_iNum];
m_pfai = new float[m_iNum];
m_pfbi = new float[m_iNum];
m_pfmi = new float[m_iNum];
m_pfAlfai = new float[m_iNum];
m_pfZi = new float[m_iNum];
m_pfXi = new float[m_iNum+1];
m_pfYi = new float[m_iNum];
X1 = new float[m_iNum+1];
X2 = new float[m_iNum+1];
X3 = new float[m_iNum+1];
m_iPMP_Num = 111;
m_pfPMP_Press = new float[m_iPMP_Num];
m_pfPMP_MP = new float[m_iPMP_Num];
m_pfV = new float[m_iPMP_Num]; //气相摩尔分数
m_pfSg = new float[m_iPMP_Num]; //气相饱和度
m_pfDo = new float[m_iPMP_Num]; //油相密度
m_pfDg = new float[m_iPMP_Num]; //气相密度
m_pfZo = new float[m_iPMP_Num]; //油相偏差因子
m_pfZg = new float[m_iPMP_Num]; //气相偏差因子
m_pfUo = new float[m_iPMP_Num]; //油相粘度
m_pfUg = new float[m_iPMP_Num]; //气相粘度
m_pfMo = new float[m_iPMP_Num]; //油相分子量
m_pfMg = new float[m_iPMP_Num]; //气相分子量
}
//-----------------------------闪蒸计算析构函数--------------------------------
EOSFlashClass::~EOSFlashClass( )
{
delete []m_pfTci;
delete []m_pfPci;
delete []m_pfWi;
delete []m_pfMi;
delete []m_pfVci;
delete []m_pfZci;
for( int i=0; i<m_iNum; i++ )
free( m_pfKij[i] );
free( m_pfKij );
delete []m_pffiL;
delete []m_pffiV;
delete []m_pffi;
delete []m_pfKi;
delete []m_pfai; //状态方程系数
delete []m_pfbi; //状态方程系数
delete []m_pfmi;
delete []m_pfAlfai; //状态方程系数
delete []m_pfZi; //油相组成
delete []m_pfXi; //油相组成
delete []m_pfYi; //气相组成
delete []X1; //临时存储油相组成
delete []X2; //临时存储油相组成
delete []X3; //临时存储油相组成
delete []m_pfPMP_Press;
delete []m_pfPMP_MP;
delete []m_pfV; //气相摩尔分数
delete []m_pfSg; //气相饱和度
delete []m_pfDo; //油相密度
delete []m_pfDg; //气相密度
delete []m_pfZo; //油相偏差因子
delete []m_pfZg; //气相偏差因子
delete []m_pfUo; //油相粘度
delete []m_pfUg; //气相粘度
delete []m_pfMo; //油相分子量
delete []m_pfMg; //气相分子量
}
//------------------------------------------------------------------------------
//文件GCPVT.PVT中保存所有预定组分(18个组分)的所有有用PVT数据
//该函数根据输入组分m_pfZi从文件GCPVT.PVT中读取相应组分的PVT参数
int EOSFlashClass::InputCompFixData( float *m_pfZi )
{
int i, j, k, m;
float ffff;
char pBuf[256];
//按组分输入PVT参数
ifstream inFile( "GCPVT.PVT" );
inFile.getline( pBuf, 128 );
k = 0;
for( i=0; i<TOTALCOMNUM-1; i++ )
{
inFile.getline( pBuf, 128 );
inFile.getline( pBuf, 128 );
if( fabs(m_pfZi[i]) < 1e-8 )
continue;
istrstream pvtData( pBuf );
pvtData>>m_pfMi[k] >>m_pfTci[k] >>m_pfPci[k]
>>m_pfWi[k] >>m_pfVci[k] >>m_pfZci[k];
k++;
}
//按组分输入交互作用系数
inFile.getline( pBuf, 256 );
inFile.getline( pBuf, 256 );
k = 0;
for( i=0; i<TOTALCOMNUM-1; i++ )
{ //k
inFile.getline( pBuf, 256 );
inFile.getline( pBuf, 256 );
if( fabs(m_pfZi[i]) < 1e-8 )
continue;
istrstream pvtData( pBuf );
m = 0;
for( j=0; j<TOTALCOMNUM-1; j++ )
{ //m
if( fabs(m_pfZi[j]) < 1e-8 )
pvtData>>ffff;
else
pvtData>>m_pfKij[k][m++];
}
k++;
}
inFile.close( );
m_iNum = k;
//输出PVT参数
ofstream outFile( "pppp.dat" );
for( k=0; k<m_iNum; k++ )
{
outFile<<m_pfMi[k]<<" "<<m_pfTci[k]<<" "<<m_pfPci[k]
<<" "<<m_pfWi[k]<<" "<<m_pfVci[k]<<" "<< m_pfZci[k] << endl;
}
for( k=0; k<m_iNum; k++ )
{
for( j=0; j<m_iNum; j++ )
outFile << m_pfKij[k][j] << " ";
outFile << endl;
}
outFile.close( );
return 0;
}
//------------------------------------------------------------------------------
//Willsion 公式求平衡常数-----fP != 0 fT != 0
int EOSFlashClass::InitKi( float fT, float fP )
{
for( int i=0; i<m_iNum; i++ )
m_pfKi[i] = float( exp(5.37*(1+m_pfWi[i])*(1-m_pfTci[i]/fT))/(fP/m_pfPci[i]) );
return 0; //成功
}
//------------------------------------------------------------------------------
//计算PR_EOS-----mi-ai-bi
int EOSFlashClass::CalEOSCoeff_mi_ai_bi( )
{
for( int i=0; i<m_iNum; i++ )
{
m_pfmi[i] = float( 0.37464 + 1.54226*m_pfWi[i] - 0.26992*pow(m_pfWi[i],2) );
m_pfai[i] = float( 0.45724*pow(R*m_pfTci[i],2) / m_pfPci[i] );
m_pfbi[i] = float( 0.0778*R*m_pfTci[i] / m_pfPci[i] );
}
return 0;
}
//------------------------------------------------------------------------------
//计算状态方程系数 PR_EOS-------_InputData->fTemp != 0
int EOSFlashClass::CalEOSCoeff_Alfa_am_bm_Am_Bm( float fP, float fT, float *pfEi )
{
int i, j;
double fTmp;
for( i=0; i<m_iNum; i++ )
{
fTmp = 1.0 + m_pfmi[i] * (1.0-sqrt(fT/m_pfTci[i]));
m_pfAlfai[i] = float (fTmp*fTmp );
}
double famb = 0;
double fbmb = 0;
for( i=0; i<m_iNum; i++ )
{
for( j=0; j<m_iNum; j++ )
{
fTmp = pfEi[i]*pfEi[j]*sqrt(m_pfai[i]*m_pfai[j]*m_pfAlfai[i]*m_pfAlfai[j]);
famb = famb + fTmp*(1.0-m_pfKij[i][j]);
}
fbmb = fbmb + pfEi[i]*m_pfbi[i];
}
m_fam = (float)famb;
m_fbm = (float)fbmb;
m_fAm = float(m_fam*fP/pow(R*fT,2));
m_fBm = float(m_fbm*fP/(R*fT));
return 0; //成功
}
//------------------------------------------------------------------------------
//计算逸度 PR_EOS
int EOSFlashClass::Calfi( float fZ, float fP, float *pfEi, float *pffi )
{
int i, j;
double fAA;
double ftemp;
for( i=0; i<m_iNum; i++ )
{
fAA = 0;
for( j=0; j<m_iNum; j++ )
fAA = fAA + pfEi[j]*sqrt(m_pfai[i]*m_pfai[j]*m_pfAlfai[i]*m_pfAlfai[j])*(1.0-m_pfKij[i][j]);
if( fabs(fZ-m_fBm)<1e-18 )
{
char *pMsg = new char[80];
sprintf( pMsg, "fZ=%f fBm=%f", fZ, m_fBm );
MessageBox( 0, pMsg, "", MB_OK );
delete pMsg;
}
ftemp = (m_pfbi[i]/m_fbm)*(fZ-1.0) - log(fabs(fZ-m_fBm));
ftemp = ftemp - m_fAm/(2.828*m_fBm)*(2.0*fAA/m_fam-m_pfbi[i]/m_fbm)
*log(fabs((fZ+2.414*m_fBm)/(fZ-0.414*m_fBm)));
pffi[i] = float(exp(ftemp)*fP*pfEi[i]);
}
return 0; //成功
}
//------------------------------------------------------------------------------
//计算系数am与温度T的偏导数
int EOSFlashClass::CaldamdT( float fT, float *pfEi )
{
int i, j;
double fItem;
m_fdamdT = 0;
for( i=0; i<m_iNum; i++ )
{
for( j=0; j<m_iNum; j++ )
{
fItem = pfEi[i]*pfEi[j]*sqrt(m_pfai[i]*m_pfai[j])*(1.0-m_pfKij[i][j]);
fItem = fItem * ( m_pfmi[i]*sqrt(m_pfAlfai[j])/sqrt(m_pfTci[i]*fT)
+ m_pfmi[j]*sqrt(m_pfAlfai[i])/sqrt(m_pfTci[j]*fT) );
m_fdamdT = float( m_fdamdT+fItem );
}
}
m_fdamdT = -m_fdamdT/2;
return 0;
}
//------------------------------------------------------------------------------
//计算系数am与温度T的二阶偏导数
int EOSFlashClass::Cald2amdT( float fT, float *pfEi )
{
int i, j;
double fItem;
m_f2damdT = 0;
for( i=0; i<m_iNum; i++ )
{
for( j=0; j<m_iNum; j++ )
{
fItem = pfEi[i]*pfEi[j]*sqrt(m_pfai[i]*m_pfai[j])*(1.0-m_pfKij[i][j] );
fItem = fItem * ( m_pfmi[i]*(1+m_pfmi[j])/(4*fT*sqrt(m_pfTci[i]*fT))
+ m_pfmi[j]*(1+m_pfmi[i])/(4*fT*sqrt(m_pfTci[j]*fT)) );
m_f2damdT = float( m_f2damdT+fItem);
}
}
return 0;
}
//------------------------------------------------------------------------------
//计算偏差系数Z与温度T的偏导数
int EOSFlashClass::CaldZdT( float fP, float fT, float fZ, float *pfEi )
{
double fdAmdT, fdBmdT;
double fTmp1, fTmp2, fTmp3;
CaldamdT ( fT, pfEi );
Cald2amdT( fT, pfEi );
fdAmdT = fP*(m_fdamdT-2*m_fam/fT)/pow(R*fT,2);
fdBmdT = -m_fbm*fP/(R*fT*fT);
fTmp1 = fZ*fZ*fdBmdT + fZ*(fdAmdT-2*fdBmdT-6*m_fBm*fdBmdT);
fTmp2 = m_fBm*fdAmdT + (m_fAm-2*m_fBm-3*pow(m_fBm,2))*fdBmdT;
fTmp3 = 2*fZ*(1-m_fBm) - 3*fZ*fZ + 2*m_fBm+3*pow(m_fBm,2)-m_fAm;
m_fdZdT = float( (fTmp1-fTmp2)/fTmp3 );
return 0;
}
//------------------------------------------------------------------------------
//计算偏差系数Z与压力P的偏导数
int EOSFlashClass::CaldZdP( float fT, float fZ )
{
double fdAmdP, fdBmdP, fdZdP_b;
fdAmdP = m_fam/pow(R*fT,2);
fdBmdP = m_fbm/(R*fT);
fdZdP_b = (m_fAm-2*m_fBm-3*m_fBm*m_fBm)*fdBmdP + m_fBm*fdAmdP;
fdZdP_b = fdZdP_b - (fdAmdP-(2+6*m_fBm)*fdBmdP)*fZ - fdBmdP*fZ*fZ;
fdZdP_b = fdZdP_b /(3*fZ*fZ - 2*(1-m_fBm)*fZ + m_fAm-2*m_fBm-3*m_fBm*m_fBm);
m_fdZdP = float( fdZdP_b );
return 0;
}
//------------------------------------------------------------------------------
//设置气体混合物的组成
void EOSFlashClass::SetComposition( float *pfZi )
{
int k = 0;
for( int i=0; i<TOTALCOMNUM-1; i++ )
if( fabs(pfZi[i]) > 1e-8 )
m_pfZi[k++] = pfZi[i];
}
//------------------------------------------------------------------------------
//由物质平衡公式求气相分数 PR_EOS
int EOSFlashClass::InitV( float *pfKi )
{
int i, iMaxNum = 10000;
int KK = 0;
double fF1, fF2, fTmp;
double fV_b = m_fV;
double fEps = 0.0001;
while( KK < iMaxNum )
{
fF1 = 0; fF2 = 0;
for( i=0; i<m_iNum; i++ )
{
fTmp = 1.0+(pfKi[i]-1.0)*fV_b;
fF1 = fF1 + m_pfZi[i]*(pfKi[i]-1.0) / fTmp;
fF2 = fF2 + m_pfZi[i]*pow((pfKi[i]-1.0),2) / (fTmp*fTmp);
}
fV_b = fV_b + fF1/fF2;
if( fV_b > 1.0 )
{
fV_b = fV_b - fF1/fF2;
fV_b = (1.0-fV_b)/2.0 + fV_b;
}
else if( fV_b < 0.0 )
{
fV_b = fV_b - fF1/fF2;
fV_b = fV_b * 0.5;
}
else
{
if( fabs( fF1 ) < fEps )
{
m_fV = (float)fV_b;
return 0; //成功
}
}
KK++;
}
m_fV = 1;
return 1; //不成功
}
//------------------------------------------------------------------------------
//计算气液相的组成
int EOSFlashClass::InitXi_Yi( float *pfKi )
{
int i;
double fX = 0;
double fY = 0;
double fTmp;
for( i=0; i<m_iNum; i++ )
{
fTmp = 1.0+(pfKi[i]-1.0)*m_fV;
m_pfXi[i] = float(m_pfZi[i]/fTmp);
m_pfYi[i] = float(m_pfZi[i]*pfKi[i]/fTmp);
fX = fX + m_pfXi[i];
fY = fY + m_pfYi[i];
}
for( i=0; i<m_iNum; i++ )
{
m_pfXi[i] = float(m_pfXi[i]/fX);
m_pfYi[i] = float(m_pfYi[i]/fY);
}
return 0; //成功
}
//------------------------------------------------------------------------------
//计算Z因子 PR_EOS
int EOSFlashClass::CalZ_Factor( float fAm, float fBm, float &fInitZ )
{
int i=0, iMax = 10000;
double fEps = 0.001;
double fF1, fF2;
double fZ, fZ1;
fZ = fInitZ;
while( i < iMax )
{
fF1 = pow(fZ,3) - (1.0-fBm)*pow(fZ,2) + (fAm-2.0*fBm-3.0*pow(fBm,2))*fZ
- (fAm*fBm-pow(fBm,2)-pow(fBm,3));
fF2 = 3.0*pow(fZ,2) - 2.0*(1.0-fBm )*fZ + (fAm-2.0*fBm-3.0*pow(fBm,2));
fZ1 = fabs(fZ - fF1/fF2);
if( fabs(fZ1-fZ) < fEps )
{
fInitZ = (float)fZ1;
return 0; //成功
}
fZ = (float)fZ1;
i++;
}
fInitZ = 1;//fZ1;
return 1; //不成功
}
//------------------------------------------------------------------------------
//计算分子量&密度------kg/mol
void EOSFlashClass::CalMolecularWeight( )
{
m_fMg = 0;
m_fMo = 0;
for( int i=0; i<m_iNum; i++ )
{
m_fMg = m_fMg + m_pfMi[i]*m_pfYi[i];
m_fMo = m_fMo + m_pfMi[i]*m_pfXi[i];
}
m_fMg = m_fMg/1000;
m_fMo = m_fMo/1000;
}
//------------------------------------------------------------------------------
//计算质量密度------kg/m^3
void EOSFlashClass::CalQualityDensity( )
{
m_fDo = float(m_fPress*m_fMo / (m_fZo*R*m_fTemp) * pow(10,6));
m_fDg = float(m_fPress*m_fMg / (m_fZg*R*m_fTemp) * pow(10,6));
}
//------------------------------------------------------------------------------
//计算Mole密度------mol/m^3
void EOSFlashClass::CalMoleDensity( )
{
m_fDo = float(m_fPress / (m_fZo*R*m_fTemp));
m_fDg = float(m_fPress / (m_fZg*R*m_fTemp));
}
//------------------------------------------------------------------------------
//计算体积分数
void EOSFlashClass::CalVolumeFraction( )
{
m_fSg = m_fV*m_fZg/((m_fV*m_fZg)+((1-m_fV)*m_fZo));
if( m_fV < 0 )
MessageBox( 0, "fSg<0", "", MB_OK );
}
//------------------------------------------------------------------------------
void EOSFlashClass::CalGasSatuation( )
{
double fMoleDo, fMoleDg;
fMoleDo = m_fPress / (m_fZo*R*m_fTemp);
fMoleDg = m_fPress / (m_fZg*R*m_fTemp);
m_fSg = float(m_fV*fMoleDo/((m_fV*fMoleDo)+((1-m_fV)*fMoleDg)));
}
//------------------------------------------------------------------------------
//计算粘度
int EOSFlashClass::CalViscosity( float *pfEi, float fDen, float &fViscosity )
{
int i;
double fTpc=0, fPpc=0, fZpc=0, fVpc=0, fMpc=0;
double fTr, fUi, fKsai, fDr, fYita, fTmp;
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -