flash.cpp

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

					 for( i=0; i<m_iNum+1; i++ )
							X1[i] = m_pfXi[i];
					 fP1 = m_fPress;
				}
				else if( DS <= 2 )
				{
					 for( i=0; i<m_iNum+1; i++ )
							X2[i] = m_pfXi[i];
					 fP2 = m_fPress;
				}
				else if( DS < 4 )
				{
					 for( i=0; i<m_iNum+1; i++ )
							X3[i] = m_pfXi[i];
					 fP3 = m_fPress;
				}
				else
				{
					 for( i=0; i<m_iNum+1; i++ )
					 {
							X1[i] = X2[i];
							X2[i] = X3[i];
							X3[i] = m_pfXi[i];
					 }
					 fP1 = fP2;
					 fP2 = fP3;
					 fP3 = m_fPress;
				}
				iRes = 0;
				break;     //成功
		  }
	 }
	 if( IC == 1000 && iRes != 0 )
		  iResFlag = 10;

	 delete PP;
	 delete U;
	 delete Y1;
	 for( i=0; i<m_iNum+1; i++ )
			free( H[i] );
	 free( H );

	 return iResFlag;	    //不成功
}*/
//------------------------------------------------------------------------------
//处理下一步的闪蒸数据
int EOSFlashClass::InitNextFlashSpot(  )
{
	 int             i;
	 float           L;

	 DS++;
	 if( DS > 100 )
		  DS = 5;

	 m_fV = m_pfXi[m_iNum];
	 L = 1 - m_fV;
	 if( L < 0.4 )
		  m_fPress = m_fPress - m_fPStep;
	 else if( L < 0.6 )
		  m_fPress = m_fPress - 1;
	 else
		  m_fPress = float(m_fPress-0.1);

	 if( DS > 4 )
		  for( i=0; i<m_iNum+1; i++ )
		  {
				 m_pfXi[i] =  ( (m_fPress-fP2)*(m_fPress-fP3)*X1[i] ) / ( (fP1-fP2)*(fP1-fP3) );
				 m_pfXi[i] += ( (m_fPress-fP1)*(m_fPress-fP3)*X2[i] ) / ( (fP2-fP1)*(fP2-fP3) );
				 m_pfXi[i] += ( (m_fPress-fP1)*(m_fPress-fP2)*X3[i] ) / ( (fP3-fP1)*(fP3-fP2) );
		  }
	 m_fV = m_pfXi[m_iNum];
	 for( i=0; i<m_iNum; i++ )
			m_pfYi[i] = ( m_pfZi[i] - (1-m_fV)*m_pfXi[i] ) / m_fV;
	 return 0;
}

//------------------------------------------------------------------------------
int EOSFlashClass::BroyDen(  )
{
	 int      LB,     iResFlag, iIterNum, iMaxIterNum=100;
	 int      i,      j,        K,        iRes;
	 double   SA = 0, SB = 0,   EPS=0.0000001;
	 float    *F;

	 float    *PP, *U, *Y1, **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 ) );

	 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;

	 iResFlag =  ObjectEquation(  );
	 F        =  m_pffi;

	 K        = 0;
	 iIterNum = 0;
	 while( iIterNum < iMaxIterNum && iResFlag == 0 )
	 {   if( K < m_iNum+1 )
		  {   //3460
				PP[K] = float(0.001);
				LB    = -3;
				//goto 3860
		  }
		  else
		  {   //3520
				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] = float(SA);
				}
				LB = -2;
				//goto 38600
		  }

		  //3860  program
		  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 )
				break;  //不成功

		  iIterNum++;
		  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] = float(SB-PP[i]);
				 SA = SA + SB * PP[i];
		  }

		  if( fabs(SA) < 1e-17 )
		  {
				m_fV = m_pfXi[m_iNum];
				iRes = 10;
				continue;
				//break;    //?不成功
		  }

		  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] = float(H[i][j]-SB*U[i]);
		  }

		  if( LB == -3 )
		  {   //3490
				PP[K] = 0;
				K++;
				//goto 3450
		  }
		  else
		  {    //3610
				 SA = 0;
				 for( i=0; i<m_iNum+1; i++ )
						SA = SA + F[i] * F[i];
				 if( SA > EPS )
				 {
					  if( iIterNum > iMaxIterNum )
					  {   //3840
							iResFlag = 10;
							break;   //不成功
					  }
					  else
					  {
						  //3520
					  }
				 }
				 else
				 {
					  m_fV  = m_pfXi[m_iNum];
					  for( i=0; i<m_iNum; i++ )
							 m_pfYi[i] = ( m_pfZi[i] - (1-m_fV)*m_pfXi[i] ) / m_fV;
					  if( DS <= 1 )
					  {
							for( i=0; i<m_iNum+1; i++ )
								  X1[i] = m_pfXi[i];
							fP1 = m_fPress;
					  }
					  else if( DS <= 2 )
					  {
							for( i=0; i<m_iNum+1; i++ )
								  X2[i] = m_pfXi[i];
							fP2 = m_fPress;
					  }
					 else if( DS < 4 )
					 {
							for( i=0; i<m_iNum+1; i++ )
								  X3[i] = m_pfXi[i];
							fP3 = m_fPress;
					 }
					 else
					 {
							for( i=0; i<m_iNum+1; i++ )
							{
								  X1[i] = X2[i];
								  X2[i] = X3[i];
								  X3[i] = m_pfXi[i];
							}
							fP1 = fP2;
							fP2 = fP3;
							fP3 = m_fPress;
					 }
					 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 = 0;
					 break;     //成功
				}
		  }
	 }
	 if( iIterNum == iMaxIterNum && iResFlag != 0 )
		  iResFlag = 10;

	 delete PP;
	 delete U;
	 delete Y1;
	 for( i=0; i<m_iNum+1; i++ )
			free( H[i] );
	 free( H );

	 return iResFlag;	    //不成功
}

//==========================设置拟压力计算的压力段数据=========================
void EOSFlashClass::SetPMP_Pressure( )
{
	float   fPressure=(float)0.05, fStep;
	double  fDewPressure = m_fDewPressure * 0.101325;
	
	m_pfPMP_Press[0] = fPressure;

	for( int i=1; i<m_iPMP_Num; 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;
		if( fPressure >= fDewPressure && m_pfPMP_Press[i-1] < fDewPressure-0.05 )
			m_pfPMP_Press[i] = float(fDewPressure);
		else
			m_pfPMP_Press[i] = fPressure;
	}
}



int EOSFlashClass::FlashCalculate( float* pfComposition, float fTemp, float fPress )
{   //摩尔组分
	/*pfComposition[0] = float(0.807);     //C1
	pfComposition[1] = (float)0.091;     //C2
	pfComposition[2] = (float)0.0179;    //C3
	pfComposition[3] = (float)0.0078;    //iC4
	pfComposition[4] = (float)0.005;     //nC4
	pfComposition[5] = (float)0.0008;    //iC5
	pfComposition[6] = 0;                //nC5
	pfComposition[7] = (float)0.0555;    //C6
	pfComposition[8] = (float)0.005;     //C7+
	pfComposition[9] = (float)0.005;     //H2
	pfComposition[10] = (float)0.005;    //H2O
	pfComposition[11] = 0;               //H2S
	pfComposition[12] = 0;               //N2
	pfComposition[13] = 0;               //O2
	pfComposition[14] = 0;               //CO
	pfComposition[15] = 0;               //CO2
	pfComposition[16] = 0;               //He
	pfComposition[17] = 0;               //Air*/



	//设置初值
	m_fTemp  = fTemp;
	m_fPress = fPress;
	m_fPStep = 5;
	m_fV     = (float)0.999;
	SetComposition( pfComposition );
	InputCompFixData( pfComposition );

	//计算露点压力
	char pMsg[128];
	int  iRes = CalDewPressure(  );
	sprintf( pMsg, "露点压力 = %f  iRes=%d", m_fDewPressure, iRes );
	//AfxMessageBox( pMsg );
	if( iRes )
	{
		AfxMessageBox( "计算露点压力错误" );
		return 1;
	}

	SetPMP_Pressure( );


	//高于露点压力------------单相区
	//根据总组分计算临界压力、临界温度、分子量
	float fTc = 0, fPc = 0, fMg = 0;
	for( int i=0; i<m_iNum; i++ )
	{
		m_pfYi[i] = m_pfZi[i];
		m_pfXi[i] = 0;
		fTc = fTc + m_pfYi[i]*m_pfTci[i];
		fPc = fPc + m_pfYi[i]*m_pfPci[i];
		fMg = fMg + m_pfYi[i]*m_pfMi[i];
	}
	sprintf( pMsg, "单相期 Pc=%f  Tc=%f  rg=%f", fPc, fTc, fMg/29.8 );
	//AfxMessageBox( pMsg );

	
	float  fU  = 0, fDewPressure = float(m_fDewPressure*0.101325);
	i = m_iPMP_Num - 1;
	while( m_pfPMP_Press[i] >= fDewPressure )
	{
		m_fPress = float(m_pfPMP_Press[i]*9.86923);
	    //单相计算气体物性参数，忽略油相

		m_pfMg[i] = fMg;            //气相分子量
		m_pfZg[i] = Zg_Calculate( fTc, float(fPc/9.86923), 
		                                    m_fTemp,
										    float(m_pfPMP_Press[i]),
										    ZG_DEFAULT );
		m_pfDg[i] = float(m_pfPMP_Press[i]*fMg / (m_pfZg[i]*0.008314*m_fTemp));
		m_pfUg[i] = Ug_Lee_Gonzalez_Eakin( float(m_pfPMP_Press[i]), m_fTemp, float(fMg/28.96), m_pfZg[i] );
		m_pfUg[i] = Ug_Calculate( fTc, float(fPc/9.86923), 
		                                    m_fTemp,
										    float(m_pfPMP_Press[i]),
										    m_pfZg[i],
										    float(fMg/28.96),
										    0,//pWellResFluidParm->GetN2()/100,
										    0,//pWellResFluidParm->GetCO2()/100,
										    0,//pWellResFluidParm->GetH2S()/100,
										    UG_DEFAULT );

		
		m_pfV[i]  = 1;             //气相摩尔分数
		m_pfSg[i] = 1;            //气相饱和度
		m_pfDo[i] = m_pfDg[i];            //油相密度
		m_pfZo[i] = m_pfZg[i];            //油相偏差因子
		m_pfUo[i] = m_pfUg[i];            //油相粘度
		m_pfMo[i] = m_pfMg[i];            //油相分子量
		i--;
	}

	//低于露点压力------------两相区
	if( i >= 0 )
		m_fPress = float(m_pfPMP_Press[i]*9.86923);
	if( m_fPress < m_fDewPressure )
	{
		m_fV = (float)0.9999;
		iRes = 1;
		while( iRes )
		{
			iRes = CalZV(  );
			m_fPress -= 0.5;
			if( m_fPress < 5 )
				break;
		}
		if( iRes )
		{
			AfxMessageBox( "计算气相Mol分数错误" );
			return 1;
		}
		while( m_fPress > 5 )
		{
		    iRes = BroyDen(  );
			if( iRes )
			{
			    sprintf( pMsg,"拟牛顿迭代错误  iRes=%d", iRes );
				AfxMessageBox( pMsg );
				break;
			}


			/*
			//将气相和液相分别按单相气体处理，计算物性参数
			fTc = 0;
			fPc = 0;
			fMg = 0;
			for( i=0; i<pCoeffData->iNum; i++ )
			{
				fTc = fTc + pFlashData->m_pfYi[i]*pCoeffData->m_pfTci[i];
				fPc = fPc + pFlashData->m_pfYi[i]*pCoeffData->m_pfPci[i];
				fMg = fMg + pFlashData->m_pfYi[i]*pCoeffData->m_pfMi[i];
			}
			pFlashData->fZg = Zg_Calculate( fTc, fPc/9.86923, 
		                                    pFlashData->fTemp,
										    pFlashData->fPress/9.86923,
											ZG_DEFAULT );
			pFlashData->fUg = Ug_Calculate( fTc, fPc/9.86923, 
		                                  pFlashData->fTemp,
										  pFlashData->fPress/9.86923,
										  pFlashData->fZg,
										  fMg/28.96,
										  pFlashData->m_pfZi[COMP_H2S],
										  pFlashData->m_pfZi[COMP_CO2],
										  pFlashData->m_pfZi[COMP_N2],
										  UG_DEFAULT );
			fCg = Cg_Calculate( fTc, fPc/9.86, pFlashData->fTemp, pFlashData->fPress/9.86, pFlashData->fZg, CG_DEFAULT );
			fBg = Bg_Calculate( pFlashData->fTemp, pFlashData->fPress/9.86923, pFlashData->fZg );


			fTc = 0;
			fPc = 0;
			fMg = 0;
			for( i=0; i<pCoeffData->iNum; i++ )
			{
				fTc = fTc + pFlashData->m_pfXi[i]*pCoeffData->m_pfTci[i];
				fPc = fPc + pFlashData->m_pfXi[i]*pCoeffData->m_pfPci[i];
				fMg = fMg + pFlashData->m_pfXi[i]*pCoeffData->m_pfMi[i];
			}
			pFlashData->fUo = Ug_Calculate( fTc, fPc/9.86923, 
		                                  pFlashData->fTemp,
										  pFlashData->fPress/9.86923,
										  pFlashData->fZo,
										  fMg/28.96,
										  pFlashData->m_pfZi[COMP_H2S],
										  pFlashData->m_pfZi[COMP_CO2],
										  pFlashData->m_pfZi[COMP_N2],
										  UG_DEFAULT );
			fCo = Cg_Calculate( fTc, fPc/9.86, pFlashData->fTemp, pFlashData->fPress/9.86, pFlashData->fZo, CG_DEFAULT );
			*/

			CalMoleDensity(  );
			CalViscosity( m_pfXi, m_fDo, fU );
			m_fUo = fU;
			CalViscosity( m_pfYi, m_fDg, fU );
			m_fUg = fU;
			m_pfUo[i] = m_fUo;            //油相粘度
			m_pfUg[i] = m_fUg;            //气相粘度
			m_pfZo[i] = m_fZo;            //油相偏差因子
			m_pfZg[i] = m_fZg;            //气相偏差因子


			CalMolecularWeight(  );
			m_pfMo[i] = m_fMo*1000;            //油相分子量
			m_pfMg[i] = m_fMg*1000;            //气相分子量

			CalQualityDensity( );
			m_pfDo[i] = m_fDo;            //油相密度
			m_pfDg[i] = m_fDg;            //气相密度

			CalVolumeFraction( );
			m_pfSg[i] = m_fSg;            //气相饱和度

			if( m_fV > 1 )
				m_fV = 1;
			m_pfV[i]  = m_fV;             //气相摩尔分数

			i--;
			m_fPStep = float((m_pfPMP_Press[i+1]-m_pfPMP_Press[i])*9.86923);
			InitNextFlashSpot(  );
			//AfxMessageBox( "下一次拟牛顿迭代" );
		}
	}
	return 0;
}