📄 cal_zsqj.cpp
字号:
// Dout[ISAVE][9] = QIN;
// Change to, 保留三位小数点
Dout[ISAVE][9] = int(QIN*1000)/1000.0;
Dout[ISAVE][10] = TOP;
Dout[ISAVE][11] = TDN;
g_iProgress = 10 + int(ISAVE*dQujian); // 计算进度
ENDCOMPUTE;
}
I = 2;
while( Dout[I][7] >= Dout[I - 1][7] && I < 9)
I++;
Dout_12 = I - 1;
g_iProgress = 95; // 计算进度
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;
}
CString strTemp;
try
{
CString kg="\n";
f3.WriteString(kg);
templs1="\t\t\t$$柱塞气举优化设计$0\n";
f3.WriteString(templs1);
f3.WriteString(kg);
templs1="\t\t\t\t[设计计算结果]$0\n";
f3.WriteString(templs1);
char* cTitle1=" 产水量$7\t最大套压$7\t最小油压$7\t举升周期数$7\t需气水比$7\t"
"日排水量$7\t 日产气量$7\t 日注气量$7\t柱塞上行$7\t柱塞下行$7\n";
f3.WriteString(cTitle1);
char* cTitle2="m^3/cycle$11\t MPa$11\t MPa$11\t cycle/d$11\tm^3/m^3$11\t"
" m^3/d$11\t10^4 m^3/d$11\t10^4 m^3/d$11\t min$11\t min$11\n";
f3.WriteString(cTitle2);
for(int i=1;i<10;i++)
{
templs1.Empty();
for(int j=1; j<12; j++)
{
if(j==6) continue;
strTemp.Format("%f", Dout[i][j]);
templs1 += strTemp;
templs1 += j==11?"\n":"\t";
}
f3.WriteString(templs1);
}
f3.WriteString(kg);
templs1="\t\t\t\t[设计推荐方案]$0\n";
f3.WriteString(templs1);
f3.WriteString(cTitle1);
f3.WriteString(cTitle2);
templs1.Empty();
for(int j=1; j<12; j++)
{
if(j==6) continue;
strTemp.Format("%f", Dout[Dout_12][j]);
templs1 += strTemp;
templs1 += j==11?"\n":"\t";
}
f3.WriteString(templs1);
}
catch (CFileException exception3)
{
AfxMessageBox("Error reading file");
return;
}
f3.Close();
// 打开保存计算的图形结果的文件
CStdioFile f4;
CFileException exception4;
int iColorIndex = 0; // 颜色索引 lzj
templs1 = ((CCVenusApp*)AfxGetApp())->GetMyAppPath();
templs1 += "\\tempfile.gra";
status1 = f4.Open(templs1,CFile::modeCreate|CFile::modeWrite);
if (!status1)
{
char s[100];
sprintf(s, "Error opening file for reading. Code:%d",
exception4.m_cause);
AfxMessageBox(s);
return;
}
try
{
templs1="$$柱塞气举优化设计[产气和注气曲线]\n";
f4.WriteString(templs1);
f4.WriteString("气量最小值，10^3 m^3/d：\n");
f4.WriteString("气量最大值，10^3 m^3/d：\n");
f4.WriteString("周期产水量最小值，m^3/cycle：\n");
f4.WriteString("周期产水量最小值，m^3/cycle：\n");
f4.WriteString("气量\n");
f4.WriteString("周期产水量\n");
f4.WriteString("气量__10^3 m^3/d\n");
f4.WriteString("周期产水量__m^3/cycle\n");
f4.WriteString(" 1\n");
for(i=1;i<10;i++)
{
templs1.Format(" %f\t%f\n", Dout[i][1], Dout[i][8]);
f4.WriteString(templs1);
}
f4.WriteString("-1\t -1\n");
f4.WriteString(cGraphColor[iColorIndex++]);
iColorIndex %= iGraphColorNum;
f4.WriteString("产气曲线\n");
for(i=1;i<10;i++)
{
templs1.Format(" %f\t%f\n", Dout[i][1], Dout[i][9]);
f4.WriteString(templs1);
}
f4.WriteString("-1\t -1\n");
f4.WriteString(cGraphColor[iColorIndex++]);
iColorIndex %= iGraphColorNum;
f4.WriteString("注气曲线\n");
f4.WriteString("-10\t-10\n");
templs1="$$柱塞气举优化设计[产水曲线]\n";
f4.WriteString(templs1);
f4.WriteString("日排水量最小值，m^3/d：\n");
f4.WriteString("日排水量最大值，m^3/d：\n");
f4.WriteString("周期产水量最小值，m^3/cycle：\n");
f4.WriteString("周期产水量最小值，m^3/cycle：\n");
f4.WriteString("日排水量\n");
f4.WriteString("周期产水量\n");
f4.WriteString("日排水量__m^3/d\n");
f4.WriteString("周期产水量__m^3/cycle\n");
f4.WriteString(" 1\n");
for(i=1;i<10;i++)
{
templs1.Format(" %f\t%f\n", Dout[i][1], Dout[i][7]);
f4.WriteString(templs1);
}
f4.WriteString("-1\t -1\n");
f4.WriteString(cGraphColor[iColorIndex++]);
iColorIndex %= iGraphColorNum;
f4.WriteString("产水量曲线\n");
f4.WriteString("-10\t-10\n");
}
catch (CFileException exception4)
{
AfxMessageBox("Error reading file");
return;
}
f4.Close();
}
double CCal_Zsqj::JAIN(double RUO, double D, double V, double MU)
{
double FG;
double RE, E, A;
//-- Computing Gas's Damping Factor
RE = RUO * D * V / MU;
E = 0.00001524;
if( RE > 2300 )
{
A = 1.14 - 2 * log10(E / D + 21.25 / pow(RE, 0.9) );
FG = 1 / (A * A);
}
else if(RE > 64)
{
FG = 1 / RE;
}
else
FG = 0;
return FG;
}
// Fluid's Damping Factor
double CCal_Zsqj::MODY(double RUO, double D, double V, double MU)
{
double FF;
double RE;
RE = RUO * D * V / MU;
if( RE > 100000 )
FF = .0032 + .221 * pow(RE, -.237);
else if(RE > 2300)
FF = .3164 / sqrt(sqrt(RE));
else if( RE > 10 )
FF = 64 / RE;
else
FF = 0;
return FF;
}
double CCal_Zsqj::ZFACTR(double P, double T)
{
double Z;
double TPC, PPC, TPR, PPR, ZO, ROR;
M = 28.96 * RG;
TPC = 92.2 + 176.6 * RG;
if(RG >= .7)
PPC = 4.881 - .3861 * RG;
else
PPC = 4.778 - .248 * RG;
TPR = T / TPC;
PPR = .000001 * P / PPC;
ZO = 1;
while(1)
{
ROR = .27 * PPR / (ZO * TPR);
Z = .6815 * ROR * ROR / pow(TPR , 3);
Z = Z + (.5353 - .6123 / TPR) * ROR * ROR;
Z = 1 + Z + (.31506 - 1.0467 / TPR - .5783 / pow(TPR,3)) * ROR;
if(fabs(Z - ZO) <= .0001)
break;
else
ZO = Z;
}
Z = .94;
return Z;
}
double CCal_Zsqj::VISCO(double Z, double P, double T)
{
double UG;
double X, Y, D, DTG;
M = 29.86 * RG;
X = 3.5 + 986 / (1.8 * T) + .01 * M;
Y = 2.4 - .2 * X;
D = (9.4 + .02 * M) * pow(1.8 * T, 1.5);
D = D / (209 + 19 * M + 1.8 * T);
DTG = 3.4844 * RG * P / 1000000 / (Z * T);
UG = .0001 * D * exp(X * pow(DTG, Y));
UG = .016;
UG = UG / 1000;
return UG;
}
void CCal_Zsqj:: TUP (double& PCMAX, double& PCMIN, double& PTMAX,
double& PTMIN, double& HLI, double& TU, double& MGAS)
{
static double PTS[KC+1], PTSK[KMAX+1], XTG[KC+1], XTGK[KMAX+1];
static double PCS[KC+1], PCSK[KMAX+1], XCL[KC+1], XCLK[KMAX+1];
static double XTP[KC+1], XTPK[KMAX+1], XTK[KC+1], XTKK[KMAX+1];
static double XTL[KC+1], XTLK[KMAX+1], XCH[KC+1], XCHK[KMAX+1];
double TAV, PAV, ZAV, ZS, DNTG, FA, FB, FD, FT;
TAV = .5 * (TMPS + TMPB) + 273.15;
PAV = .5 * (PTMIN + PTMAX);
ZAV = ZFACTR(PAV, TAV);
ZS = 1;
DNTG = 1.293 * RG;
FA = G * TN * DNTG / (PN * TAV * ZAV);
FB = M * AT / (ZAV * R * TAV);
FD = DAMP / PM;;
FT = AT / PM;
PTSK[0] = PTMAX;
PCSK[0] = PCMAX;
XCLK[0] = 0;
// HLI * AT / (AT + AC);
XTGK[0] = HLI; // HC - 2569
XTPK[0] = 0;
XTKK[0] = 0;
XCHK[0] = 0;
PTS[0] = PTSK[0];
PCS[0] = PCSK[0];
XCL[0] = XCLK[0];
XTG[0] = XTGK[0];
XTP[0] = XTPK[0];
XTK[0] = XTKK[0];
XCH[0] = XCHK[0];
double DXTP, DXTK, DXCH, DPCS, MLTP, DT;
int K;
DXTP = 0, DXTK = 0; // Initial value of dXtp/dt
DXCH = 0, DPCS = 0;
MLTP = .02;
MGAS = 0;
TU = 0;
DT = .001, K = 1;
int I;
double KK, VO, PT, PC, CL, TP, TG, TK, CH, HL, PCG, PTG, PHC, PWF, DLAV=0;
double HLX=0, MGS=0, FF=0, FHU=0, FHD=0, UG=0, DTG=0, FG=0, FGU=0, FGD=0;
double PPG=0, DPHC=0, DXCL=0, DPPG, PPU=0, PPD=0, XTPR=0, KD=0, DXTG=0;
double MG=0, MLC=0, ML=0, DPCG=0, MGC=0, MGTP=0, DPTG=0, PPCG=0, PPTG=0;
double DP=0, DPPG1=0, DPPG2=0;
double ACG, ATP=0, ATG=0, ATK=0, BCH=0, AV=0, ASD=0, EPC, EPT, FP;
double TUPC1=0, TUPC2=0, TUPT1, TUPT2;
while(1)
{
if( TU > 30 ) DT = .004;
KK = 1 / DT;
VO = DXTP;
for(I=1; I<=KK;I++)
{
TU = TU + DT;
PT = PTS[I - 1];
PC = PCS[I - 1];
CL = XCL[I - 1];
TP = XTP[I - 1];
TG = XTG[I - 1];
TK = XTK[I - 1];
CH = XCH[I - 1];
HL = XTG[I - 1] - XTP[I - 1]; // FLUID'S LENGTH
PCG = PC * exp(FA * (HC - CL) / 2);
PTG = PT * exp(FA * (HC - TG) / 2);
if(PT <= PTMIN) PT = PTMIN;
PHC = PC * exp(FA * (HC - CL)) + DNTL * G * CL;
PWF = PHC + DLAV * G * (H - HC);
PREWF(PHC, PWF, ZAV, TAV);
// ---- Computing The Mass Rate in Ckoke ----
if( TG < (HC - .5) )
{
HLX = AT * HL / ACH;
MGFLOW(PT, ZAV, TAV, MGS);
MGAS = MGAS + MGS * DT;
// --- Resolving The Danamic Model of Plunger -----
FF = MODY(DNTL, DTI, DXTP, MUL);
FHU = FF * HL * DNTL * DXTP * DXTP / (2 * DTI);
FF = MODY(DNTL, DTI, DXTK, MUL);
FHD = FF * TK * DNTL * DXTK * DXTK / (2 * DTI);
// ----Gas Damping Factor -----
UG = VISCO(ZAV, PT, TAV);
DTG = M * PT / (ZAV * R * TAV);
FG = JAIN(DTG, DTI, DXTP, UG);
FGU = FF * (HC - TG) * DTG * DXTP * DXTP / (2 * DTI);
FG = JAIN(DNTL, DTI, DXTK, UG);
FGD = FF * TK * DTG * DXTK * DXTK / (2 * DTI);
PPG = (PHC - DNTL * G * TK - FHD - FGD) * exp(-.5 * FA * (TP - TK));
DPHC = (DNTL * G - FA * PHC) * DXCL;
DPHC = DPHC + DPCS * exp(FA * (HC - CL));
DPPG = -.5 * FA * PPG * (DXTP - DXTK);
DPPG = exp(FA * (TP - TK)) * (DPHC - DNTL * G * DXTK);
PPU = PT * exp(FA * (HC - TG)) + DNTL * HL * G + FHU + FGU;
PPD = (PHC - DNTL * G * TK - FHD - FGD) / exp(FA * (TP - TK));
// --- Resolving The DXTP by Runge-Kutta Method ----
RUNGE(VO, FD, FT, PPD, PPU, DT, DXTP, XTPR);
if( CL > 0 )
{
KD = 1;
DXTG = DXTP - MLTP / (AT * DNTL);
DXTK = DXTP - (MG / FB - (TP - TK) * DPPG) / PPG;
MLC = AT * DNTL * DXTK - ML;
DXCL = -MLC / (AC * DNTL); // dXcg/dt = ?
DPCG = PCG * DXCL / (HC - CL); // dPcg/dt = ?
}
else
{
KD = 2;
HLX = AT * HL / ACH;
CL = 0;
DXCL = 0;
DXTG = DXTP - MLTP / (AT * DNTL);
DXTK = ML / (AT * DNTL);
MGC = (TP - TK) * DPPG + PPG * (DXTP - DXTK);
MGC = FB * MGC + MGTP - MG;
DPCG = -R * ZAV * TAV * MGC / (M * AC * HC);
}
DPTG = PTG * DXTP + (MLTP - MGS) / FB; // dPtg/dt = ?
DPTG = DPTG / (HC - TG);
PPCG = PCG + DPCG * DT;
PPTG = PTG + DPTG * DT;
PCS[I] = PPCG / exp(.5 * FA * (HC - CL));
PTS[I] = PPTG / exp(.5 * FA * (HC - TG));
XTP[I] = XTP[I - 1] + DXTP * DT;
XTG[I] = XTG[I - 1] + DXTG * DT;
XCL[I] = XCL[I - 1] + DXCL * DT;
XTK[I] = XTK[I - 1] + DXTK * DT;
XCH[I] = XCH[I - 1] + DXCH * DT;
DPCS = (PCS[I] - PCS[I - 1]) / DT;
}
else
{
DP = .5 * DNTL * (ACH / AK + EC) * DXTP * DXTP;
PT = PTMIN + DP;
PPU = PT + DNTL * G * (HC - TP) - FHU;
PPD = (PHC - DNTL * G * TK - FHD) / exp(FA * (TP - TK));
PPG = (PHC - DNTL * G * TK - FHD) / exp(.5 * FA * (TP - TK));
DPPG1 = -.5 * FA * PPG * (DXTP - DXTK);
DPPG2 = exp(-.5 * FA * (TP - TK));
DPPG = DPPG1 + .5 * FA * (TP - TK) * DXTK * DPPG2;
DXTK = ML / (AT * DNTL);
if( TP < HC - .01 )
{
KD = 3;
// --- dXtp/dt Resultin by Runge-Kutta Method ---
RUNGE(VO, FD, FT, PPD, PPU, DT, DXTP, XTPR);
MGC = (TP - TK) * DPPG + PPG * (DXTP - DXTK);
MGC = FB * MGC - MG;
DPCG = -R * ZAV * TAV * MGC / (M * AC * HC);
DXCL = 0;
XCL[I] = 0;
DXTG = 0;
XTG[I] = HC;
XTP[I] = XTP[I - 1] + DXTP * DT; // Modified Factor (1.5)
XTK[I] = XTK[I - 1] + DXTK * DT;
XCH[I] = XCH[I - 1] + DXCH * DT;
PPCG = PCG + DPCG * DT;
PCS[I] = PPCG / exp(.5 * FA * HC);
PTS[I] = PT;
}
else
{
KD = 4;
return;
}
}
}// End For
ACG = int(100 * CL) / 100.0;
ATP = int(100 * TP) / 100.0;
ATG = int(100 * TG) / 100.0;
ATK = int(100 * TK) / 100.0;
BCH = int(100 * CH) / 100.0;
AV = int(1000 * DXTP) / 1000.0;
ASD = K - int(K / 15) * 15 + 1;
EPC = 6650000;
EPT = 1000000;
FP = .00005;
TUPC1 = 350 - (PCS[K - 1] - EPC) * FP;
TUPC2 = 350 - (PCS[K] - EPC) * FP;
TUPT1 = 450 - (PTS[K - 1] - EPT) * FP;
TUPT2 = 450 - (PTS[K] - EPT) * FP;
K = K + 1;
XTK[0] = XTKK[K - 1];
XTP[0] = XTPK[K - 1];
XTG[0] = XTGK[K - 1];
XCL[0] = XCLK[K - 1];
XCH[0] = XCHK[K - 1];
PCS[0] = PCSK[K - 1];
PTS[0] = PTSK[K - 1];
}// End While(1)
}
void CCal_Zsqj::PREWF(double &PHC, double &PWF, double &Z, double &T)
{
double ERRP = .001;
double X1, X2, EG, EL, DLAV, PWF1;
FLOW:
// ---- Computing The Flow Pressure -----
X1 = PS / 1000000;
X2 = PWF / 1000000;
// EG = MG / (1.293 * RG);
EG = 0;
EL = 1 - EG;
DLAV = EL * DNTL + EG * (M * PWF / (Z * R * T));
PWF1 = PHC + DLAV * G * (H - HC);
if( fabs(PWF - PWF1) / PWF > ERRP)
{
PWF = PWF1;
goto FLOW;
}
else
PWF = .5 * (PWF1 + PWF);
}
// Computing The Mass Rate in Ckoke ----
void CCal_Zsqj::MGFLOW(double &PT, double &ZAV, double &TAV, double &MGS)
{
double PR, A, QH;
PR = PTMIN / PT;
if( PR < .5266 ) PR = .5266;
A = pow(PR , 2 / CPV) - pow(PR , (CPV + 1) / CPV);
A = sqrt(CPV * A / (CPV - 1));
QH = 4066 * PT * DD * DD * A;
QH = QH / sqrt(RG * TAV * ZAV) / 86400;
MGS = 1.293 * RG * QH / .865;
}
void CCal_Zsqj::RUNGE(double &VO, double &FD, double &FT, double &PPD, double &PPU, double &DT, double &DXTP, double &XTPR)
{
double X1, Y1, X2, Y2, X3, Y3, X4, Y4;
double TP = 0;
X1 = VO;
Y1 = -FD * X1 + FT * (PPD - PPU) - G;
X2 = VO + .5 * Y1 * DT;
Y2 = -FD * X2 + FT * (PPD - PPU) - G;
X3 = VO + .5 * Y2 * DT;
Y3 = -FD * X3 + FT * (PPD - PPU) - G;
X4 = VO + Y3 * DT;
Y4 = -FD * X4 + FT * (PPD - PPU) - G;
XTPR = TP + (X1 + 2 * (X2 + X3) + X4) * DT / 6;
DXTP = VO + (Y1 + 2 * (Y2 + Y3) + Y4) * DT / 6;
}
void CCal_Zsqj::TDOWN(double &PCMAX, double &PCMIN, double &PTMAX, double &PTMIN, double &TDG, double &TDL, double &HLI)
{
double VDG = 2.702;
double VDL = .344;
TDG = (HC - HLI) / VDG;
TDL = HLI / VDL;
}
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -