qjyljs.cpp

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

// qjyljs.cpp : implementation file
//

#include "stdafx.h"
#include "qjyljs.h"

#include "yqsxz.h"
#include "math.h"

#ifdef _DEBUG
#define new DEBUG_NEW
#undef THIS_FILE
static char THIS_FILE[] = __FILE__;
#endif

/////////////////////////////////////////////////////////////////////////////
// qjyljs
extern double ptidu[1000],dtidu[1000];      
qjyljs::qjyljs()
{
}

qjyljs::~qjyljs()
{
}



/////////////////////////////////////////////////////////////////////////////
// qjyljs message handlers

double qjyljs::qjyl(double P0,double T0,double TG,double L,double DOL,double DIL,double RR,double RG,double GG,double LGR,double GL,double H2S,double N2,double CO2,double NaCl,int NP,int CD,int FD)//气井压力
{     //Sub 气井压力(Pcom, P0,       T0,       TG,       L,       DOL,       DIL,       RR,       RG,       GG,       LGR,       GL,       H2S,       N2,       CO2,                NP,    CD,    FD)
	  //P起点      [[[[  pqidian
	  //T起点      [[[[  tqidian
      //T梯度      [[[[  ttidu
	  //Q气        [[[[  qqi
      //G液        [[[[  gye
	  //G气        [[[[  gqi
	  //L管        [[[[  lguan
	  //D环管外    [[[[  dhgw
	  //D环管内    [[[[  dhgn
	  //RR管       [[[[  rrguan
	  //N管段      [[[[  nguanduan
	  //D计算      [[[[  djisuan
      //D流动      [[[[  dliudong
	  //P梯度(3000)    [[[[  ptidu[3000]
	  //D梯度      [[[[  dtidu[3000]
	  //D管        [[[[  dguan
	  //Q气1       [[[[  qqi1
	  //G气1       [[[[  gqi1

	  double Pcom,pqidian,tqidian,ttidu,qqi,LGRatio,gye,gqi;
      double lguan,dhgw,dhgn,rrguan,YN2,YCO2,YH2S,YNaCl;
	  double djisuan,dliudong,dguan,PC,TC,GLR,ML,GE,P2;
	  double DLTD,qqi1,gqi1,T,pold,P;
	  double PR,TR,Z,UG,RE,D5,PF,XY,FF;
	  int nguanduan,i,KK;
	  yqsxz yqsxz;

	  pqidian = P0;
	  tqidian = T0;
	  ttidu = TG;
	  qqi = RG;
	  LGRatio = LGR;
	  gye = GL;
	  gqi = GG;
	  lguan = L;
	  dhgw = DOL;
	  dhgn = DIL;
	  rrguan = RR;
	  YN2 = N2;
	  YCO2 = CO2;
	  YH2S = H2S;
	  YNaCl = NaCl;
	  nguanduan = NP; //流动管分段数
	  djisuan = CD; //从上往下为1，从下往上为-1
	  dliudong = FD; //从上往下为-1，从下往上为1
	  if( nguanduan <= 0)
	  {
        ldiv_t lguanshang;
        lguanshang=ldiv(long(lguan),500);
		if(lguanshang.quot== lguan/50)
		{
			nguanduan = lguanshang.quot;
		}
		else
		{
			nguanduan = lguanshang.quot + 1;
		}
	  }
	  else
	  {
		for(i=1;i<=999;i++)
		{
		  ptidu[i] = 0;
		  dtidu[i] = 0;
		}
	  }
	  if( nguanduan <= 2)
			nguanduan=2;
	  DLTD = lguan / nguanduan;
	  dguan = dhgw - dhgn;
	  if( rrguan <= 0)
	  {
		  rrguan = 0.0008;
	  }
      yqsxz.trqnjcs(GG, H2S, CO2,1);//临界参数
	  PC=yqsxz.Pc;//对比压力
	  TC=yqsxz.Tc;//对比温度

	  if( DLTD > 0)
	  {
		if( LGRatio > 0)
		{
		  GLR = 1 / LGRatio;
		  if( gye == 0)
		  {
			  gye = 1;
		  }
		  ML = 6084 / (141.5 / gye - 137.4);
		  GE = 8.3093 * gye / ML * 293 / 0.1015;
		  qqi1 = qqi * (1 + GE / GLR);
		  gqi1 = (GLR * gqi + 815.95 * gye) / (GLR + 23638.39 * gye / ML);
		}
		else
		{
		  qqi1 = qqi;
		  gqi1 = gqi;
		}
		if( qqi1 == 0)
		{
			qqi1 = pow(10,-20);
		}
		if( rrguan == 0)
		{
			rrguan = 0.0008;
		}
		P2 = pqidian;
	    for(i=1;i<=nguanduan;i++)
		{
		  pqidian = P2;
		  P2 = pqidian + 0.5 * djisuan;
		  T = tqidian + ttidu * DLTD / 2 * djisuan;
		  pold = P2;
		  for(KK=1;KK<=50;KK++)
		  {
		  	P2 = (P2 + pold) / 2;
			pold = P2;
			P = (pqidian + P2) / 2;
			if( P < 0)
			{
				P = 0.001;
			}
			PR = P / PC;
			TR = T / TC;
	        Z=yqsxz.trqpcxs(PR, TR);//天然气偏差因子
			if( qqi > 0.01)
			{
              UG=yqsxz.trqnd(PR, TR, YCO2, YH2S, YN2, GG,T-273.15);  	
			  RE = 0.0178 * gqi1 * qqi1 / dguan / UG;
			  FF=yqsxz.mczlxs(RE, rrguan);
			}
			else
			{
			  FF = 0;
			}
			XY = 0.0683 * gqi1 * DLTD / Z / T;
			D5 = (dhgw - dhgn) * pow((pow(dhgw,2) - pow(dhgn,2)),2);
			PF = (1.324E-18 * FF * pow((qqi1 * T * Z),2) * (exp(XY) - 1)) / D5;
			if( dliudong == 1)
			{
			  if( djisuan == 1)
			  {
				P2 = (pow(pqidian,2) * exp(XY) + PF);
			  }
			  else
			  {
				P2 = (pow(pqidian,2) - PF) / exp(XY);
			  }
			}
			else
			{
			  if( djisuan == 1)
			  {
				P2 = (pow(pqidian,2) * exp(XY) - PF);
			  }
			  else
			  {
				P2 = (pow(pqidian,2) + PF) / exp(XY);
			  }
			}
			if( P2 > 0)
			{
			  P2 = pow(P2,0.5);
			}
			else
			{
			  P2 = 0;
			  break;
			}
			if(fabs(pold - P2) < 0.005)
			{
				break;
			}
		  }
		  ptidu[i] = int(1000 * P2 + 0.5);
		  ptidu[i]=ptidu[i]/ 1000;
		  dtidu[i] = dtidu[i - 1] + djisuan * DLTD;
		  tqidian = tqidian + ttidu * DLTD * djisuan;
		  if( P2 == 0)
		  {
			  break;
		  }
		}
		Pcom = P2;
	  }
	  else
	  {
		Pcom = pqidian;
	  }         
	  return (Pcom);
}