📄 cal_zsqj.cpp
字号:
// Cal_Zsqj.cpp: implementation of the CCal_Zsqj class.
//
//////////////////////////////////////////////////////////////////////
#include "stdafx.h"
#include "cvenus.h"
#include "Cal_Zsqj.h"
#include <math.h>
#ifdef _DEBUG
#undef THIS_FILE
static char THIS_FILE[]=__FILE__;
#define new DEBUG_NEW
#endif
extern int g_iProgress ; // 计算进度
extern bool bFinishThread;
extern char cGraphColor[][16];
extern int iGraphColorNum;
//////////////////////////////////////////////////////////////////////
// Construction/Destruction
//////////////////////////////////////////////////////////////////////
#define DAMP0 110
#define KMAX 2000
#define MAIR 28.96
#define R 8314
#define KC 2001
#define EC 14.6
#define PI 3.1416
#define PN 101325
#define TN 298.15
#define G 9.81
CCal_Zsqj::CCal_Zsqj()
{
/*
DAMP0 = 110;
KMAX = 2000;
MAIR = 28.96;
R = 8314;
KC = 2001;
EC = 14.6;
PI = 3.1416;
PN = 101325;
TN = 298.15;
G = 9.81;
*/
}
CCal_Zsqj::~CCal_Zsqj()
{
}
// 举升分析
void CCal_Zsqj::Jsfx()
{
double DIN[21], Dout[12];
bool Dout_12;
CF1Book* m_pGrid;
m_pGrid = &(((CCVenusApp*)AfxGetApp())->m_pMainDlg->m_grid);
int i;
CString str;
for(i=1; i<21; i++)
{
str = m_pGrid->GetTextRC(i, 5);
DIN[i] = atof(str);
}
g_iProgress = 20; // 计算进度
ENDCOMPUTE;
// 开始计算,Sub MMMain (DIN!(), Dout!())
double DWT, H2S, CO2, vw, PCAS, TAV;
H = DIN[1];
HC = DIN[2];
DC = DIN[3];
DTO = DIN[4];
DTI = DIN[5];
DWT = DIN[6];
PS = DIN[7];
DNTL = DIN[8];
MUL = DIN[9];
RG = DIN[10];
H2S = DIN[11];
CO2 = DIN[12];
CPV = DIN[13];
TMPS = DIN[14];
TMPB = TMPS + DWT * H;
TMPS = 273 + TMPS;
TMPB = 273 + TMPB;
PTMIN = DIN[15];
GLR = DIN[16];
QWD = DIN[17];
PM = DIN[18];
vw = DIN[19];
PCAS = DIN[20];
M = RG * MAIR;
double C1 = 6890;
double C2 = 10000;
AC = .25 * PI * (DC * DC - DTO * DTO);
AT = .25 * PI * DTI * DTI;
AK = .25 * PI * DD * DD;
ACH = .25 * PI * DCH * DCH;
TAV = 273.15 + .5 * (TMPS + TMPB);
PS = 1000000 * PS;
PTMIN = 1000000 * PTMIN;
double HLI, PLH, VUP, VDL, VDG, N1, FF, PLF, PHC0, PHC1, CR1;
Dout[1] = vw;
HLI = vw / AT;
DAMP = DAMP0 + 4.154 * HLI;
PLH = DNTL * G * HLI;
VUP = 5.08;
VDG = 0.85;
VDL = 0.2;
N1 = HC / VUP + (HC - HLI) / VDG + HLI / VDL;
N1 = 86400 / N1;
FF = MODY(DNTL, DTI, VUP, MUL);
PLF = .5 * FF * HLI * DNTL * VUP * VUP / DTI;
PHC0 = PTMIN + PLH + PLF + PM * G / AT;
PHC1 = PHC0;
CR1 = 1;
double PAV, PGHT, PGFT, ZAVT, UG, DNTGT, FG, PHC2;
PGHT = PGFT = 0;
while(CR1 > .5)
{
PAV = PHC0 + .5 * (PGHT + PGFT);
ZAVT = ZFACTR(PAV, TAV); // Gas Divation Factor in Tub.
UG = VISCO(ZAVT, PAV, TAV); // Gas Viscosity in Tub.
DNTGT = M * PAV / (R * ZAVT * TAV); // Gas Density in Tub.
FG = JAIN(DNTGT, DTI, VUP, UG); // Damping Factor in Tub.
PGHT = DNTGT * G * (HC - HLI);
PGFT = .5 * FG * (HC - HLI) * DNTGT * VUP * VUP / DTI;
PHC2 = PHC0 + PGFT + PGHT; // Pressure at the
if( fabs((PHC1 - PHC2) / PHC2) > .00001)
PHC1 = PHC2;
else
CR1 = 0;
}
double PHC, PGHC, PCMIN1, CR2, ZAVC, DNTGC, DEQ, PCMIN2;
PHC = .5 * (PHC2 + PHC1); // Pressure at the
PGHC = DNTGT * G * HC; // Initial Hydro-gas Pressure at Casing
PCMIN1 = PHC - PGHC; // Minimum Pressure at Cas.
CR2 = 1;
while(CR2 > .5)
{
PAV = .5 * (PHC + PCMIN1); // Average Pressure
ZAVC = ZFACTR(PAV, TAV); // Gas Deviation Factor in Cas.
UG = VISCO(ZAVC, PAV, TAV); // Gas Viscosity in Cas.
DNTGC = M * PAV / (R * ZAVC * TAV); // Gas Density in Cas.
DEQ = DC - DTO; // Equverlent Dia.
FG = JAIN(DNTGC, DEQ, VUP, UG); // Damping Factor
PGHC = DNTGC * G * HC;
PCMIN2 = PHC - PGHC; // Minimum Pressure at Cas.
if( fabs((PCMIN2 - PCMIN1) / PCMIN2) > .0001)
PCMIN1 = PCMIN2;
else
CR2 = 0;
}
double PCMIN, PCMAX, FA, PTMAX, VG;
PCMIN = .5 * (PCMIN1 + PCMIN2);
PCMIN = 1.1 * PCMIN;
PCMAX = (1 + AT / AC) * PCMIN;
FA = 1.293 * G * TN * RG / (PN * TAV * ZAVC);
PTMAX = PCMAX * exp(FA * HC) - DNTL * G * HLI;
PTMAX = PTMAX / exp(FA * (HC - HLI));
VG = PCMAX * TN * AT * (HC - HLI);
VG = VG * 1.15 * (1 + .02 * HC / 305);
VG = VG / (PN * TAV * ZAVT);
VG = int(VG);
double QG, MGAS, TU, TDG, TDL, TPL, TOP, N2;
// -- Dynamic Model----------------------
// MGAS = TUP(PCMAX, PCMIN, PTMAX, PTMIN, HLI, TU);
// HLI = TDOWN(PCMAX, PCMIN, PTMAX, PTMIN, TDG, TDL);
TDL = TDG = TU = MGAS = 0;
QG = 1.25 * MGAS / (1.293 * RG); //Modified Factor 1.25
QG = int(QG + .5);
TU = int(TU + .5);
TDG = int(TDG + .5);
TDL = int(TDL + .5);
double DGL = 0;
TPL = TU + TDG + DGL;
TOP = 15 + 5 * vw;
TPL = 60 * TOP + (HC - HLI) / VDG + HLI / VDL;
N2 = int(864000 / TPL + 5);
N2 = N2 / 10;
double PCAVG, DTAP, QPI, P1, P2, P3, P4, TF, N, WD;
PCAVG = .5 * (PCMAX + PCMIN) + DNTGC * G * HC + DNTL * G * (H - HC);
DTAP = 1000000;
QPI = QWD / DTAP / 86400;
P1 = int(PCMAX / 1000) / 1000.0;
Dout[2] = P1;
P2 = int(PCMIN / 1000) / 1000.0;
P3 = int(PTMAX / 1000) / 1000.0;
Dout[3] = P3;
P4 = .5 * (P1 + P2);
TF = vw / (QPI * (PS - PCAVG));
TF = int(TF);
N = int(864000 / TF + 5);
N = N / 10;
if(N > N2) N = N2;
if((N - int(N)) > .7)
N = int(N) + 1;
else if(N - int(N) < .3)
N = int(N);
else
N = int(N) + .5;
g_iProgress = 50; // 计算进度
ENDCOMPUTE;
double VIN, TDN, GDD, GDF, QIN;
WD = N * vw;
VIN = int(VG / vw - GLR);
if( N <= 0 )
{
N = 0;
TOP = 0;
TDN = 0;
}
else
TDN = int(1440 / N) - TOP;
WD = N * vw;
VIN = int(VG / vw - GLR);
WD = N * vw;
GDD = int(N * GLR * vw) / C2;
GDF = int(N * VG) / C2;
QIN = N * VIN * vw / C2;
Dout[4] = N;
Dout[5] = int(VG / vw);
Dout[6] = VIN;
Dout[7] = WD;
Dout[8] = GDD;
Dout[9] = QIN;
Dout[10] = TOP;
Dout[11] = TDN;
if( PCAS >= Dout[2])
Dout_12 = true;
else
Dout_12 = false;
g_iProgress = 80; // 计算进度
ENDCOMPUTE;
// 打开保存计算结果的文件,保存计算结果
CStdioFile f3;
CFileException exception3;
BOOL status1;
CString templs1;
templs1 = ((CCVenusApp*)AfxGetApp())->GetMyAppPath();
templs1 += "\\tempfile.txt";
status1 = f3.Open(templs1,CFile::modeCreate|CFile::modeWrite);
if (!status1)
{
char s[100];
sprintf(s, "Error opening file for reading. Code:%d",
exception3.m_cause);
AfxMessageBox(s);
return;
}
try
{
char* cOutput_str[12]
= {"产水量\tm^3/cycle\t",
"最大套压\tMPa\t",
"最小油压\tMPa\t",
"举升周期数\tcycle/d\t",
"需气水比\tm^3/m^3\t",
"周期注气量\tm^3/cycle\t",
"日排水量\tm^3/d\t",
"日产气量\t10^4m^3/d\t",
"日注气量\t10^4m^3/d\t",
"柱塞上行\tmin\t",
"柱塞下行\tmin\t",
"分析结果\t$16\t"
};
CString kg="\n";
f3.WriteString(kg);
templs1="$$柱塞气举分析计算结果$0\n";
f3.WriteString(templs1);
f3.WriteString(kg);
for(int i=0;i<11;i++)
{
templs1.Format("%f\n", Dout[i+1]);
templs1 = cOutput_str[i] + templs1;
f3.WriteString(templs1);
}
templs1 = Dout[12]?"能够举升":"不能举升";
templs1 = cOutput_str[11] + templs1;
f3.WriteString(templs1);
}
catch (CFileException exception3)
{
AfxMessageBox("Error reading file");
return;
}
f3.Close();
}
// 优化设计
void CCal_Zsqj::Yhsj()
{
int i;
double DIN[19], Dout[10][12];
BOOL bDtmx;
int Dout_12; // 推荐方案
for(i=0;i<10;i++)
{
for(int j=0;j<12;j++)
Dout[i][j]=0;
}
CF1Book* m_pGrid;
m_pGrid = &(((CCVenusApp*)AfxGetApp())->m_pMainDlg->m_grid);
CString str;
for(i=1; i<19; i++)
{
str = m_pGrid->GetTextRC(i, 5);
DIN[i] = atof(str);
}
bDtmx = ((CButton*)((CCVenusApp*)AfxGetApp())->m_pMainDlg->m_Page18.GetDlgItem(IDC_CHECK1))->GetCheck();
g_iProgress = 10; // 计算进度
ENDCOMPUTE;
// 开始计算,Sub MMAin (DIN!(), Dout!())
double DWT, H2S, CO2;
H = DIN[1];
HC = DIN[2];
DC = DIN[3];
DTO = DIN[4];
DTI = DIN[5];
DWT = DIN[6];
PS = DIN[7];
DNTL = DIN[8];
MUL = DIN[9];
RG = DIN[10];
H2S = DIN[11];
CO2 = DIN[12];
CPV = DIN[13];
TMPS = DIN[14];
TMPB = TMPS + DWT * H;
TMPS = 273 + TMPS;
TMPB = 273 + TMPB;
PTMIN = DIN[15];
GLR = DIN[16];
QWD = DIN[17];
PM = DIN[18];
double C1, C2, TAV;
M = RG * MAIR;
C1 = 6890;
C2 = 10000;
AC = .25 * PI * (DC * DC - DTO * DTO);
AT = .25 * PI * DTI * DTI;
AK = .25 * PI * DD * DD;
ACH = .25 * PI * DCH * DCH;
TAV = 273.15 + .5 * (TMPS + TMPB);
PS = 1000000 * PS;
PTMIN = 1000000 * PTMIN;
// -- Foss & Gaul Method ---
int IBEGIN=6, IEND=14, ISAVE;
double vw, HLI, PLH, VUP, VDG, VDL, N1, FF, PLF, PHC0, PHC1, CR1;
double PAV, PGHT=0, PGFT=0, ZAVT, UG, DNTGT, FG, PHC2;
double PHC, PGHC, PCMIN1, CR2, ZAVC, DNTGC, DEQ, PCMIN2;
double PCMIN, PCMAX, FA, PTMAX, VG;
double TU=0, MGAS=0, TDG=0, TDL=0;
double QG=0, TPL=0, DGL=0, TOP=0, N2=0, PCAVG=0, DTAP=0, QPI=0;
double P1=0, P2=0, P3=0, P4=0, TF=0, N=0,WD=0, VIN=0, TDN=0, GDD=0;
double GDF=0, QIN=0;
double dQujian;
dQujian = 90.0 / (IEND - IBEGIN + 1);
for(int I = IBEGIN; I <= IEND ; I++)
{
vw = I * .1;
ISAVE = I - IBEGIN + 1;
Dout[ISAVE][1] = vw;
HLI = vw / AT;
DAMP = DAMP0 + 4.154 * HLI;
PLH = DNTL * G * HLI;
VUP = 5.08;
VDG = .85;
VDL = .2;
N1 = HC / VUP + (HC - HLI) / VDG + HLI / VDL;
N1 = 86400 / N1;
FF = MODY(DNTL, DTI, VUP, MUL);
PLF = .5 * FF * HLI * DNTL * VUP * VUP / DTI;
PHC0 = PTMIN + PLH + PLF + PM * G / AT;
PHC1 = PHC0;
CR1 = 1;
// LOOP the TUBINING
while( CR1 > .5)
{
PAV = PHC0 + .5 * (PGHT + PGFT);
ZAVT = ZFACTR(PAV, TAV); // Gas Divation Factor in Tub.
UG = VISCO(ZAVT, PAV, TAV); // Gas Viscosity in Tub.
DNTGT = M * PAV / (R * ZAVT * TAV); // Gas Density in Tub.
FG = JAIN(DNTGT, DTI, VUP, UG); // Damping Factor in Tub.
PGHT = DNTGT * G * (HC - HLI);
PGFT = .5 * FG * (HC - HLI) * DNTGT * VUP * VUP / DTI;
PHC2 = PHC0 + PGFT + PGHT; // Pressure at the
if (fabs((PHC1 - PHC2) / PHC2) > .00001)
PHC1 = PHC2;
else
CR1 = 0;
}
PHC = .5 * (PHC2 + PHC1); // Pressure at the
PGHC = DNTGT * G * HC; // Initial Hydro-gas Pressure at Casing
PCMIN1 = PHC - PGHC; // Minimum Pressure at Cas.
CR2 = 1;
while( CR2 > .5)
{
PAV = .5 * (PHC + PCMIN1); // Average Pressure
ZAVC = ZFACTR(PAV, TAV); // Gas Deviation Factor in Cas.
UG = VISCO(ZAVC, PAV, TAV); // Gas Viscosity in Cas.
DNTGC = M * PAV / (R * ZAVC * TAV); // Gas Density in Cas.
DEQ = DC - DTO; // quverlent Dia.
FG = JAIN(DNTGC, DEQ, VUP, UG); // Damping Factor
PGHC = DNTGC * G * HC;
PCMIN2 = PHC - PGHC; // Minimum Pressure at Cas.
if(fabs((PCMIN2 - PCMIN1) / PCMIN2) > .0001)
PCMIN1 = PCMIN2;
else
CR2 = 0;
}
PCMIN = .5 * (PCMIN1 + PCMIN2);
PCMIN = 1.1 * PCMIN;
PCMAX = (1 + AT / AC) * PCMIN;
FA = 1.293 * G * TN * RG / (PN * TAV * ZAVC);
// PCAVG = .5 * (PHC + PCMAX * EXP(FA * HC));
PTMAX = PCMAX * exp(FA * HC) - DNTL * G * HLI;
PTMAX = PTMAX / exp(FA * (HC - HLI));
VG = PCMAX * TN * AT * (HC);
VG = VG * 1.25 * (1 + .02 * HC / 305);
VG = VG / (PN * TAV * ZAVT);
VG = int(VG + 150);
// -- Dynamic Model----------------------
if (bDtmx)
{
TUP(PCMAX, PCMIN, PTMAX, PTMIN, HLI, TU, MGAS);
TDOWN(PCMAX, PCMIN, PTMAX, PTMIN, TDG, TDL, HLI);
}
QG = 1.25 * MGAS / (1.293 * RG); // Modified Factor 1.25
QG = int(QG + .5);
TU = int(TU + .5);
TDG = int(TDG + .5);
TDL = int(TDL + .5);
TPL = TU + TDG + DGL;
TOP = 15 + 5 * vw;
TPL = 60 * TOP + (HC - HLI) / VDG + HLI / VDL;
N2 = int(864000 / TPL + 5);
N2 = N2 / 10;
PCAVG = .5 * (PCMAX + PCMIN) + DNTGC * G * HC + DNTL * G * (H - HC);
DTAP = 1000000;
QPI = QWD / DTAP / 86400;
P1 = int(PCMAX / 1000) / 1000.0;
Dout[ISAVE][2] = P1;
P2 = int(PCMIN / 1000) / 1000.0;
P3 = int(PTMAX / 1000) / 1000.0;
Dout[ISAVE][3] = P3;
P4 = .5 * (P1 + P2);
TF = vw / (QPI * (PS - PCAVG));
TF = int(TF);
N = int(864000 / TF + 5);
N = N / 10;
if( N > N2) N = N2;
if (N - int(N) > .7)
N = int(N) + 1;
else if (N - int(N) < .3)
N = int(N);
else
N = int(N) + .5;
WD = N * vw;
VIN = int(VG / vw - GLR);
if( N <= 0 )
{
N = 0;
TOP = 0;
TDN = 0;
}
else
TDN = int(1440 / N) - TOP;
WD = N * vw;
// VIN = int(VG / vw - GLR);
// 这里与VB计算的有差异,
// 因为这里VIN的计算结果是负数,假如是-556.4,同样是int函数
// VB与C的计算结果不同,VB是-557,而C的计算结果是-556。
// 所以计算结果中的6项有差异,这又直接导致第9项QIN有差异
// 为了弥补此项差异,更改如下
VIN = int(VG / vw - GLR - 1);
WD = N * vw;
GDD = int(N * GLR * vw) / C2;
GDF = int(N * VG) / C2;
QIN = N * VIN * vw / C2;
Dout[ISAVE][4] = N;
Dout[ISAVE][5] = int(VG / vw);
Dout[ISAVE][6] = VIN;
Dout[ISAVE][7] = WD;
Dout[ISAVE][8] = GDD;
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -