📄 flash.cpp
字号:
for( i=0; i<m_iNum; i++ )
{
fTpc = fTpc + pfEi[i]*m_pfTci[i];
fZpc = fZpc + pfEi[i]*m_pfZci[i];
fVpc = fVpc + pfEi[i]*m_pfVci[i];
fMpc = fMpc + pfEi[i]*m_pfMi[i];
}
fPpc = R*fZpc*fTpc/fVpc;
fKsai = pow(fTpc, 1.0/6.0)/(sqrt(fMpc)*pow(fPpc,2.0/3.0));
fTr = m_fTemp/fTpc;
if( fTr < 1.5 )
fUi = 3.4*pow(fTr,8.0/9.0)/fKsai;
else
fUi = 16.68*pow((0.1338*fTr-0.0932),5.0/9.0)/fKsai;
fYita = pow(fTpc,1.0/6.0)/(sqrt(fMpc)*pow(fPpc,2.0/3.0));
fDr = fDen*fZpc*R*fTpc/fPpc;
fTmp = 1.08*(exp(1.439*fDr)-exp(-1.111*pow(fDr,1.858)));
fViscosity = float((fUi + fTmp/fYita)/10000.0);
return 0;
}
/*int EOSFlashClass::CalViscosity( float *pfEi, float fDen, float &fViscosity )
{
int i;
float fTmp1=0, fTmp2=0, fTmp3=0, fTmp4=0, fTmp5=0, fTmp6=0;
float fTri, fTmp, fUi, fUXin, fKsai, fDr;
for( i=0; i<m_iNum; i++ )
{
fTmp1 = fTmp1 + pfEi[i]*m_pfTci[i];
fTmp2 = fTmp2 + pfEi[i]*m_pfPci[i];
fTmp3 = fTmp3 + pfEi[i]*m_pfMi[i];
fTri = fTemp/m_pfTci[i];
fTmp = pow(m_pfTci[i],1.0/6.0)
/( sqrt(m_pfMi[i])
*pow(m_pfPci[i],2.0/3.0) );
if( fTri < 1.5 )
fUi = 34*pow(10,-5)*pow(fTri,0.94)/fTmp;
else
fUi = 17.78*pow(10,-5)*pow((4.58*fTri-1.67),5.0/8.0)/fTmp;
fTmp4 = fTmp4 + pfEi[i]*sqrt(m_pfMi[i])*fUi;
fTmp5 = fTmp5 + pfEi[i]*sqrt(m_pfMi[i]);
fTmp6 = fTmp6 + pfEi[i]*m_pfVci[i];
}
fUXin = fTmp4/fTmp5;
fKsai = pow(fTmp1,1.0/6.0)/(sqrt(fTmp3)*pow(fTmp2,2.0/3.0));
fDr = fDen/fTmp6;
fDr = fDen*fTmp6;
fTmp = 0.1023+0.023364*fDr+0.058533*fDr*fDr-0.40758*pow(fDr,3)+0.0093324*pow(fDr,4);
fViscosity = fUXin + (pow(fTmp,4)-pow(10,-4))/fKsai;
return 0;
}
*/
//------------------------------------------------------------------------------
// PR_EOS
int EOSFlashClass::CalZV( )
{
int i, LL = 0;
int iIterNum,iMaxIterNum=100;
float fF1;
double fF2;
double fT0 = 0, fAB0 = 0, fAB;
float fBackV = m_fV;
iIterNum = 0;
InitKi( m_fTemp, m_fPress );
CalEOSCoeff_mi_ai_bi( );
while ( iIterNum < iMaxIterNum )
{
m_fV = fBackV;
if( InitV( m_pfKi ) ) //计算气相mol分数
return 1; //不成功
if( InitXi_Yi( m_pfKi ) ) //计算油,气相的组成
return 2; //不成功
while( LL < 1000 )
{
CalEOSCoeff_Alfa_am_bm_Am_Bm( m_fPress, m_fTemp, m_pfYi );
fF1 = 1.0;
if( CalZ_Factor( m_fAm, m_fBm, fF1 ) )
return 3;
m_fZg = fF1;
//outFile << "Zg=" << fZg << endl;
if( Calfi( m_fZg, m_fPress, m_pfYi, m_pffiV ) )
return 4;
//for( i=0; i<iNum; i++ )
// outFile << "Fivi=" << m_pffiV[i] << endl;
CalEOSCoeff_Alfa_am_bm_Am_Bm( m_fPress, m_fTemp, m_pfXi );
fF1 = 0;
if( CalZ_Factor( m_fAm, m_fBm, fF1 ) )
return 5;
m_fZo = fF1;
//outFile << "Zg=" << fZo << endl;
if( Calfi( m_fZo, m_fPress, m_pfXi, m_pffiL ) )
return 6;
//for( i=0; i<iNum; i++ )
// outFile << "Fili=" << m_pffiL[i] << endl;
fF1 = 0;
for( i=0; i<m_iNum; i++ )
fF1 = float(fF1+pow(m_pffiL[i]/m_pffiV[i]-1.0,2));
if( fF1 < 0.000001 )
{ //4610
fAB = 0;
for( i=0; i<m_iNum; i++ )
fAB = fAB + m_pfZi[i] * ( m_pfKi[i] - 1.0 )
/ ( 1.0 + ( m_pfKi[i] - 1.0 ) * m_fV );
if( fabs( fAB ) < 0.0001 )
{ //4700
return 0; //成功
}
if( LL > 1 )
{ //4660
fF2 = -( m_fTemp - fT0 ) * fAB / ( fAB - fAB0 );
fT0 = m_fTemp;
fAB0 = fAB;
m_fTemp = float(m_fTemp+fF2);
//goto 4650
}
else
{
fT0 = m_fTemp;
fAB0 = fAB;
m_fTemp = float(0.98*m_fTemp);
//goto 4650
}
//4650
LL++; //goto 4430
}
else
{ //4590
for( i=0; i<m_iNum; i++ )
m_pfKi[i] = m_pfKi[i] * m_pffiL[i] / m_pffiV[i];
iIterNum++;
//goto 4230
break;
}
}
}
return 10;
}
//------------------------------------------------------------------------------
//计算露点压力 PR_EOS
int EOSFlashClass::CalDewPressure( )
{
int i, iResFlag;
int LOOP1, Line900;
int iLoopNum1, iLoopNum2, iLoopNum3;
int LC = 0;
int LK, LT;
float AB;
double P1, P0;
double SS, XX;
double SG, SGYZ = 1.05;
iLoopNum1 = 0;
P1 = m_fPress;
CalEOSCoeff_mi_ai_bi( );
do
{ //500
iLoopNum1++; LK = 1; LT = 0;
iResFlag = InitKi( m_fTemp, m_fPress );
if( iResFlag ) return 1;
for( i=0; i<m_iNum; i++ )
m_pfYi[i] = m_pfZi[i];
iLoopNum2 = 0;
do
{
//530
iLoopNum2++;
CalEOSCoeff_Alfa_am_bm_Am_Bm( m_fPress, m_fTemp, m_pfYi );
AB = 1.0;
if( CalZ_Factor( m_fAm, m_fBm, AB ) )
return 3;
m_fZg = AB;
if( m_fZg-m_fBm <= 0 )
{
m_fPress = m_fPress - m_fPStep;
P1 = m_fPress;
//if( P1 <= 0.0001 ) return 1;
//910输入压力值偏高,现已降到P,并重新 拟合理
//做气层循环 goto 500
break;
}
//600
if( Calfi( m_fZg, m_fPress, m_pfYi, m_pffiV ) )
return 4;
//630
iLoopNum3 = 0;
do
{
iLoopNum3++;
if( LT == 0 )
{
for( i=0; i<m_iNum; i++ )
m_pfXi[i] = m_pfZi[i] / m_pfKi[i];
if( LK < 1.5 )
{ XX = 0;
for( i=0; i<m_iNum; i++ )
XX = XX + m_pfXi[i];
for( i=0; i<m_iNum; i++ )
m_pfXi[i] = float(m_pfXi[i]/XX);
}
}
//660 Program
CalEOSCoeff_Alfa_am_bm_Am_Bm( m_fPress, m_fTemp, m_pfXi );
AB = 0;
if( CalZ_Factor( m_fAm, m_fBm, AB ) )
{
m_fPress = m_fPress - m_fPStep;
P1 = m_fPress;
//if( P1 <= 0.0001 ) return 1;
//910
LOOP1 = 1;
break; //goto 500 做外层气相循环
}
m_fZo = AB;
if( m_fZo-m_fBm <= 0 )
{
m_fPress = m_fPress - m_fPStep;
P1 = m_fPress;
//if( P1 <= 0.0001 ) return 1;
//910 输入压力值偏高,现已降到P,并重新 拟合理"<<P1<<endl;
LOOP1 = 1;
break;
}
//740
if( Calfi( m_fZo, m_fPress, m_pfXi, m_pffiL ) )
return 6;
SS = 0;
for( i=0; i<m_iNum; i++ )
SS = SS + pow( m_pffiL[i]/m_pffiV[i] - 1.0, 2 );
if( SS > 1e-8 )
{
for( i=0; i<m_iNum; i++ )
m_pfKi[i] = m_pfKi[i]*m_pffiL[i]/m_pffiV[i];
LT = 0; LC++; LK = 2;
//goto 630
continue;
}
//800
AB = 0; SG = 0;
for( i=0; i<m_iNum; i++ )
AB = AB + m_pfZi[i] * ( m_pfKi[i] - 1 ) / m_pfKi[i];
for( i=0; i<m_iNum; i++ )
SG = SG + m_pfZi[i] / m_pfKi[i];
if( fabs(AB) < 0.001 )
Line900 = 1;
else
{
Line900 = 0;
P0 = m_fPress;
if( SG > SGYZ )
SG = SGYZ;
m_fPress = float(P0*SG);
if( fabs( m_fPress - P0 ) > 0.01 )
{
LT = 1;
//goto 530
LOOP1 = 0;
break;
}
else
Line900 = 1;
}
if( Line900 == 1 )
{ //900
if( m_fPress != P1 )
{ //930
m_fDewPressure = m_fPress;
//MessageBox( 0, "Dew Pressure", "", MB_OK );
return 0;
}
else
{
m_fPress = m_fPress - m_fPStep;
P1 = m_fPress;
//if( P1 <= 0.0001 ) return 1;
//goto 500
LOOP1 = 1; break;
}
}
}while( iLoopNum3 < 1000 );//end 630
if( LOOP1 == 1 ) break;
}while( iLoopNum2 < 100 );//end 530
}while( iLoopNum1 < 100 ); //end 500
if( iLoopNum1 == 100 )
{
//MessageBox( 0, "Dew Pressure", "", MB_OK );
return 1;
}
else
{
m_fDewPressure = m_fPress;
return 0;
}
}
//------------------------------------------------------------------------------
//构造目标函数
int EOSFlashClass::ObjectEquation( )
{
int i;
double fAB = 0;
float fTmp;
//计算气相
CalEOSCoeff_Alfa_am_bm_Am_Bm( m_fPress, m_fTemp, m_pfYi );
fTmp = 1.0;
if( CalZ_Factor( m_fAm, m_fBm, fTmp ) )
return 1;
m_fZg = fTmp;
if( Calfi( m_fZg, m_fPress, m_pfYi, m_pffiV ) )
return 2;
//计算油相
CalEOSCoeff_Alfa_am_bm_Am_Bm( m_fPress, m_fTemp, m_pfXi );
fTmp = 0;
if( CalZ_Factor( m_fAm, m_fBm, fTmp ) )
return 3;
m_fZo = fTmp;
if( Calfi( m_fZo, m_fPress, m_pfXi, m_pffiL ) )
return 4;
for( i=0; i<m_iNum; i++ )
{
fTmp = (m_pffiL[i]/m_pfXi[i]) / (m_pffiV[i]/m_pfYi[i]);
fAB = fAB + m_pfZi[i] * ( fTmp - 1.0 ) / ( 1.0 + ( fTmp - 1.0 ) * m_fV );
m_pffi[i] = m_pffiL[i] - m_pffiV[i];
}
m_pffi[m_iNum] = (float)fAB;
m_pfXi[m_iNum] = m_fV;
return 0; //成功
}
//------------------------------------------------------------------------------
//单步叠代
int EOSFlashClass::Step( float *PP, float **H )
{
int i, j, iResFlag;
float SA = 0, SB = 0;
float *F;
float *U, *Y1;
F = m_pffi;
U = new float[m_iNum+1];
Y1 = new float[m_iNum+1];
for( i=0; i<m_iNum+1; i++ )
{
m_pfXi[i] = m_pfXi[i] + PP[i];
U[i] = F[i];
}
m_fV = m_pfXi[m_iNum];
for( i=0; i<m_iNum; i++ )
m_pfYi[i] = ( m_pfZi[i] - m_pfXi[i]*(1-m_fV) ) / m_fV;
iResFlag = ObjectEquation( );
if( iResFlag )
return iResFlag; //不成功
for( i=0; i<m_iNum+1; i++ )
Y1[i] = F[i] - U[i];
SA = 0;
for( i=0; i<m_iNum+1; i++ )
{
SB = 0;
for( j=0; j<m_iNum+1; j++ )
SB = SB + H[i][j] * Y1[j];
U[i] = SB - PP[i];
SA = SA + SB * PP[i];
}
if( SA == 0 )
{
iResFlag = 2;
return iResFlag;
}
for( j=0; j<m_iNum+1; j++ )
{
SB = 0;
for( i=0; i<m_iNum+1; i++ )
SB = SB + PP[i] * H[i][j];
SB = SB /SA;
for( i=0; i<m_iNum+1; i++ )
H[i][j] = H[i][j] - SB * U[i];
}
return 0;
}
//------------------------------------------------------------------------------
/*int EOSFlashClass::BroyDen( )
{
int IC = 0, iResFlag, iRes;
int MAXF=1000;
int i, j;
float SA = 0, EPS=float(0.00001);
float *F, fV;
float *PP;
float *U, *Y1;
float **H;
PP = new float[m_iNum+1];
U = new float[m_iNum+1];
Y1 = new float[m_iNum+1];
H = (float **)calloc( m_iNum+1, sizeof( float* ) );
for( i=0; i<m_iNum+1; i++ )
H[i] = (float*)calloc( m_iNum+1, sizeof( float ) );
F = m_pffi;
iResFlag = ObjectEquation( );
if( iResFlag )
{
delete PP;
delete U;
delete Y1;
for( i=0; i<m_iNum+1; i++ )
free( H[i] );
free( H );
return 5;
}
for( i=0; i<m_iNum+1; i++ )
{
PP[i] = 0;
for( j=0; j<m_iNum+1; j++ )
H[i][j] = 0;
}
for( i=0; i<m_iNum+1; i++ )
H[i][i] = 1;
for( i=0; i<m_iNum+1; i++ )
{
PP[i] = float(0.001);
iResFlag = Step( PP, H );
if( iResFlag )
{
delete PP;
delete U;
delete Y1;
for( i=0; i<m_iNum+1; i++ )
free( H[i] );
free( H );
return 2;
}
PP[i] = 0;
}
while( IC < MAXF )
{
for( i=0; i<m_iNum+1; i++ )
{
SA = 0;
for( j=0; j<m_iNum+1; j++ )
SA = SA - H[i][j]*F[j];
PP[i] = SA;
}
iResFlag = Step( PP, H );
if( iResFlag )
break;
SA = 0;
for( i=0; i<m_iNum+1; i++ )
SA = SA + F[i] * F[i];
if( SA < EPS )
{
fV = m_pfXi[m_iNum];
for( i=0; i<m_iNum; i++ )
m_pfYi[i] = ( m_pfZi[i] - (1-fV)*m_pfXi[i] ) / fV;
if( DS <= 1 )
{
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -