Compare commits

..

No commits in common. "c88990815a729d2d15d3c14e4cd411c081843f17" and "da4481b9d156d64a2041ace163c140b1dac1da8a" have entirely different histories.

12 changed files with 589 additions and 359 deletions

File diff suppressed because one or more lines are too long

View File

@ -16,7 +16,7 @@
<TargetCommonOption>
<Device>STM32H750XBHx</Device>
<Vendor>STMicroelectronics</Vendor>
<PackID>Keil.STM32H7xx_DFP.2.6.0</PackID>
<PackID>Keil.STM32H7xx_DFP.4.0.0</PackID>
<PackURL>https://www.keil.com/pack/</PackURL>
<Cpu>IRAM(0x20000000,0x00020000) IRAM2(0x24000000,0x00080000) IROM(0x08000000,0x00020000) CPUTYPE("Cortex-M7") FPU3(DFPU) CLOCK(12000000) ELITTLE</Cpu>
<FlashUtilSpec></FlashUtilSpec>

View File

@ -10,19 +10,31 @@ Detail *Detail_Construct(void) {
if (!pDetail) {
return NULL;
}
memset(pDetail, 0, sizeof(Detail));
pDetail->iNCC = 0;
for (int i = 0; i < NUMBEROFCOMPONENTS; i++) {
pDetail->aiCID[i] = -1;
}
pDetail->dOldMixID = 0.0;
pDetail->dOldPb = 0.0;
pDetail->dOldTb = 0.0;
pDetail->dOldPf = 0.0;
pDetail->dOldTf = 0.0;
for (int i = 0; i < NUMBEROFCOMPONENTS; i++) {
pDetail->dXi[i] = 0.0;
}
if (Detail_table(pDetail) != 0) {
free(pDetail);
return NULL;
@ -36,23 +48,26 @@ void Detail_Destroy(Detail *pDetail) {
}
}
int Detail_compositionchange(Detail *pDetail, const NGParSTRUCT *ptNGPar) {
int Detail_compositionchange(Detail *pDetail, NGParSTRUCT *ptNGPar) {
double dMixID = 0.0;
for (int i = 0; i < NUMBEROFCOMPONENTS; i++)
int i;
for (i = 0; i < NUMBEROFCOMPONENTS; i++)
dMixID += ((i + 2) * ptNGPar->adMixture[i]);
if (dMixID != pDetail->dOldMixID) {
pDetail->dOldMixID = dMixID;
return 1;
return true;
} else {
return false;
}
return 0;
}
void Detail_Run(Detail *pDetail, NGParSTRUCT *ptNGPar) {
const int bCompChange = Detail_compositionchange(pDetail, ptNGPar);
int i;
bool bCompChange = Detail_compositionchange(pDetail, ptNGPar);
ptNGPar->bForceUpdate = ptNGPar->bForceUpdate || bCompChange;
if (ptNGPar->bForceUpdate) {
pDetail->iNCC = 0;
for (int i = 0; i < NUMBEROFCOMPONENTS; i++) {
for (i = 0; i < NUMBEROFCOMPONENTS; i++) {
if (ptNGPar->adMixture[i] > 0.0) {
pDetail->aiCID[pDetail->iNCC] = i;
pDetail->dXi[pDetail->iNCC] = ptNGPar->adMixture[i];
@ -62,9 +77,9 @@ void Detail_Run(Detail *pDetail, NGParSTRUCT *ptNGPar) {
Detail_paramdl(pDetail);
Detail_chardl(pDetail, ptNGPar);
}
if (fabs(ptNGPar->dPb - pDetail->dOldPb) > P_CHG_TOL ||
fabs(ptNGPar->dTb - pDetail->dOldTb) > T_CHG_TOL ||
ptNGPar->bForceUpdate) {
if ((fabs(ptNGPar->dPb - pDetail->dOldPb) > P_CHG_TOL) ||
(fabs(ptNGPar->dTb - pDetail->dOldTb) > T_CHG_TOL) ||
(ptNGPar->bForceUpdate)) {
pDetail->dP = ptNGPar->dPb * 1.0e-6;
pDetail->dT = ptNGPar->dTb;
Detail_temp(pDetail);
@ -76,13 +91,13 @@ void Detail_Run(Detail *pDetail, NGParSTRUCT *ptNGPar) {
ptNGPar->dRhob = pDetail->dRhoTP;
pDetail->dOldTb = ptNGPar->dTb;
pDetail->dOldPb = ptNGPar->dPb;
ptNGPar->bForceUpdate = 1;
ptNGPar->bForceUpdate = true;
}
pDetail->dP = ptNGPar->dPf * 1.0e-6;
pDetail->dT = ptNGPar->dTf;
if ((fabs(ptNGPar->dTf - pDetail->dOldTf) > T_CHG_TOL) || (ptNGPar->bForceUpdate)) {
Detail_temp(pDetail);
ptNGPar->bForceUpdate = 1;
ptNGPar->bForceUpdate = true;
}
if ((fabs(ptNGPar->dPf - pDetail->dOldPf) > P_CHG_TOL) || (ptNGPar->bForceUpdate)) {
Detail_ddetail(pDetail, ptNGPar);
@ -97,7 +112,7 @@ void Detail_Run(Detail *pDetail, NGParSTRUCT *ptNGPar) {
ptNGPar->dFpv = sqrt(ptNGPar->dZb / ptNGPar->dZf);
} else {
ptNGPar->lStatus = GENERAL_CALCULATION_FAILURE;
ptNGPar->bForceUpdate = 0;
ptNGPar->bForceUpdate = false;
}
}
@ -105,6 +120,7 @@ int Detail_table(Detail *pDetail) {
if (!pDetail) {
return -1;
}
double adAn[58];
double adUn[58];
adAn[0] = 0.153832600;
@ -166,6 +182,7 @@ int Detail_table(Detail *pDetail) {
adAn[56] = -0.001226752;
adAn[57] = 0.002850908;
adUn[0] = 0.0;
adUn[1] = 0.5;
adUn[2] = 1.0;
@ -224,8 +241,10 @@ int Detail_table(Detail *pDetail) {
adUn[55] = 0.0;
adUn[56] = 1.0;
adUn[57] = 0.0;
memcpy(pDetail->adAn, adAn, sizeof(adAn));
memcpy(pDetail->adUn, adUn, sizeof(adUn));
for (int j = 0; j < NUMBEROFCOMPONENTS; j++) {
for (int k = j; k < NUMBEROFCOMPONENTS; k++) {
pDetail->adTable6Eij[j][k] = 1.0;
@ -368,7 +387,7 @@ int Detail_table(Detail *pDetail) {
pDetail->adTable6Gij[2][3] = 0.370296;
pDetail->adTable6Gij[2][5] = 1.673090;
const double adTableHhvMol[4][NUMBEROFCOMPONENTS] = {
double adTableHhvMol[4][NUMBEROFCOMPONENTS] = {
{
892.97, 0, 0, 1564.34, 2224.01, 45.074, 562.94, 286.63, 282.8, 0, 2874.2, 2883.82, 3535.98, 3542.89,
4203.23, 4862.87, 5522.4, 6182.91, 6842.69, 0, 0
@ -387,7 +406,7 @@ int Detail_table(Detail *pDetail) {
}
};
memcpy(pDetail->adTableHhvMol, adTableHhvMol, sizeof(adTableHhvMol));
const double adTableLhvMol[4][NUMBEROFCOMPONENTS] = {
double adTableLhvMol[4][NUMBEROFCOMPONENTS] = {
{
802.82, 0, 0, 1429.12, 2043.71, 0, 517.87, 241.56, 282.8, 0, 2648.83, 2658.45, 3265.54, 3272.45, 3887.71,
4502.28, 5116.73, 5732.17, 6346.88, 0, 0
@ -410,6 +429,24 @@ int Detail_table(Detail *pDetail) {
}
void Detail_paramdl(Detail *pDetail) {
const double adTable5Mri[21] = {
16.0430, 28.0135, 44.0100, 30.0700, 44.0970, 18.0153, 34.0820, 2.0159, 28.0100, 31.9988, 58.1230, 58.1230,
72.1500, 72.1500, 86.1770, 100.2040, 114.2310, 128.2580, 142.2850, 4.0026, 39.9480
};
const double adTable5Ei[21] = {
151.318300, 99.737780, 241.960600, 244.166700, 298.118300, 514.015600, 296.355000, 26.957940, 105.534800,
122.766700, 324.068900, 337.638900, 365.599900, 370.682300, 402.636293, 427.722630, 450.325022, 470.840891,
489.558373, 2.610111, 119.629900
};
const double adTable5Ki[21] = {
0.4619255, 0.4479153, 0.4557489, 0.5279209, 0.5837490, 0.3825868, 0.4618263, 0.3514916, 0.4533894, 0.4186954,
0.6406937, 0.6341423, 0.6738577, 0.6798307, 0.7175118, 0.7525189, 0.7849550, 0.8152731, 0.8437826, 0.3589888,
0.4216551
};
const double adTable5Gi[21] = {
0.000000, 0.027815, 0.189065, 0.079300, 0.141239, 0.332500, 0.088500, 0.034369, 0.038953, 0.021000, 0.256692,
0.281835, 0.332267, 0.366911, 0.289731, 0.337542, 0.383381, 0.427354, 0.469659, 0.000000, 0.000000
};
for (int j = 0; j < NUMBEROFCOMPONENTS; j++) {
pDetail->adTable5Qi[j] = 0.0;
pDetail->adTable5Fi[j] = 0.0;
@ -424,28 +461,10 @@ void Detail_paramdl(Detail *pDetail) {
pDetail->adTable5Si[6] = 0.3900;
pDetail->adTable5Wi[5] = 1.0000;
for (int j = pDetail->iNCC - 1; j >= 0; j--) {
const double adTable5Ki[21] = {
0.4619255, 0.4479153, 0.4557489, 0.5279209, 0.5837490, 0.3825868, 0.4618263, 0.3514916, 0.4533894, 0.4186954,
0.6406937, 0.6341423, 0.6738577, 0.6798307, 0.7175118, 0.7525189, 0.7849550, 0.8152731, 0.8437826, 0.3589888,
0.4216551
};
const double adTable5Mri[21] = {
16.0430, 28.0135, 44.0100, 30.0700, 44.0970, 18.0153, 34.0820, 2.0159, 28.0100, 31.9988, 58.1230, 58.1230,
72.1500, 72.1500, 86.1770, 100.2040, 114.2310, 128.2580, 142.2850, 4.0026, 39.9480
};
pDetail->dMri[j] = adTable5Mri[pDetail->aiCID[j]];
pDetail->dKi[j] = adTable5Ki[pDetail->aiCID[j]];
}
for (int j = 0; j < pDetail->iNCC; j++) {
const double adTable5Gi[21] = {
0.000000, 0.027815, 0.189065, 0.079300, 0.141239, 0.332500, 0.088500, 0.034369, 0.038953, 0.021000, 0.256692,
0.281835, 0.332267, 0.366911, 0.289731, 0.337542, 0.383381, 0.427354, 0.469659, 0.000000, 0.000000
};
const double adTable5Ei[21] = {
151.318300, 99.737780, 241.960600, 244.166700, 298.118300, 514.015600, 296.355000, 26.957940, 105.534800,
122.766700, 324.068900, 337.638900, 365.599900, 370.682300, 402.636293, 427.722630, 450.325022, 470.840891,
489.558373, 2.610111, 119.629900
};
pDetail->dGi[j] = adTable5Gi[pDetail->aiCID[j]];
pDetail->dEi[j] = adTable5Ei[pDetail->aiCID[j]];
}
@ -488,8 +507,7 @@ void Detail_chardl(Detail *pDetail, NGParSTRUCT *ptNGPar)
ptNGPar->dLhvMol = 0.0;
for (int i = 0; i < pDetail->iNCC; i++) {
ptNGPar->dMrx += pDetail->dXi[i] * pDetail->dMri[i];
switch (ptNGPar->dCbtj)
{
switch (ptNGPar->dCbtj) {
case 2:
ptNGPar->dHhvMol += pDetail->adTableHhvMol[0][i] * ptNGPar->adMixture[i];
ptNGPar->dLhvMol += pDetail->adTableHhvMol[0][i] * ptNGPar->adMixture[i];
@ -502,7 +520,6 @@ void Detail_chardl(Detail *pDetail, NGParSTRUCT *ptNGPar)
ptNGPar->dHhvMol += pDetail->adTableHhvMol[2][i] * ptNGPar->adMixture[i];
ptNGPar->dLhvMol += pDetail->adTableLhvMol[2][i] * ptNGPar->adMixture[i];
break;
default: ;
}
k2p5 += pDetail->dXi[i] * pow(pDetail->dKi[i], 2.5);
u2p5 += pDetail->dXi[i] * pow(pDetail->dEi[i], 2.5);
@ -511,33 +528,33 @@ void Detail_chardl(Detail *pDetail, NGParSTRUCT *ptNGPar)
pDetail->dF += pDetail->dXi[i] * pDetail->dXi[i] * pDetail->dFi[i];
for (int j = i; j < pDetail->iNCC; j++) {
const double Xij = (i == j) ? pDetail->dXi[i] * pDetail->dXi[j] : 2.0 * pDetail->dXi[i] * pDetail->dXi[j];
double Xij = (i == j) ? pDetail->dXi[i] * pDetail->dXi[j] : 2.0 * pDetail->dXi[i] * pDetail->dXi[j];
if (pDetail->dKij[i][j] != 1.0) {
const double term = pow(pDetail->dKi[i] * pDetail->dKi[j], 2.5);
double term = pow(pDetail->dKi[i] * pDetail->dKi[j], 2.5);
k5p0 += Xij * (pow(pDetail->dKij[i][j], 5.0) - 1.0) * term;
}
if (pDetail->dUij[i][j] != 1.0) {
const double term = pow(pDetail->dEi[i] * pDetail->dEi[j], 2.5);
double term = pow(pDetail->dEi[i] * pDetail->dEi[j], 2.5);
u5p0 += Xij * (pow(pDetail->dUij[i][j], 5.0) - 1.0) * term;
}
if (pDetail->dGij[i][j] != 1.0) {
const double avgG = (pDetail->dGi[i] + pDetail->dGi[j]) / 2.0;
double avgG = (pDetail->dGi[i] + pDetail->dGi[j]) / 2.0;
pDetail->dW += Xij * (pDetail->dGij[i][j] - 1.0) * avgG;
}
const double Eij = pDetail->dEij[i][j] * sqrt(pDetail->dEi[i] * pDetail->dEi[j]);
const double Gij = pDetail->dGij[i][j] * (pDetail->dGi[i] + pDetail->dGi[j]) / 2.0;
const double e0p5 = sqrt(Eij);
const double e2p0 = Eij * Eij;
const double e3p0 = Eij * e2p0;
const double e3p5 = e3p0 * e0p5;
const double e4p5 = Eij * e3p5;
const double e6p0 = e3p0 * e3p0;
const double e7p5 = e4p5 * Eij * e2p0;
const double e9p5 = e7p5 * e2p0;
const double e11p0 = e4p5 * e4p5 * e2p0;
const double e12p0 = e11p0 * Eij;
const double e12p5 = e12p0 * e0p5;
const double s3 = Xij * pow(pow(pDetail->dKi[i], 3.0) * pow(pDetail->dKi[j], 3.0), 0.5);
double Eij = pDetail->dEij[i][j] * sqrt(pDetail->dEi[i] * pDetail->dEi[j]);
double Gij = pDetail->dGij[i][j] * (pDetail->dGi[i] + pDetail->dGi[j]) / 2.0;
double e0p5 = sqrt(Eij);
double e2p0 = Eij * Eij;
double e3p0 = Eij * e2p0;
double e3p5 = e3p0 * e0p5;
double e4p5 = Eij * e3p5;
double e6p0 = e3p0 * e3p0;
double e7p5 = e4p5 * Eij * e2p0;
double e9p5 = e7p5 * e2p0;
double e11p0 = e4p5 * e4p5 * e2p0;
double e12p0 = e11p0 * Eij;
double e12p5 = e12p0 * e0p5;
double s3 = Xij * pow(pow(pDetail->dKi[i], 3.0) * pow(pDetail->dKi[j], 3.0), 0.5);
pDetail->adBcoef[0] += s3;
pDetail->adBcoef[1] += s3 * e0p5;
pDetail->adBcoef[2] += s3 * Eij;
@ -559,8 +576,8 @@ void Detail_chardl(Detail *pDetail, NGParSTRUCT *ptNGPar)
}
}
ptNGPar->dHhvMol = ptNGPar->dHhvMol ;
ptNGPar->dLhvMol = ptNGPar->dLhvMol ;
ptNGPar->dHhvMol = ptNGPar->dHhvMol / ptNGPar->dMrx;
ptNGPar->dLhvMol = ptNGPar->dLhvMol / ptNGPar->dMrx;
for (int i = 0; i < 18; i++) {
pDetail->adBcoef[i] *= pDetail->adAn[i];
@ -572,20 +589,20 @@ void Detail_chardl(Detail *pDetail, NGParSTRUCT *ptNGPar)
void Detail_bvir(Detail *pDetail) {
pDetail->dB = pDetail->ddBdT = pDetail->dd2BdT2 = 0.0;
const double t = pDetail->dT;
const double t0p5 = sqrt(t);
const double t2p0 = t * t;
const double t3p0 = t * t2p0;
const double t3p5 = t3p0 * t0p5;
const double t4p5 = t * t3p5;
const double t6p0 = t3p0 * t3p0;
const double t11p0 = t4p5 * t4p5 * t2p0;
const double t7p5 = t6p0 * t * t0p5;
const double t9p5 = t7p5 * t2p0;
const double t12p0 = t9p5 * t0p5 * t2p0;
const double t12p5 = t12p0 * t0p5;
// double t1p5 = t * t0p5;
// double t4p0 = t2p0 * t2p0;
double t = pDetail->dT;
double t0p5 = sqrt(t);
double t2p0 = t * t;
double t3p0 = t * t2p0;
double t3p5 = t3p0 * t0p5;
double t4p5 = t * t3p5;
double t6p0 = t3p0 * t3p0;
double t11p0 = t4p5 * t4p5 * t2p0;
double t7p5 = t6p0 * t * t0p5;
double t9p5 = t7p5 * t2p0;
double t12p0 = t9p5 * t0p5 * t2p0;
double t12p5 = t12p0 * t0p5;
double t1p5 = t * t0p5;
double t4p0 = t2p0 * t2p0;
double Bx[18];
Bx[0] = pDetail->adBcoef[0];
Bx[1] = pDetail->adBcoef[1] / t0p5;
@ -633,22 +650,25 @@ void Detail_bvir(Detail *pDetail) {
void Detail_temp(Detail *pDetail) {
Detail_bvir(pDetail);
const double tr = pDetail->dT / pDetail->dU;
const double tr0p5 = sqrt(tr);
const double tr1p5 = tr * tr0p5;
const double tr2p0 = tr * tr;
const double tr3p0 = tr * tr2p0;
const double tr4p0 = tr * tr3p0;
const double tr5p0 = tr * tr4p0;
const double tr6p0 = tr * tr5p0;
const double tr7p0 = tr * tr6p0;
const double tr8p0 = tr * tr7p0;
const double tr9p0 = tr * tr8p0;
const double tr11p0 = tr6p0 * tr5p0;
const double tr13p0 = tr6p0 * tr7p0;
const double tr21p0 = tr9p0 * tr9p0 * tr3p0;
const double tr22p0 = tr * tr21p0;
const double tr23p0 = tr * tr22p0;
double tr = pDetail->dT / pDetail->dU;
double tr0p5, tr1p5, tr2p0, tr3p0, tr4p0, tr5p0, tr6p0;
double tr7p0, tr8p0, tr9p0, tr11p0, tr13p0, tr21p0;
double tr22p0, tr23p0 ;
tr0p5 = sqrt(tr);
tr1p5 = tr * tr0p5;
tr2p0 = tr * tr;
tr3p0 = tr * tr2p0;
tr4p0 = tr * tr3p0;
tr5p0 = tr * tr4p0;
tr6p0 = tr * tr5p0;
tr7p0 = tr * tr6p0;
tr8p0 = tr * tr7p0;
tr9p0 = tr * tr8p0;
tr11p0 = tr6p0 * tr5p0;
tr13p0 = tr6p0 * tr7p0;
tr21p0 = tr9p0 * tr9p0 * tr3p0;
tr22p0 = tr * tr21p0;
tr23p0 = tr * tr22p0;
pDetail->adFn[12] = pDetail->adAn[12] * pDetail->dF * tr6p0;
pDetail->adFn[13] = pDetail->adAn[13] / tr2p0;
@ -700,26 +720,30 @@ void Detail_temp(Detail *pDetail) {
}
void Detail_ddetail(Detail *pDetail, NGParSTRUCT *ptNGPar) {
double xnumer, xdenom;
const int imax = 150;
int imax, i;
double epsp, epsr, epsmin;
double x1, x2, x3, y1, y2, y3;
double delx, delprv, delmin, delbis, xnumer, xdenom, sgndel;
double y2my3, y3my1, y1my2, boundn;
imax = 150;
epsp = 1.e-6;
epsr = 1.e-6;
epsmin = 1.e-7;
pDetail->dRho = 0.0;
Detail_braket(pDetail, ptNGPar);
if (ptNGPar->lStatus == MAX_NUM_OF_ITERATIONS_EXCEEDED ||
ptNGPar->lStatus == NEGATIVE_DENSITY_DERIVATIVE) {
return;
}
double x1 = pDetail->dRhoL;
double x2 = pDetail->dRhoH;
double y1 = pDetail->dPRhoL - pDetail->dP;
double y2 = pDetail->dPRhoH - pDetail->dP;
double delx = x1 - x2;
double delprv = delx;
double x3 = x1;
double y3 = y1;
for (int i = 0; i < imax; i++) {
const double epsmin = 1.e-7;
const double epsr = 1.e-6;
const double epsp = 1.e-6;
x1 = pDetail->dRhoL;
x2 = pDetail->dRhoH;
y1 = pDetail->dPRhoL - pDetail->dP;
y2 = pDetail->dPRhoH - pDetail->dP;
delx = x1 - x2;
delprv = delx;
x3 = x1;
y3 = y1;
for (i = 0; i < imax; i++) {
if (y2 * y3 > 0.0) {
x3 = x1;
y3 = y1;
@ -734,16 +758,16 @@ void Detail_ddetail(Detail *pDetail, NGParSTRUCT *ptNGPar) {
y2 = y3;
y3 = y1;
}
const double delmin = epsmin * fabs(x2);
const double delbis = 0.5 * (x3 - x2);
delmin = epsmin * fabs(x2);
delbis = 0.5 * (x3 - x2);
if (fabs(delprv) < delmin || fabs(y1) < fabs(y2)) {
delx = delbis;
delprv = delbis;
} else {
if (x3 != x1) {
const double y2my3 = y2 - y3;
const double y3my1 = y3 - y1;
const double y1my2 = y1 - y2;
y2my3 = y2 - y3;
y3my1 = y3 - y1;
y1my2 = y1 - y2;
xdenom = -(y1my2) * (y2my3) * (y3my1);
xnumer = x1 * y2 * y3 * (y2my3)
+ x2 * y3 * y1 * (y3my1)
@ -765,11 +789,11 @@ void Detail_ddetail(Detail *pDetail, NGParSTRUCT *ptNGPar) {
return;
}
if (fabs(delx) < delmin) {
const double sgndel = delbis / fabs(delbis);
sgndel = delbis / fabs(delbis);
delx = 1.0000009 * sgndel * delmin;
delprv = delx;
}
const double boundn = delx * (x2 + delx - x3);
boundn = delx * (x2 + delx - x3);
if (boundn > 0.0) {
delx = delbis;
delprv = delbis;
@ -785,21 +809,23 @@ void Detail_ddetail(Detail *pDetail, NGParSTRUCT *ptNGPar) {
}
void Detail_braket(Detail *pDetail, NGParSTRUCT *ptNGPar) {
double rho2;
const int imax = 200;
double rho1 = 0.0;
double p1 = 0.0;
double rhomax = 1.0 / pDetail->dKp3;
int imax, it;
double del, rhomax, videal;
double rho1, rho2, p1, p2;
imax = 200;
rho1 = 0.0;
p1 = 0.0;
rhomax = 1.0 / pDetail->dKp3;
if (pDetail->dT > 1.2593 * pDetail->dU)
rhomax = 20.0 * rhomax;
const double videal = RGASKJ * pDetail->dT / pDetail->dP;
videal = RGASKJ * pDetail->dT / pDetail->dP;
if (fabs(pDetail->dB) < (0.167 * videal)) {
rho2 = 0.95 / (videal + pDetail->dB);
} else {
rho2 = 1.15 / videal;
}
double del = rho2 / 20.0;
for (int it = 0; it < imax; it++) {
del = rho2 / 20.0;
for (it = 0; it < imax; it++) {
if (rho2 > rhomax && ptNGPar->lStatus != MAX_DENSITY_IN_BRAKET_EXCEEDED) {
ptNGPar->lStatus = MAX_DENSITY_IN_BRAKET_EXCEEDED;
del = 0.01 * (rhomax - rho1) + (pDetail->dP / (RGASKJ * pDetail->dT)) / 20.0;
@ -807,7 +833,7 @@ void Detail_braket(Detail *pDetail, NGParSTRUCT *ptNGPar) {
continue;
}
Detail_pdetail(pDetail, rho2);
const double p2 = pDetail->dPCalc;
p2 = pDetail->dPCalc;
if (p2 > pDetail->dP) {
pDetail->dRhoL = rho1;
pDetail->dPRhoL = p1;
@ -830,6 +856,7 @@ void Detail_braket(Detail *pDetail, NGParSTRUCT *ptNGPar) {
}
ptNGPar->lStatus = MAX_NUM_OF_ITERATIONS_EXCEEDED;
pDetail->dRho = rho2;
return;
}
void Detail_pdetail(Detail *pDetail, double dD) {
@ -837,19 +864,20 @@ void Detail_pdetail(Detail *pDetail, double dD) {
}
double Detail_zdetail(Detail *pDetail, double d) {
const double D1 = pDetail->dKp3 * d;
const double D2 = D1 * D1;
const double D3 = D2 * D1;
const double D4 = D3 * D1;
const double D5 = D4 * D1;
const double D6 = D5 * D1;
const double D7 = D6 * D1;
const double D8 = D7 * D1;
const double D9 = D8 * D1;
const double exp1 = exp(-D1);
const double exp2 = exp(-D2);
const double exp3 = exp(-D3);
const double exp4 = exp(-D4);
double D1, D2, D3, D4, D5, D6, D7, D8, D9, exp1, exp2, exp3, exp4;
D1 = pDetail->dKp3 * d;
D2 = D1 * D1;
D3 = D2 * D1;
D4 = D3 * D1;
D5 = D4 * D1;
D6 = D5 * D1;
D7 = D6 * D1;
D8 = D7 * D1;
D9 = D8 * D1;
exp1 = exp(-D1);
exp2 = exp(-D2);
exp3 = exp(-D3);
exp4 = exp(-D4);
pDetail->dZ = 1.0 + pDetail->dB * d
+ pDetail->adFn[12] * D1 * (exp3 - 1.0 - 3.0 * D3 * exp3)
+ (pDetail->adFn[13] + pDetail->adFn[14] + pDetail->adFn[15]) * D1 * (exp2 - 1.0 - 2.0 * D2 * exp2)
@ -880,19 +908,22 @@ double Detail_zdetail(Detail *pDetail, double d) {
}
double Detail_dZdT(Detail *pDetail, double d) {
const double D1 = pDetail->dKp3 * d;
const double D2 = D1 * D1;
const double D3 = D2 * D1;
const double D4 = D3 * D1;
const double D5 = D4 * D1;
const double D6 = D5 * D1;
const double D7 = D6 * D1;
const double D8 = D7 * D1;
const double exp1 = exp(-D1);
const double exp2 = exp(-D2);
const double exp3 = exp(-D3);
const double exp4 = exp(-D4);
for (int i = 12; i < 58; i++) {
double tmp;
int i;
double D1, D2, D3, D4, D5, D6, D7, D8, exp1, exp2, exp3, exp4;
D1 = pDetail->dKp3 * d;
D2 = D1 * D1;
D3 = D2 * D1;
D4 = D3 * D1;
D5 = D4 * D1;
D6 = D5 * D1;
D7 = D6 * D1;
D8 = D7 * D1;
exp1 = exp(-D1);
exp2 = exp(-D2);
exp3 = exp(-D3);
exp4 = exp(-D4);
for (i = 12; i < 58; i++) {
if (pDetail->adUn[i] && pDetail->adFn[i]) {
pDetail->fx[i] = (pDetail->adFn[i] * pDetail->adUn[i] * D1) / pDetail->dT;
} else {
@ -902,7 +933,7 @@ double Detail_dZdT(Detail *pDetail, double d) {
pDetail->ddZdT = d * pDetail->ddBdT;
if (pDetail->dF)
pDetail->ddZdT += pDetail->fx[12] - (pDetail->fx[12] * (1.0 - 3.0 * D3) * exp3);
double tmp = (1.0 - 2.0 * D2) * exp2;
tmp = (1.0 - 2.0 * D2) * exp2;
pDetail->ddZdT += (pDetail->fx[13] - (pDetail->fx[13] * tmp));
pDetail->ddZdT += pDetail->fx[14] - (pDetail->fx[14] * tmp);
pDetail->ddZdT += pDetail->fx[15] - (pDetail->fx[15] * tmp);
@ -935,19 +966,22 @@ double Detail_dZdT(Detail *pDetail, double d) {
}
double Detail_d2ZdT2(Detail *pDetail, double d) {
const double D1 = pDetail->dKp3 * d;
const double D2 = D1 * D1;
const double D3 = D2 * D1;
const double D4 = D3 * D1;
const double D5 = D4 * D1;
const double D6 = D5 * D1;
const double D7 = D6 * D1;
const double D8 = D7 * D1;
const double exp1 = exp(-D1);
const double exp2 = exp(-D2);
const double exp3 = exp(-D3);
const double exp4 = exp(-D4);
for (int i = 12; i < 58; i++) {
double tmp;
int i;
double D1, D2, D3, D4, D5, D6, D7, D8, exp1, exp2, exp3, exp4;
D1 = pDetail->dKp3 * d;
D2 = D1 * D1;
D3 = D2 * D1;
D4 = D3 * D1;
D5 = D4 * D1;
D6 = D5 * D1;
D7 = D6 * D1;
D8 = D7 * D1;
exp1 = exp(-D1);
exp2 = exp(-D2);
exp3 = exp(-D3);
exp4 = exp(-D4);
for (i = 12; i < 58; i++) {
if (pDetail->adUn[i] && pDetail->adFn[i]) {
pDetail->fx[i] = (pDetail->adFn[i] * D1 * pDetail->adUn[i] * (pDetail->adUn[i] + 1.0)) / (
pDetail->dT * pDetail->dT);
@ -958,13 +992,13 @@ double Detail_d2ZdT2(Detail *pDetail, double d) {
pDetail->dd2ZdT2 = d * pDetail->dd2BdT2;
if (pDetail->dF)
pDetail->dd2ZdT2 += pDetail->fx[12] - (pDetail->fx[12] * (1.0 - 3.0 * D3) * exp3);
double tmp = (1.0 - 2.0 * D2) * exp2;
tmp = (1.0 - 2.0 * D2) * exp2;
pDetail->dd2ZdT2 += -pDetail->fx[13] + (pDetail->fx[13] * tmp);
pDetail->dd2ZdT2 += -pDetail->fx[14] + (pDetail->fx[14] * tmp);
pDetail->dd2ZdT2 += -pDetail->fx[15] + (pDetail->fx[15] * tmp);
tmp = (1.0 - 4.0 * D4) * exp4;
pDetail->dd2ZdT2 += -pDetail->fx[16] + (pDetail->fx[16] * tmp);
pDetail->dd2ZdT2 += -pDetail->fx[17] + pDetail->fx[17] * tmp;
pDetail->dd2ZdT2 += -pDetail->fx[17] + (pDetail->fx[17] * tmp);
pDetail->dd2ZdT2 = pDetail->dd2ZdT2 + (pDetail->fx[18] + pDetail->fx[19]) * D1 * 2.0
+ (pDetail->fx[21] + pDetail->fx[22]) * D1 * (2.0 - 2.0 * D2) * exp2
+ (pDetail->fx[23] + pDetail->fx[24] + pDetail->fx[25]) * D1 * (2.0 - 4.0 * D4) * exp4
@ -992,19 +1026,21 @@ double Detail_d2ZdT2(Detail *pDetail, double d) {
double Detail_dZdD(Detail *pDetail, double d) {
double temp, temp1, temp2, temp3;
const double D1 = pDetail->dKp3 * d;
const double D2 = D1 * D1;
const double D3 = D2 * D1;
const double D4 = D3 * D1;
const double D5 = D4 * D1;
const double D6 = D5 * D1;
const double D7 = D6 * D1;
const double D8 = D7 * D1;
const double exp1 = exp(-D1);
const double exp2 = exp(-D2);
const double exp3 = exp(-D3);
const double exp4 = exp(-D4);
for (int i = 12; i < 58; i++) {
int i;
double D1, D2, D3, D4, D5, D6, D7, D8, exp1, exp2, exp3, exp4;
D1 = pDetail->dKp3 * d;
D2 = D1 * D1;
D3 = D2 * D1;
D4 = D3 * D1;
D5 = D4 * D1;
D6 = D5 * D1;
D7 = D6 * D1;
D8 = D7 * D1;
exp1 = exp(-D1);
exp2 = exp(-D2);
exp3 = exp(-D3);
exp4 = exp(-D4);
for (i = 12; i < 58; i++) {
pDetail->fx[i] = pDetail->adFn[i];
}
pDetail->ddZdD = pDetail->dB / pDetail->dKp3;
@ -1120,10 +1156,11 @@ double Detail_dZdD(Detail *pDetail, double d) {
return pDetail->ddZdD;
}
void Detail_relativedensity(const Detail *pDetail, NGParSTRUCT *ptNGPar) {
void Detail_relativedensity(Detail *pDetail, NGParSTRUCT *ptNGPar) {
double dBX, dZa;
const double dMWair = 28.96256;
const double dBX = -0.12527 + 5.91e-4 * ptNGPar->dTb - 6.62e-7 * ptNGPar->dTb * ptNGPar->dTb;
const double dZa = 1.0 + (dBX * pDetail->dP) / (RGASKJ * ptNGPar->dTb);
dBX = -0.12527 + 5.91e-4 * ptNGPar->dTb - 6.62e-7 * ptNGPar->dTb * ptNGPar->dTb;
dZa = 1.0 + (dBX * pDetail->dP) / (RGASKJ * ptNGPar->dTb);
ptNGPar->dRD_Ideal = ptNGPar->dMrx / dMWair;
ptNGPar->dRD_Real = ptNGPar->dRD_Ideal * (dZa / ptNGPar->dZb);
}

View File

@ -6,6 +6,7 @@
#define _DETAIL_H
#include "NGCal.h"
#include <stdbool.h>
typedef struct Detail {
@ -87,7 +88,7 @@ Detail *Detail_Construct(void);
void Detail_Destroy(Detail *pDetail);
int Detail_compositionchange(Detail *pDetail, const NGParSTRUCT *pAGA10);
int Detail_compositionchange(Detail *pDetail, NGParSTRUCT *pAGA10);
int Detail_table(Detail *pDetail);
@ -105,7 +106,7 @@ void Detail_pdetail(Detail *pDetail, double dRho);
void Detail_ddetail(Detail *pDetail, NGParSTRUCT *pAGA10);
void Detail_relativedensity(const Detail *pDetail, NGParSTRUCT *pAGA10);
void Detail_relativedensity(Detail *pDetail, NGParSTRUCT *pAGA10);
double Detail_zdetail(Detail *pDetail, double dRho);

View File

@ -3,36 +3,6 @@
//
#include "NGCal.h"
#include "FlowCal.h"
#include "math.h"
double format_double(double value, int digits) {
// 处理默认位数4位
if (digits == 0) {
digits = 4;
}
// 验证位数有效性
if (digits < 1 || digits > 5) {
fprintf(stderr, "Error: Invalid digit value (must be 1-5 or 0 for default)\n");
return NAN; // 返回 NaN 表示错误
}
char format_str[10];
char buffer[50];
double result;
// 生成动态格式字符串(如 "%.4f"
snprintf(format_str, sizeof(format_str), "%%.%df", digits);
// 格式化数值到字符串
snprintf(buffer, sizeof(buffer), format_str, value);
// 转换回 double
result = strtod(buffer, NULL);
return result;
}
void OFlowCal(FlowParSTRUCT *ptFlowPar, NGParSTRUCT *ptNGPar) {
@ -80,7 +50,6 @@ void OFlowCal(FlowParSTRUCT *ptFlowPar, NGParSTRUCT *ptNGPar) {
ptFlowPar->dPb_M = (101325);
ptFlowPar->dTb_M = (293.15);
break;
default: ;
}
double ngArray[NUMBEROFCOMPONENTS];
@ -92,24 +61,24 @@ void OFlowCal(FlowParSTRUCT *ptFlowPar, NGParSTRUCT *ptNGPar) {
Crit(ptNGPar, 0);
ptFlowPar->dFpv =format_double( ptNGPar->dFpv,4);
ptFlowPar->dFpv = ptNGPar->dFpv;
ptFlowPar->dOrificeD =format_double( ptFlowPar->dOrificeD * (
1 + 0.000001 * CaiLiaoPzxs(ptFlowPar->dOrificeMaterial) * (ptFlowPar->dTf - 293.15)),2);
ptFlowPar->dPipeD =format_double( ptFlowPar->dPipeD * (1 + 0.000001 * CaiLiaoPzxs(ptFlowPar->dPipeMaterial) * (
ptFlowPar->dTf - 293.15)),2);
ptFlowPar->dBeta =format_double( ptFlowPar->dOrificeD / ptFlowPar->dPipeD,4);
ptFlowPar->dOrificeD = ptFlowPar->dOrificeD * (
1 + 0.000001 * CaiLiaoPzxs(ptFlowPar->dOrificeMaterial) * (ptFlowPar->dTf - 293.15));
ptFlowPar->dPipeD = ptFlowPar->dPipeD * (1 + 0.000001 * CaiLiaoPzxs(ptFlowPar->dPipeMaterial) * (
ptFlowPar->dTf - 293.15));
ptFlowPar->dBeta = ptFlowPar->dOrificeD / ptFlowPar->dPipeD;
ptFlowPar->dE =format_double( calculateE(ptFlowPar->dBeta),4);
ptFlowPar->dFG =format_double( calculateFG(ptNGPar->dRD_Real),4);
ptFlowPar->dFT =format_double( calculateFT(ptFlowPar->dTb_M, ptFlowPar->dTf),4);
ptFlowPar->dE = calculateE(ptFlowPar->dBeta);
ptFlowPar->dFG = calculateFG(ptNGPar->dRD_Real);
ptFlowPar->dFT = calculateFT(ptFlowPar->dTb_M, ptFlowPar->dTf);
ptFlowPar->dKappa =format_double( calculateKappa(ptNGPar->dZf),4);
ptFlowPar->dDViscosity =format_double( Dlndjs(ptFlowPar->dPf / 1e6, ptFlowPar->dTf),5);
ptFlowPar->dDExpCoefficient =format_double( calculateEpsilon(ptFlowPar->dPf, ptFlowPar->dDp,
ptFlowPar->dBeta, ptFlowPar->dKappa),4);
ptFlowPar->dKappa = calculateKappa(ptNGPar->dZf);
ptFlowPar->dDViscosity = Dlndjs(ptFlowPar->dPf / 1e6, ptFlowPar->dTf);
ptFlowPar->dDExpCoefficient = calculateEpsilon(ptFlowPar->dPf, ptFlowPar->dDp,
ptFlowPar->dBeta, ptFlowPar->dKappa);
double D = ptFlowPar->dPipeD / 1000.0;
double d = ptFlowPar->dOrificeD / 1000.0;
@ -123,25 +92,38 @@ void OFlowCal(FlowParSTRUCT *ptFlowPar, NGParSTRUCT *ptNGPar) {
double Qf_initial = (C_initial * ptFlowPar->dE * ptFlowPar->dDExpCoefficient * M_PI * pow(d, 2) / 4)
* sqrt(2 * deltaP / (ptNGPar->dRhof * (1 - pow(beta, 4))));
ptFlowPar->dVFlowf = Qf_initial;
double tolerance = 1e-6;
int maxIter = 100;
double currentC = C_initial;
double currentReD = calculateReD(Qf_initial, D, ptNGPar->dRhof, ptFlowPar->dDViscosity);
int iter = 0;
double prevC = 0;
do {
int maxIter = 100;
prevC = currentC;
currentC = calculateCd(beta, currentReD, ptFlowPar->dPipeD, ptFlowPar->dPtmode);
double Qf = (currentC * ptFlowPar->dDExpCoefficient * M_PI * pow(d, 2) / 4)
* sqrt(2 * deltaP / (ptNGPar->dRhof * (1 - pow(beta, 4))));
ptFlowPar->dVFlowf = Qf;
currentReD = calculateReD(Qf, D, ptNGPar->dRhof, ptFlowPar->dDViscosity);
iter++;
if (iter > maxIter) {
fprintf(stderr, "\n");
fprintf(stderr, "<EFBFBD><EFBFBD><EFBFBD><EFBFBD>δ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>\n");
}
} while (fabs(currentC - prevC) / currentC > tolerance);
double K = calculateK(ptFlowPar->dPipeType);
double G_me = calculateRoughnessFactor(ptFlowPar->dPipeD, K, currentC);
double C_corrected = currentC * G_me;
@ -149,11 +131,17 @@ void OFlowCal(FlowParSTRUCT *ptFlowPar, NGParSTRUCT *ptNGPar) {
ptFlowPar->dCd = C_corrected;
ptFlowPar->dRoughNessPipe = G_me;
ptFlowPar->dRnPipe = currentReD;
double Qn = ptFlowPar->dVFlowf * (ptFlowPar->dFpv * ptFlowPar->dFpv * P1 / ptFlowPar->dPb_M)
* (ptFlowPar->dTb_M) / Tf;
ptFlowPar->dVFlowb = Qn;
ptFlowPar->dMFlowb = ptFlowPar->dVFlowb * ptNGPar->dRhob;
ptFlowPar->dEFlowb = ptFlowPar->dVFlowb * ptNGPar->dHhvMol*ptFlowPar->dPb_M*1e-6/RGASKJ/ptFlowPar->dTb_M;
ptFlowPar->dEFlowb = ptFlowPar->dVFlowb * ptNGPar->dHhvMol;
ptFlowPar->dVelocityFlow = ptFlowPar->dVFlowf / (M_PI * pow((ptFlowPar->dPipeD / 2000), 2));
}
@ -162,7 +150,7 @@ void OFlowCal(FlowParSTRUCT *ptFlowPar, NGParSTRUCT *ptNGPar) {
double CaiLiaoPzxs(const int tempCaiLiao) {
double CaiLiaoPzxs(int tempCaiLiao) {
double CaiLiaoPzxs = 0;
switch (tempCaiLiao) {
case 0:
@ -224,7 +212,6 @@ double CaiLiaoPzxs(const int tempCaiLiao) {
case 14:
CaiLiaoPzxs = 17.2;
break;
default: ;
}
return CaiLiaoPzxs;
}
@ -331,7 +318,7 @@ double calculateKappa(double dZf) {
double calculateReD(double Qf, double D, double rho, double mu) {
return (4 * Qf * rho*1000) / (M_PI * D * mu);
return (4 * Qf * rho) / (M_PI * D * mu);
}
@ -392,7 +379,7 @@ double Dlndjs(double tempP_jy, double tempT) {
double Dlndjs_Dlnd_P[11] = {0.1, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10};
double ky, kx;
double s1, s2, ky, kx;
int i, m = 0, n = 0;
if (tempT < Dlndjs_Dlnd_T[0]) {
@ -434,7 +421,7 @@ double Dlndjs(double tempP_jy, double tempT) {
kx = 0;
}
double s1 = Dlndjs_Dlnd_Data[m][n] + (Dlndjs_Dlnd_Data[m][n + 1] - Dlndjs_Dlnd_Data[m][n]) * ky;
double s2 = Dlndjs_Dlnd_Data[m + 1][n] + (Dlndjs_Dlnd_Data[m + 1][n + 1] - Dlndjs_Dlnd_Data[m + 1][n]) * ky;
s1 = Dlndjs_Dlnd_Data[m][n] + (Dlndjs_Dlnd_Data[m][n + 1] - Dlndjs_Dlnd_Data[m][n]) * ky;
s2 = Dlndjs_Dlnd_Data[m + 1][n] + (Dlndjs_Dlnd_Data[m + 1][n + 1] - Dlndjs_Dlnd_Data[m + 1][n]) * ky;
return (s1 + (s2 - s1) * kx) / 100000.0;
}

View File

@ -1,7 +1,7 @@
#include "NGCal.h"
#include "Therm.h"
#include "Detail.h"
#include "math.h"
#include "String.h"
static Therm *ptTherm;
static Detail *ptDetail;
@ -39,6 +39,10 @@ double SOS(NGParSTRUCT *ptNGPar) {
}
double Crit(NGParSTRUCT *ptNGPar, double dPlenumVelocity) {
double DH, DDH, S, H;
double tolerance = 1.0;
double R, P, T, Z;
int i;
if (NULL == ptDetail || NULL == ptTherm) {
NGCal_UnInit();
if (NGCal_NGCal != NGCal_Init(ptNGPar)) {
@ -50,19 +54,18 @@ double Crit(NGParSTRUCT *ptNGPar, double dPlenumVelocity) {
Therm_Run(ptTherm, ptNGPar, ptDetail);
double DH = (ptNGPar->dSOS * ptNGPar->dSOS - dPlenumVelocity * dPlenumVelocity) / 2.0;
double S = ptNGPar->dS;
double H = ptNGPar->dH;
double R = ptNGPar->dRhof;
double P = ptNGPar->dPf;
double Z = ptNGPar->dZf;
double T = ptNGPar->dTf;
// DDH = 10.0;
for (int i = 1; i < MAX_NUM_OF_ITERATIONS; i++) {
double tolerance = 1.0;
Therm_HS_Mode(ptTherm, ptNGPar, ptDetail, H - DH, S, 1);
DH = (ptNGPar->dSOS * ptNGPar->dSOS - dPlenumVelocity * dPlenumVelocity) / 2.0;
S = ptNGPar->dS;
H = ptNGPar->dH;
R = ptNGPar->dRhof;
P = ptNGPar->dPf;
Z = ptNGPar->dZf;
T = ptNGPar->dTf;
DDH = 10.0;
for (i = 1; i < MAX_NUM_OF_ITERATIONS; i++) {
Therm_HS_Mode(ptTherm, ptNGPar, ptDetail, H - DH, S, true);
Therm_Run(ptTherm, ptNGPar, ptDetail);
double DDH = DH;
DDH = DH;
DH = (ptNGPar->dSOS * ptNGPar->dSOS - dPlenumVelocity * dPlenumVelocity) / 2.0;
if (fabs(DDH - DH) < tolerance) break;
}
@ -73,13 +76,14 @@ double Crit(NGParSTRUCT *ptNGPar, double dPlenumVelocity) {
return ptNGPar->dCstar;
}
double Cperf(const NGParSTRUCT *ptNGPar) {
double k = ptNGPar->dKappa;
double root = 2.0 / (k + 1.0);
double exponent = (k + 1.0) / (k - 1.0);
double Cperf(NGParSTRUCT *ptNGPar) {
double k, root, exponent;
k = ptNGPar->dKappa;
root = 2.0 / (k + 1.0);
exponent = (k + 1.0) / (k - 1.0);
return (sqrt(k * pow(root, exponent)));
}
double CRi(const NGParSTRUCT *ptNGPar) {
double CRi(NGParSTRUCT *ptNGPar) {
return (Cperf(ptNGPar) / sqrt(ptNGPar->dZf));
}

View File

@ -1,10 +1,20 @@
/*************************************************************************
* <EFBFBD>ļ<EFBFBD>: NGCal.h
**************************************************************************/
#ifndef _NGCal_H
#define _NGCal_H
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#define NORMAL 9000
#define NGCal_NGCal 9001
#define MEMORY_ALLOCATION_ERROR 9002
@ -14,49 +24,71 @@
#define MAX_DENSITY_IN_BRAKET_EXCEEDED 9006
#define FLOW_CALC_ERROR 9007
#define FLOW_CALC_DIEDAI_ERROR 9008
#define NUMBEROFCOMPONENTS 21
#define M_PI 3.1415926535897932
#define MAX_NUM_OF_ITERATIONS 100
#define P_CHG_TOL 0.001
#define T_CHG_TOL 0.001
#define P_MAX 1.379e8
#define P_MIN 0.0
#define T_MAX 473.15
#define T_MIN 143.0
#define RGASKJ 8.314510e-3
#define RGAS 8.314510
typedef struct tagNGParSTRUCT
{
long lStatus;
int bForceUpdate;
double adMixture[21];
int dCbtj;
double dPb;
double dTb;
double dPf;
double dTf;
double dMrx;
double dZb;
double dZf;
double dFpv;
double dDb;
double dDf;
double dRhob;
double dRhof;
double dRD_Ideal;
double dRD_Real;
double dHo;
double dH;
double dS;
double dCpi;
double dCp;
double dCv;
double dk;
double dKappa;
double dSOS;
double dCstar;
double dHhvMol;
double dLhvMol;
long lStatus;
int bForceUpdate;
double adMixture[21];
int dCbtj;
double dPb;
double dTb;
double dPf;
double dTf;
double dMrx;
double dZb;
double dZf;
double dFpv;
double dDb;
double dDf;
double dRhob;
double dRhof;
double dRD_Ideal;
double dRD_Real;
double dHo;
double dH;
double dS;
double dCpi;
double dCp;
double dCv;
double dk;
double dKappa;
double dSOS;
double dCstar;
double dHhvMol;
double dLhvMol;
} NGParSTRUCT;

View File

@ -259,7 +259,7 @@ void Therm_CprCvrHS(Therm *therm, NGParSTRUCT *ptNGPar, Detail *detail) {
}
void Therm_HS_Mode(Therm *therm, NGParSTRUCT *ptNGPar, Detail *detail, double H, double S, int bGuess) {
void Therm_HS_Mode(Therm *therm, NGParSTRUCT *ptNGPar, Detail *detail, double H, double S, bool bGuess) {
double s0 = S;
double h0 = H;
double t1, t2, tmin, tmax;

View File

@ -51,7 +51,7 @@ double Therm_So(Therm *therm, NGParSTRUCT *ptNGPar);
void Therm_CprCvrHS(Therm *therm, NGParSTRUCT *ptNGPar, Detail *detail);
void Therm_HS_Mode(Therm *therm, NGParSTRUCT *ptNGPar, Detail *detail, double H, double S, int bGuess);
void Therm_HS_Mode(Therm *therm, NGParSTRUCT *ptNGPar, Detail *detail, double H, double S, bool bGuess);
double Therm_H(Therm *therm, NGParSTRUCT *ptNGPar, Detail *detail);

View File

@ -9,11 +9,11 @@ int main() {
// 设置基本参数
flowParams.dPatm = 0.0981; // 标准大气压(bar)
flowParams.dPf = 1.48; // 压力(MPa)
flowParams.dPfType = 0; // 0=表压1=绝压
flowParams.dPf = 4; // 压力(MPa)
flowParams.dPfType = 1; // 0=表压1=绝压
flowParams.dDp = 12.50; // 差压(kPa)
flowParams.dTf = 15; // 温度(°C)
flowParams.dCbtj = 1; // 参比条件类型(0=标准状态)
flowParams.dTf = 10; // 温度(°C)
flowParams.dCbtj = 0; // 参比条件类型(0=标准状态)
// 设置管道参数
flowParams.dPipeD = 259.38; // 管道内径(mm)
@ -23,7 +23,7 @@ int main() {
// 设置材料参数
flowParams.dPipeMaterial = 2; // 20号钢
flowParams.dOrificeMaterial = 9; // 镍铬合金
flowParams.dOrificeMaterial = 10; // 镍铬合金
// 设置天然气组分(示例: 95%甲烷5%其他)
// 初始化天然气组分数组(GB/T 21446-2008 典型示例组成)
@ -31,21 +31,22 @@ int main() {
flowParams.dNG_Compents[i] = 0.0; // 先全部初始化为0
}
flowParams.dNG_Compents[0] = 92.47; // 甲烷(CH4)
flowParams.dNG_Compents[1] = 0.68; // 氮气(N2)
flowParams.dNG_Compents[2] = 1.75; // 二氧化碳(CO2)
flowParams.dNG_Compents[3] =3.5; // 乙烷(C2H6)
flowParams.dNG_Compents[4] = 0.98; // 丙烷(C3H8)
// 按照GB/T 21446-2008标准中典型天然气组分赋值(体积百分比)
flowParams.dNG_Compents[0] = 90.6724; // 甲烷(CH4)
flowParams.dNG_Compents[1] = 3.1284; // 氮气(N2)
flowParams.dNG_Compents[2] = 0.4676; // 二氧化碳(CO2)
flowParams.dNG_Compents[3] =4.5279; // 乙烷(C2H6)
flowParams.dNG_Compents[4] = 0.8280; // 丙烷(C3H8)
flowParams.dNG_Compents[5] = 0.00; // 水(H2O)
flowParams.dNG_Compents[6] = 0.00; // 硫化氢(H2S)
flowParams.dNG_Compents[7] = 0.0; // 氢气(H2)
flowParams.dNG_Compents[8] = 0.00; // 一氧化碳(CO)
flowParams.dNG_Compents[9] = 0.00; // 氧气(O2)
flowParams.dNG_Compents[10] = 0.34; // 异丁烷(i-C4H10)
flowParams.dNG_Compents[11] = 0.22; // 正丁烷(n-C4H10)
flowParams.dNG_Compents[12] = 0.0; // 异戊烷(i-C5H12)
flowParams.dNG_Compents[13] = 0.06; // 正戊烷(n-C5H12)
flowParams.dNG_Compents[14] = 0.0; // 己烷(C6H14)
flowParams.dNG_Compents[10] = 0.1037; // 异丁烷(i-C4H10)
flowParams.dNG_Compents[11] = 0.1563; // 正丁烷(n-C4H10)
flowParams.dNG_Compents[12] = 0.0321; // 异戊烷(i-C5H12)
flowParams.dNG_Compents[13] = 0.0443; // 正戊烷(n-C5H12)
flowParams.dNG_Compents[14] = 0.0393; // 己烷(C6H14)
flowParams.dNG_Compents[15] = 0.0; // 庚烷(C7H16)
flowParams.dNG_Compents[16] = 0.0; // 辛烷(C8H18)
flowParams.dNG_Compents[17] = 0.0; // 壬烷(C9H20)
@ -53,55 +54,6 @@ int main() {
flowParams.dNG_Compents[19] = 0.0; // 氦气(He)
flowParams.dNG_Compents[20] = 0.0; // 其他组分
// flowParams.dNG_Compents[0] = 88.36; // 甲烷(CH4)
// flowParams.dNG_Compents[1] = 0.68; // 氮气(N2)
// flowParams.dNG_Compents[2] = 1.57; // 二氧化碳(CO2)
// flowParams.dNG_Compents[3] =6.25; // 乙烷(C2H6)
// flowParams.dNG_Compents[4] = 2.4; // 丙烷(C3H8)
// flowParams.dNG_Compents[5] = 0.00; // 水(H2O)
// flowParams.dNG_Compents[6] = 0.00; // 硫化氢(H2S)
// flowParams.dNG_Compents[7] = 0.04; // 氢气(H2)
// flowParams.dNG_Compents[8] = 0.00; // 一氧化碳(CO)
// flowParams.dNG_Compents[9] = 0.00; // 氧气(O2)
// flowParams.dNG_Compents[10] = 0.15; // 异丁烷(i-C4H10)
// flowParams.dNG_Compents[11] = 0.35; // 正丁烷(n-C4H10)
// flowParams.dNG_Compents[12] = 0.05; // 异戊烷(i-C5H12)
// flowParams.dNG_Compents[13] = 0.1; // 正戊烷(n-C5H12)
// flowParams.dNG_Compents[14] = 0.01; // 己烷(C6H14)
// flowParams.dNG_Compents[15] = 0.0; // 庚烷(C7H16)
// flowParams.dNG_Compents[16] = 0.0; // 辛烷(C8H18)
// flowParams.dNG_Compents[17] = 0.0; // 壬烷(C9H20)
// flowParams.dNG_Compents[18] = 0.0; // 癸烷(C10H22)
// flowParams.dNG_Compents[19] = 0.04; // 氦气(He)
// flowParams.dNG_Compents[20] = 0.0; // 其他组分
// 按照GB/T 21446-2008标准中典型天然气组分赋值(体积百分比)
// flowParams.dNG_Compents[0] = 90.6724; // 甲烷(CH4)
// flowParams.dNG_Compents[1] = 3.1284; // 氮气(N2)
// flowParams.dNG_Compents[2] = 0.4676; // 二氧化碳(CO2)
// flowParams.dNG_Compents[3] =4.5279; // 乙烷(C2H6)
// flowParams.dNG_Compents[4] = 0.8280; // 丙烷(C3H8)
// flowParams.dNG_Compents[5] = 0.00; // 水(H2O)
// flowParams.dNG_Compents[6] = 0.00; // 硫化氢(H2S)
// flowParams.dNG_Compents[7] = 0.0; // 氢气(H2)
// flowParams.dNG_Compents[8] = 0.00; // 一氧化碳(CO)
// flowParams.dNG_Compents[9] = 0.00; // 氧气(O2)
// flowParams.dNG_Compents[10] = 0.1037; // 异丁烷(i-C4H10)
// flowParams.dNG_Compents[11] = 0.1563; // 正丁烷(n-C4H10)
// flowParams.dNG_Compents[12] = 0.0321; // 异戊烷(i-C5H12)
// flowParams.dNG_Compents[13] = 0.0443; // 正戊烷(n-C5H12)
// flowParams.dNG_Compents[14] = 0.0393; // 己烷(C6H14)
// flowParams.dNG_Compents[15] = 0.0; // 庚烷(C7H16)
// flowParams.dNG_Compents[16] = 0.0; // 辛烷(C8H18)
// flowParams.dNG_Compents[17] = 0.0; // 壬烷(C9H20)
// flowParams.dNG_Compents[18] = 0.0; // 癸烷(C10H22)
// flowParams.dNG_Compents[19] = 0.0; // 氦气(He)
// flowParams.dNG_Compents[20] = 0.0; // 其他组分
// flowParams.dNG_Compents[0] =96.5; // 甲烷(CH4)
// flowParams.dNG_Compents[1] =0.30; // 氮气(N2)
// flowParams.dNG_Compents[2] =0.6; // 二氧化碳(CO2)

217
User/main.c Normal file
View File

@ -0,0 +1,217 @@
/**
*********************************************************************
* @file main.c
* @author liaodeyun
* @version V1.0
* @date 2025-7-6
* @brief RT-Thread 3.0 + STM32
*********************************************************************
* @attention
*
* : H750 PRO STM32
* :http://www.firebbs.cn
* :https://fire-stm32.taobao.com
*
**********************************************************************
*/
/*
*************************************************************************
*
*************************************************************************
*/
#include "board.h"
#include "rtthread.h"
#include <stdio.h>
#include "NGCal.h"
#include "Therm.h"
#include "Detail.h"
#include "FlowCal.h"
/*
*************************************************************************
*
*************************************************************************
*/
/* 定义线程控制块 */
static rt_thread_t OFlowCal_thread = RT_NULL;
static rt_thread_t key_thread = RT_NULL;
/*
*************************************************************************
*
*************************************************************************
*/
static void OFlowCal_thread_entry(void* parameter);
static void key_thread_entry(void* parameter);
/*
*************************************************************************
* main
*************************************************************************
*/
/**
* @brief
* @param
* @retval
*/
int main(void)
{
/*
* RTT系统初始化已经在main函数之前完成
* component.c文件中的rtthread_startup()
* main函数中线线
*/
rt_kprintf("这是一个[野火]-STM32H750天然气流量计算机实验程序\n\n");
rt_kprintf("按下K1挂起线程按下K2恢复线程\n");
OFlowCal_thread = /* 线程控制块指针 */
rt_thread_create( "OFlowCal", /* 线程名字 */
OFlowCal_thread_entry, /* 线程入口函数 */
RT_NULL, /* 线程入口函数参数 */
512, /* 线程栈大小 */
3, /* 线程的优先级 */
20); /* 线程时间片 */
/* 启动线程,开启调度 */
if (OFlowCal_thread != RT_NULL)
rt_thread_startup(OFlowCal_thread);
else
return -1;
key_thread = /* 线程控制块指针 */
rt_thread_create( "key", /* 线程名字 */
key_thread_entry, /* 线程入口函数 */
RT_NULL, /* 线程入口函数参数 */
512, /* 线程栈大小 */
2, /* 线程的优先级 */
20); /* 线程时间片 */
/* 启动线程,开启调度 */
if (key_thread != RT_NULL)
rt_thread_startup(key_thread);
else
return -1;
}
/*
*************************************************************************
* 线
*************************************************************************
*/
static void OFlowCal_thread_entry(void* parameter)
{
while (1)
{
// 初始化结构体
FlowParSTRUCT flowParams = {0};
NGParSTRUCT ngParams = {0};
// 设置基本参数
flowParams.dPatm = 0.0981; // 标准大气压(bar)
flowParams.dPf = 1.48; // 压力(MPa)
flowParams.dPfType = 0; // 0=表压1=绝压
flowParams.dDp = 12.50; // 差压(kPa)
flowParams.dTf = 15.0; // 温度(°C)
flowParams.dCbtj = 0; // 参比条件类型(0=标准状态)
// 设置管道参数
flowParams.dPipeD = 259.38; // 管道内径(mm)
flowParams.dOrificeD = 150.25; // 孔板孔径(mm)
flowParams.dPipeType = 0; // 管道类型
flowParams.dPtmode = 0; // 取压方式(0=法兰取压1=角接取压)
// 设置材料参数
flowParams.dPipeMaterial = 2; // 20号钢
flowParams.dOrificeMaterial = 10; // 镍铬合金
// 设置天然气组分(示例: 95%甲烷5%其他)
// 初始化天然气组分数组(GB/T 21446-2008 典型示例组成)
for(int i=0; i<NUMBEROFCOMPONENTS; i++) {
flowParams.dNG_Compents[i] = 0.0; // 先全部初始化为0
}
// 按照GB/T 21446-2008标准中典型天然气组分赋值(体积百分比)
flowParams.dNG_Compents[0] = 88.36; // 甲烷(CH4)
flowParams.dNG_Compents[1] = 0.68; // 氮气(N2)
flowParams.dNG_Compents[2] = 1.57; // 二氧化碳(CO2)
flowParams.dNG_Compents[3] = 6.25; // 乙烷(C2H6)
flowParams.dNG_Compents[4] = 2.4; // 丙烷(C3H8)
flowParams.dNG_Compents[5] = 0.00; // 水(H2O)
flowParams.dNG_Compents[6] = 0.00; // 硫化氢(H2S)
flowParams.dNG_Compents[7] = 0.04; // 氢气(H2)
flowParams.dNG_Compents[8] = 0.00; // 一氧化碳(CO)
flowParams.dNG_Compents[9] = 0.00; // 氧气(O2)
flowParams.dNG_Compents[10] = 0.15; // 异丁烷(i-C4H10)
flowParams.dNG_Compents[11] = 0.35; // 正丁烷(n-C4H10)
flowParams.dNG_Compents[12] = 0.05; // 异戊烷(i-C5H12)
flowParams.dNG_Compents[13] = 0.1; // 正戊烷(n-C5H12)
flowParams.dNG_Compents[14] = 0.01; // 己烷(C6H14)
flowParams.dNG_Compents[15] = 0.0; // 庚烷(C7H16)
flowParams.dNG_Compents[16] = 0.0; // 辛烷(C8H18)
flowParams.dNG_Compents[17] = 0.0; // 壬烷(C9H20)
flowParams.dNG_Compents[18] = 0.0; // 癸烷(C10H22)
flowParams.dNG_Compents[19] = 0.04; // 氦气(He)
flowParams.dNG_Compents[20] = 0.0; // 其他组分
// 调用流量计算函数
OFlowCal(&flowParams, &ngParams);
// 打印计算结果
printf("========== 流量计算结果 ==========\n");
printf("工况流量: %.4f m3/s\n", flowParams.dVFlowf);
printf("标况流量: %.4f m3/s\n", flowParams.dVFlowb);
printf("标况质量流量: %.4f kg/s\n", flowParams.dMFlowb);
printf("标况能量流量: %.4f MJ/s\n", flowParams.dEFlowb / 1e6);
printf("管道内流速: %.2f m/s\n", flowParams.dVelocityFlow);
printf("流出系数: %.4f\n", flowParams.dCd);
printf("雷诺数: %.2f\n", flowParams.dRnPipe);
}
}
static void key_thread_entry(void* parameter)
{
rt_err_t uwRet = RT_EOK;
while (1)
{
if( Key_Scan(KEY1_GPIO_PORT,KEY1_PIN) == KEY_ON )/* K1 被按下 */
{
printf("挂起LED1线程\n");
uwRet = rt_thread_suspend(OFlowCal_thread);/* 挂起LED1线程 */
if(RT_EOK == uwRet)
{
rt_kprintf("挂起LED1线程成功\n");
}
else
{
rt_kprintf("挂起LED1线程失败失败代码0x%lx\n",uwRet);
}
}
if( Key_Scan(KEY2_GPIO_PORT,KEY2_PIN) == KEY_ON )/* K1 被按下 */
{
printf("恢复LED1线程\n");
uwRet = rt_thread_resume(OFlowCal_thread);/* 恢复LED1线程 */
if(RT_EOK == uwRet)
{
rt_kprintf("恢复LED1线程成功\n");
}
else
{
rt_kprintf("恢复LED1线程失败失败代码0x%lx\n",uwRet);
}
}
rt_thread_delay(20);
}
}
/********************************END OF FILE****************************/