commit c31b82ef9c2e275cac1f89d90905800f3608dfc5 Author: ldeyun Date: Tue Jul 8 00:11:32 2025 +0000 NG diff --git a/NG/Detail.c b/NG/Detail.c new file mode 100644 index 0000000..e05724d --- /dev/null +++ b/NG/Detail.c @@ -0,0 +1,1166 @@ +#include "NGCal.h" +#include "Detail.h" +#include +#include +#include +#define NUMBEROFCOMPONENTS 21 + +Detail *Detail_Construct(void) { + Detail *pDetail = (Detail *) malloc(sizeof(Detail)); + 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; + } + return pDetail; +} + +void Detail_Destroy(Detail *pDetail) { + if (pDetail) { + free(pDetail); + } +} + +int Detail_compositionchange(Detail *pDetail, NGParSTRUCT *ptNGPar) { + double dMixID = 0.0; + int i; + for (i = 0; i < NUMBEROFCOMPONENTS; i++) + dMixID += ((i + 2) * ptNGPar->adMixture[i]); + if (dMixID != pDetail->dOldMixID) { + pDetail->dOldMixID = dMixID; + return true; + } else { + return false; + } +} + +void Detail_Run(Detail *pDetail, NGParSTRUCT *ptNGPar) { + int i; + bool bCompChange = Detail_compositionchange(pDetail, ptNGPar); + ptNGPar->bForceUpdate = ptNGPar->bForceUpdate || bCompChange; + if (ptNGPar->bForceUpdate) { + pDetail->iNCC = 0; + for (i = 0; i < NUMBEROFCOMPONENTS; i++) { + if (ptNGPar->adMixture[i] > 0.0) { + pDetail->aiCID[pDetail->iNCC] = i; + pDetail->dXi[pDetail->iNCC] = ptNGPar->adMixture[i]; + pDetail->iNCC++; + } + } + 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)) { + pDetail->dP = ptNGPar->dPb * 1.0e-6; + pDetail->dT = ptNGPar->dTb; + Detail_temp(pDetail); + Detail_ddetail(pDetail, ptNGPar); + ptNGPar->dDb = pDetail->dRho; + ptNGPar->dZb = Detail_zdetail(pDetail, pDetail->dRho); + pDetail->dRhoTP = (pDetail->dP * ptNGPar->dMrx) / (ptNGPar->dZb * RGASKJ * pDetail->dT); + Detail_relativedensity(pDetail, ptNGPar); + ptNGPar->dRhob = pDetail->dRhoTP; + pDetail->dOldTb = ptNGPar->dTb; + pDetail->dOldPb = ptNGPar->dPb; + 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 = true; + } + if ((fabs(ptNGPar->dPf - pDetail->dOldPf) > P_CHG_TOL) || (ptNGPar->bForceUpdate)) { + Detail_ddetail(pDetail, ptNGPar); + ptNGPar->dDf = pDetail->dRho; + ptNGPar->dZf = Detail_zdetail(pDetail, pDetail->dRho); + pDetail->dRhoTP = (pDetail->dP * ptNGPar->dMrx) / (ptNGPar->dZf * RGASKJ * pDetail->dT); + ptNGPar->dRhof = pDetail->dRhoTP; + pDetail->dOldTf = ptNGPar->dTf; + pDetail->dOldPf = ptNGPar->dPf; + } + if ((ptNGPar->dZb > 0.0) && (ptNGPar->dZf > 0.0)) { + ptNGPar->dFpv = sqrt(ptNGPar->dZb / ptNGPar->dZf); + } else { + ptNGPar->lStatus = GENERAL_CALCULATION_FAILURE; + ptNGPar->bForceUpdate = false; + } +} + +int Detail_table(Detail *pDetail) { + if (!pDetail) { + return -1; + } + + double adAn[58]; + double adUn[58]; + adAn[0] = 0.153832600; + adAn[1] = 1.341953000; + adAn[2] = -2.998583000; + adAn[3] = -0.048312280; + adAn[4] = 0.375796500; + adAn[5] = -1.589575000; + adAn[6] = -0.053588470; + adAn[7] = 0.886594630; + adAn[8] = -0.710237040; + adAn[9] = -1.471722000; + adAn[10] = 1.321850350; + adAn[11] = -0.786659250; + adAn[12] = 2.29129E-09; + adAn[13] = 0.157672400; + adAn[14] = -0.436386400; + adAn[15] = -0.044081590; + adAn[16] = -0.003433888; + adAn[17] = 0.032059050; + adAn[18] = 0.024873550; + adAn[19] = 0.073322790; + adAn[20] = -0.001600573; + adAn[21] = 0.642470600; + adAn[22] = -0.416260100; + adAn[23] = -0.066899570; + adAn[24] = 0.279179500; + adAn[25] = -0.696605100; + adAn[26] = -0.002860589; + adAn[27] = -0.008098836; + adAn[28] = 3.150547000; + adAn[29] = 0.007224479; + adAn[30] = -0.705752900; + adAn[31] = 0.534979200; + adAn[32] = -0.079314910; + adAn[33] = -1.418465000; + adAn[34] = -5.99905E-17; + adAn[35] = 0.105840200; + adAn[36] = 0.034317290; + adAn[37] = -0.007022847; + adAn[38] = 0.024955870; + adAn[39] = 0.042968180; + adAn[40] = 0.746545300; + adAn[41] = -0.291961300; + adAn[42] = 7.294616000; + adAn[43] = -9.936757000; + adAn[44] = -0.005399808; + adAn[45] = -0.243256700; + adAn[46] = 0.049870160; + adAn[47] = 0.003733797; + adAn[48] = 1.874951000; + adAn[49] = 0.002168144; + adAn[50] = -0.658716400; + adAn[51] = 0.000205518; + adAn[52] = 0.009776195; + adAn[53] = -0.020487080; + adAn[54] = 0.015573220; + adAn[55] = 0.006862415; + adAn[56] = -0.001226752; + adAn[57] = 0.002850908; + + + adUn[0] = 0.0; + adUn[1] = 0.5; + adUn[2] = 1.0; + adUn[3] = 3.5; + adUn[4] = -0.5; + adUn[5] = 4.5; + adUn[6] = 0.5; + adUn[7] = 7.5; + adUn[8] = 9.5; + adUn[9] = 6.0; + adUn[10] = 12.0; + adUn[11] = 12.5; + adUn[12] = -6.0; + adUn[13] = 2.0; + adUn[14] = 3.0; + adUn[15] = 2.0; + adUn[16] = 2.0; + adUn[17] = 11.0; + adUn[18] = -0.5; + adUn[19] = 0.5; + adUn[20] = 0.0; + adUn[21] = 4.0; + adUn[22] = 6.0; + adUn[23] = 21.0; + adUn[24] = 23.0; + adUn[25] = 22.0; + adUn[26] = -1.0; + adUn[27] = -0.5; + adUn[28] = 7.0; + adUn[29] = -1.0; + adUn[30] = 6.0; + adUn[31] = 4.0; + adUn[32] = 1.0; + adUn[33] = 9.0; + adUn[34] = -13.0; + adUn[35] = 21.0; + adUn[36] = 8.0; + adUn[37] = -0.5; + adUn[38] = 0.0; + adUn[39] = 2.0; + adUn[40] = 7.0; + adUn[41] = 9.0; + adUn[42] = 22.0; + adUn[43] = 23.0; + adUn[44] = 1.0; + adUn[45] = 9.0; + adUn[46] = 3.0; + adUn[47] = 8.0; + adUn[48] = 23.0; + adUn[49] = 1.5; + adUn[50] = 5.0; + adUn[51] = -0.5; + adUn[52] = 4.0; + adUn[53] = 7.0; + adUn[54] = 3.0; + 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; + pDetail->adTable6Uij[j][k] = 1.0; + pDetail->adTable6Kij[j][k] = 1.0; + pDetail->adTable6Gij[j][k] = 1.0; + } + } + pDetail->adTable6Eij[0][1] = 0.971640; + pDetail->adTable6Eij[0][2] = 0.960644; + pDetail->adTable6Eij[0][4] = 0.994635; + pDetail->adTable6Eij[0][5] = 0.708218; + pDetail->adTable6Eij[0][6] = 0.931484; + pDetail->adTable6Eij[0][7] = 1.170520; + pDetail->adTable6Eij[0][8] = 0.990126; + pDetail->adTable6Eij[0][10] = 1.019530; + pDetail->adTable6Eij[0][11] = 0.989844; + pDetail->adTable6Eij[0][12] = 1.002350; + pDetail->adTable6Eij[0][13] = 0.999268; + pDetail->adTable6Eij[0][14] = 1.107274; + pDetail->adTable6Eij[0][15] = 0.880880; + pDetail->adTable6Eij[0][16] = 0.880973; + pDetail->adTable6Eij[0][17] = 0.881067; + pDetail->adTable6Eij[0][18] = 0.881161; + pDetail->adTable6Eij[1][2] = 1.022740; + pDetail->adTable6Eij[1][3] = 0.970120; + pDetail->adTable6Eij[1][4] = 0.945939; + pDetail->adTable6Eij[1][5] = 0.746954; + pDetail->adTable6Eij[1][6] = 0.902271; + pDetail->adTable6Eij[1][7] = 1.086320; + pDetail->adTable6Eij[1][8] = 1.005710; + pDetail->adTable6Eij[1][9] = 1.021000; + pDetail->adTable6Eij[1][10] = 0.946914; + pDetail->adTable6Eij[1][11] = 0.973384; + pDetail->adTable6Eij[1][12] = 0.959340; + pDetail->adTable6Eij[1][13] = 0.945520; + pDetail->adTable6Eij[2][3] = 0.925053; + pDetail->adTable6Eij[2][4] = 0.960237; + pDetail->adTable6Eij[2][5] = 0.849408; + pDetail->adTable6Eij[2][6] = 0.955052; + pDetail->adTable6Eij[2][7] = 1.281790; + pDetail->adTable6Eij[2][8] = 1.500000; + pDetail->adTable6Eij[2][10] = 0.906849; + pDetail->adTable6Eij[2][11] = 0.897362; + pDetail->adTable6Eij[2][12] = 0.726255; + pDetail->adTable6Eij[2][13] = 0.859764; + pDetail->adTable6Eij[2][14] = 0.855134; + pDetail->adTable6Eij[2][15] = 0.831229; + pDetail->adTable6Eij[2][16] = 0.808310; + pDetail->adTable6Eij[2][17] = 0.786323; + pDetail->adTable6Eij[2][18] = 0.765171; + pDetail->adTable6Eij[3][4] = 1.022560; + pDetail->adTable6Eij[3][5] = 0.693168; + pDetail->adTable6Eij[3][6] = 0.946871; + pDetail->adTable6Eij[3][7] = 1.164460; + pDetail->adTable6Eij[3][11] = 1.013060; + pDetail->adTable6Eij[3][13] = 1.005320; + pDetail->adTable6Eij[4][7] = 1.034787; + pDetail->adTable6Eij[4][11] = 1.004900; + pDetail->adTable6Eij[6][14] = 1.008692; + pDetail->adTable6Eij[6][15] = 1.010126; + pDetail->adTable6Eij[6][16] = 1.011501; + pDetail->adTable6Eij[6][17] = 1.012821; + pDetail->adTable6Eij[6][18] = 1.014089; + pDetail->adTable6Eij[7][8] = 1.100000; + pDetail->adTable6Eij[7][10] = 1.300000; + pDetail->adTable6Eij[7][11] = 1.300000; + pDetail->adTable6Uij[0][1] = 0.886106; + pDetail->adTable6Uij[0][2] = 0.963827; + pDetail->adTable6Uij[0][4] = 0.990877; + pDetail->adTable6Uij[0][6] = 0.736833; + pDetail->adTable6Uij[0][7] = 1.156390; + pDetail->adTable6Uij[0][11] = 0.992291; + pDetail->adTable6Uij[0][13] = 1.003670; + pDetail->adTable6Uij[0][14] = 1.302576; + pDetail->adTable6Uij[0][15] = 1.191904; + pDetail->adTable6Uij[0][16] = 1.205769; + pDetail->adTable6Uij[0][17] = 1.219634; + pDetail->adTable6Uij[0][18] = 1.233498; + pDetail->adTable6Uij[1][2] = 0.835058; + pDetail->adTable6Uij[1][3] = 0.816431; + pDetail->adTable6Uij[1][4] = 0.915502; + pDetail->adTable6Uij[1][6] = 0.993476; + pDetail->adTable6Uij[1][7] = 0.408838; + pDetail->adTable6Uij[1][11] = 0.993556; + pDetail->adTable6Uij[2][3] = 0.969870; + pDetail->adTable6Uij[2][6] = 1.045290; + pDetail->adTable6Uij[2][8] = 0.900000; + pDetail->adTable6Uij[2][14] = 1.066638; + pDetail->adTable6Uij[2][15] = 1.077634; + pDetail->adTable6Uij[2][16] = 1.088178; + pDetail->adTable6Uij[2][17] = 1.098291; + pDetail->adTable6Uij[2][18] = 1.108021; + pDetail->adTable6Uij[3][4] = 1.065173; + pDetail->adTable6Uij[3][6] = 0.971926; + pDetail->adTable6Uij[3][7] = 1.616660; + pDetail->adTable6Uij[3][10] = 1.250000; + pDetail->adTable6Uij[3][11] = 1.250000; + pDetail->adTable6Uij[3][12] = 1.250000; + pDetail->adTable6Uij[3][13] = 1.250000; + pDetail->adTable6Uij[6][14] = 1.028973; + pDetail->adTable6Uij[6][15] = 1.033754; + pDetail->adTable6Uij[6][16] = 1.038338; + pDetail->adTable6Uij[6][17] = 1.042735; + pDetail->adTable6Uij[6][18] = 1.046966; + pDetail->adTable6Kij[0][1] = 1.003630; + pDetail->adTable6Kij[0][2] = 0.995933; + pDetail->adTable6Kij[0][4] = 1.007619; + pDetail->adTable6Kij[0][6] = 1.000080; + pDetail->adTable6Kij[0][7] = 1.023260; + pDetail->adTable6Kij[0][11] = 0.997596; + pDetail->adTable6Kij[0][13] = 1.002529; + pDetail->adTable6Kij[0][14] = 0.982962; + pDetail->adTable6Kij[0][15] = 0.983565; + pDetail->adTable6Kij[0][16] = 0.982707; + pDetail->adTable6Kij[0][17] = 0.981849; + pDetail->adTable6Kij[0][18] = 0.980991; + pDetail->adTable6Kij[1][2] = 0.982361; + pDetail->adTable6Kij[1][3] = 1.007960; + pDetail->adTable6Kij[1][6] = 0.942596; + pDetail->adTable6Kij[1][7] = 1.032270; + pDetail->adTable6Kij[2][3] = 1.008510; + pDetail->adTable6Kij[2][6] = 1.007790; + pDetail->adTable6Kij[2][14] = 0.910183; + pDetail->adTable6Kij[2][15] = 0.895362; + pDetail->adTable6Kij[2][16] = 0.881152; + pDetail->adTable6Kij[2][17] = 0.867520; + pDetail->adTable6Kij[2][18] = 0.854406; + pDetail->adTable6Kij[3][4] = 0.986893; + pDetail->adTable6Kij[3][6] = 0.999969; + pDetail->adTable6Kij[3][7] = 1.020340; + pDetail->adTable6Kij[6][14] = 0.968130; + pDetail->adTable6Kij[6][15] = 0.962870; + pDetail->adTable6Kij[6][16] = 0.957828; + pDetail->adTable6Kij[6][17] = 0.952441; + pDetail->adTable6Kij[6][18] = 0.948338; + pDetail->adTable6Gij[0][2] = 0.807653; + pDetail->adTable6Gij[0][7] = 1.957310; + pDetail->adTable6Gij[1][2] = 0.982746; + pDetail->adTable6Gij[2][3] = 0.370296; + pDetail->adTable6Gij[2][5] = 1.673090; + + 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 + }, + { + 891.56, 0, 0, 1562.14, 2221.1, 44.433, 562.38, 286.15, 282.91, 0, 2870.58, 2879.76, 3531.68, 3538.6, + 4198.24, 4857.18, 5516.01, 6175.82, 6834.9, 0, 0 + }, + { + 891.09, 0, 0, 1561.41, 2220.13, 44.224, 562.19, 285.99, 282.95, 0, 2869.38, 2878.57, 3530.24, 3537.17, + 4196.58, 4855.29, 5513.88, 6173.46, 6832.31, 0, 0 + }, + { + 890.63, 0, 0, 1560.69, 2219.17, 44.016, 562.01, 285.83, 282.98, 0, 2868.2, 2877.4, 3528.83, 3535.77, + 4194.95, 4853.43, 5511.8, 6171.15, 6829.77, 0, 0 + } + }; + memcpy(pDetail->adTableHhvMol, adTableHhvMol, sizeof(adTableHhvMol)); + 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 + }, + { + 802.69, 0, 0, 1428.84, 2043.37, 0, 517.95, 241.72, 282.91, 0, 2648.42, 2657.6, 3265.08, 3272, 3887.21, + 4501.72, 5116.11, 5731.49, 6346.14, 0, 0 + }, + { + 802.65, 0, 0, 1428.74, 2043.23, 0, 517.97, 241.76, 282.95, 0, 2648.26, 2657.45, 3264.89, 3271.83, 3887.01, + 4501.49, 5115.87, 5731.22, 6345.85, 0, 0 + }, + { + 802.6, 0, 0, 1428.64, 2043.11, 0, 517.99, 241.81, 282.98, 0, 2648.12, 2657.32, 3264.73, 3271.67, 3886.84, + 4501.3, 5115.66, 5730.99, 6345.59, 0, 0 + } + }; + memcpy(pDetail->adTableLhvMol, adTableLhvMol, sizeof(adTableLhvMol)); + return 0; +} + +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; + pDetail->adTable5Si[j] = 0.0; + pDetail->adTable5Wi[j] = 0.0; + } + pDetail->adTable5Qi[2] = 0.690000; + pDetail->adTable5Qi[5] = 1.067750; + pDetail->adTable5Qi[6] = 0.633276; + pDetail->adTable5Fi[7] = 1.0000; + pDetail->adTable5Si[5] = 1.5822; + pDetail->adTable5Si[6] = 0.3900; + pDetail->adTable5Wi[5] = 1.0000; + for (int j = pDetail->iNCC - 1; j >= 0; j--) { + pDetail->dMri[j] = adTable5Mri[pDetail->aiCID[j]]; + pDetail->dKi[j] = adTable5Ki[pDetail->aiCID[j]]; + } + for (int j = 0; j < pDetail->iNCC; j++) { + pDetail->dGi[j] = adTable5Gi[pDetail->aiCID[j]]; + pDetail->dEi[j] = adTable5Ei[pDetail->aiCID[j]]; + } + for (int j = 0; j < pDetail->iNCC; j++) { + pDetail->dQi[j] = pDetail->adTable5Qi[pDetail->aiCID[j]]; + pDetail->dFi[j] = 0.0; + if (pDetail->aiCID[j] == 7) pDetail->dFi[j] = pDetail->adTable5Fi[7]; + pDetail->dSi[j] = pDetail->adTable5Si[pDetail->aiCID[j]]; + pDetail->dWi[j] = pDetail->adTable5Wi[pDetail->aiCID[j]]; + } + for (int j = 0; j < pDetail->iNCC; j++) { + for (int k = j; k < pDetail->iNCC; k++) { + pDetail->dUij[j][k] = pDetail->adTable6Uij[pDetail->aiCID[j]][pDetail->aiCID[k]]; + pDetail->dKij[j][k] = pDetail->adTable6Kij[pDetail->aiCID[j]][pDetail->aiCID[k]]; + pDetail->dEij[j][k] = pDetail->adTable6Eij[pDetail->aiCID[j]][pDetail->aiCID[k]]; + pDetail->dGij[j][k] = pDetail->adTable6Gij[pDetail->aiCID[j]][pDetail->aiCID[k]]; + } + } +} + +void Detail_chardl(Detail *pDetail, NGParSTRUCT *ptNGPar) +{ + double tmfrac = 0.0; + + for (int j = 0; j < pDetail->iNCC; j++) { + tmfrac += pDetail->dXi[j]; + } + for (int j = 0; j < pDetail->iNCC; j++) { + pDetail->dXi[j] /= tmfrac; + } + for (int j = 0; j < 18; j++) { + pDetail->adBcoef[j] = 0.0; + } + double k5p0 = 0.0, k2p5 = 0.0, u5p0 = 0.0, u2p5 = 0.0; + pDetail->dW = 0.0; + double q1p0 = 0.0; + pDetail->dF = 0.0; + ptNGPar->dMrx = 0.0; + ptNGPar->dHhvMol = 0.0; + ptNGPar->dLhvMol = 0.0; + for (int i = 0; i < pDetail->iNCC; i++) { + ptNGPar->dMrx += pDetail->dXi[i] * pDetail->dMri[i]; + switch (ptNGPar->dCbtj) { + case 2: + ptNGPar->dHhvMol += pDetail->adTableHhvMol[0][i] * ptNGPar->adMixture[i]; + ptNGPar->dLhvMol += pDetail->adTableHhvMol[0][i] * ptNGPar->adMixture[i]; + break; + case 1: + ptNGPar->dHhvMol += pDetail->adTableHhvMol[1][i] * ptNGPar->adMixture[i]; + ptNGPar->dLhvMol += pDetail->adTableLhvMol[1][i] * ptNGPar->adMixture[i]; + break; + case 0: + ptNGPar->dHhvMol += pDetail->adTableHhvMol[2][i] * ptNGPar->adMixture[i]; + ptNGPar->dLhvMol += pDetail->adTableLhvMol[2][i] * ptNGPar->adMixture[i]; + break; + } + k2p5 += pDetail->dXi[i] * pow(pDetail->dKi[i], 2.5); + u2p5 += pDetail->dXi[i] * pow(pDetail->dEi[i], 2.5); + pDetail->dW += pDetail->dXi[i] * pDetail->dGi[i]; + q1p0 += pDetail->dXi[i] * pDetail->dQi[i]; + pDetail->dF += pDetail->dXi[i] * pDetail->dXi[i] * pDetail->dFi[i]; + + for (int j = i; j < pDetail->iNCC; 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) { + 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) { + 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) { + double avgG = (pDetail->dGi[i] + pDetail->dGi[j]) / 2.0; + pDetail->dW += Xij * (pDetail->dGij[i][j] - 1.0) * avgG; + } + 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; + pDetail->adBcoef[3] += s3 * e3p5; + pDetail->adBcoef[4] += s3 * Gij / e0p5; + pDetail->adBcoef[5] += s3 * Gij * e4p5; + pDetail->adBcoef[6] += s3 * pDetail->dQi[i] * pDetail->dQi[j] * e0p5; + pDetail->adBcoef[7] += s3 * pDetail->dSi[i] * pDetail->dSi[j] * e7p5; + pDetail->adBcoef[8] += s3 * pDetail->dSi[i] * pDetail->dSi[j] * e9p5; + pDetail->adBcoef[9] += s3 * pDetail->dWi[i] * pDetail->dWi[j] * e6p0; + pDetail->adBcoef[10] += s3 * pDetail->dWi[i] * pDetail->dWi[j] * e12p0; + pDetail->adBcoef[11] += s3 * pDetail->dWi[i] * pDetail->dWi[j] * e12p5; + pDetail->adBcoef[12] += s3 * pDetail->dFi[i] * pDetail->dFi[j] / e6p0; + pDetail->adBcoef[13] += s3 * e2p0; + pDetail->adBcoef[14] += s3 * e3p0; + pDetail->adBcoef[15] += s3 * pDetail->dQi[i] * pDetail->dQi[j] * e2p0; + pDetail->adBcoef[16] += s3 * e2p0; + pDetail->adBcoef[17] += s3 * e11p0; + } + } + + 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]; + } + pDetail->dKp3 = pow(k5p0 + pow(k2p5, 2.0), 0.6); + pDetail->dU = pow(u5p0 + pow(u2p5, 2.0), 0.2); + pDetail->dQp2 = q1p0 * q1p0; +} + +void Detail_bvir(Detail *pDetail) { + pDetail->dB = pDetail->ddBdT = pDetail->dd2BdT2 = 0.0; + 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; + Bx[2] = pDetail->adBcoef[2] / t; + Bx[3] = pDetail->adBcoef[3] / t3p5; + Bx[4] = pDetail->adBcoef[4] * t0p5; + Bx[5] = pDetail->adBcoef[5] / t4p5; + Bx[6] = pDetail->adBcoef[6] / t0p5; + Bx[7] = pDetail->adBcoef[7] / t7p5; + Bx[8] = pDetail->adBcoef[8] / t9p5; + Bx[9] = pDetail->adBcoef[9] / t6p0; + Bx[10] = pDetail-> adBcoef[10] / t12p0; + Bx[11] = pDetail-> adBcoef[11] / t12p5; + Bx[12] = pDetail-> adBcoef[12] * t6p0; + Bx[13] = pDetail-> adBcoef[13] / t2p0; + Bx[14] = pDetail-> adBcoef[14] / t3p0; + Bx[15] = pDetail-> adBcoef[15] / t2p0; + Bx[16] = pDetail-> adBcoef[16] / t2p0; + Bx[17] = pDetail-> adBcoef[17] / t11p0; + for (int i = 0; i < 18; i++) { + pDetail->dB += Bx[i]; + } + for (int i = 0; i < 18; i++) { + if (pDetail->adUn[i] != 0.0) { + Bx[i] *= pDetail->adUn[i] ; + } + } + for (int i = 0; i < 18; i++) { + if (pDetail->adUn[i] != 0.0) { + pDetail->ddBdT += Bx[i]/t; + } + } + pDetail->ddBdT = -pDetail->ddBdT; + for (int i = 0; i < 18; i++) { + if (pDetail->adUn[i] != 0.0 && pDetail->adUn[i] != -1.0) { + Bx[i] *= (pDetail->adUn[i] + 1.0) ; + } + } + for (int i = 0; i < 18; i++) { + if (pDetail->adUn[i] != 0.0 && pDetail->adUn[i] != -1.0) { + pDetail->dd2BdT2 += Bx[i]/ t2p0; + } + } +} + +void Detail_temp(Detail *pDetail) { + Detail_bvir(pDetail); + 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; + pDetail->adFn[14] = pDetail->adAn[14] / tr3p0; + pDetail->adFn[15] = pDetail->adAn[15] * pDetail->dQp2 / tr2p0; + pDetail->adFn[16] = pDetail->adAn[16] / tr2p0; + pDetail->adFn[17] = pDetail->adAn[17] / tr11p0; + pDetail->adFn[18] = pDetail->adAn[18] * tr0p5; + pDetail->adFn[19] = pDetail->adAn[19] / tr0p5; + pDetail->adFn[20] = pDetail->adAn[20]; + + pDetail->adFn[21] = pDetail->adAn[21] / tr4p0; + pDetail->adFn[22] = pDetail->adAn[22] / tr6p0; + pDetail->adFn[23] = pDetail->adAn[23] / tr21p0; + pDetail->adFn[24] = pDetail->adAn[24] * pDetail->dW / tr23p0; + pDetail->adFn[25] = pDetail->adAn[25] * pDetail->dQp2 / tr22p0; + pDetail->adFn[26] = pDetail->adAn[26] * pDetail->dF * tr; + pDetail->adFn[27] = pDetail->adAn[27] * pDetail->dQp2 * tr0p5; + pDetail->adFn[28] = pDetail->adAn[28] * pDetail->dW / tr7p0; + pDetail->adFn[29] = pDetail->adAn[29] * pDetail->dF * tr; + pDetail->adFn[30] = pDetail->adAn[30] / tr6p0; + pDetail->adFn[31] = pDetail->adAn[31] * pDetail->dW / tr4p0; + pDetail->adFn[32] = pDetail->adAn[32] * pDetail->dW / tr; + pDetail->adFn[33] = pDetail->adAn[33] * pDetail->dW / tr9p0; + pDetail->adFn[34] = pDetail->adAn[34] * pDetail->dF * tr13p0; + pDetail->adFn[35] = pDetail->adAn[35] / tr21p0; + pDetail->adFn[36] = pDetail->adAn[36] * pDetail->dQp2 / tr8p0; + pDetail->adFn[37] = pDetail->adAn[37] * tr0p5; + pDetail->adFn[38] = pDetail->adAn[38]; + pDetail->adFn[39] = pDetail->adAn[39] / tr2p0; + pDetail->adFn[40] = pDetail->adAn[40] / tr7p0; + pDetail->adFn[41] = pDetail->adAn[41] * pDetail->dQp2 / tr9p0; + pDetail->adFn[42] = pDetail->adAn[42] / tr22p0; + pDetail->adFn[43] = pDetail->adAn[43] / tr23p0; + pDetail->adFn[44] = pDetail->adAn[44] / tr; + pDetail->adFn[45] = pDetail->adAn[45] / tr9p0; + pDetail->adFn[46] = pDetail->adAn[46] * pDetail->dQp2 / tr3p0; + pDetail->adFn[47] = pDetail->adAn[47] / tr8p0; + pDetail->adFn[48] = pDetail->adAn[48] * pDetail->dQp2 / tr23p0; + pDetail->adFn[49] = pDetail->adAn[49] / tr1p5; + pDetail->adFn[50] = pDetail->adAn[50] * pDetail->dW / tr5p0; + pDetail->adFn[51] = pDetail->adAn[51] * pDetail->dQp2 * tr0p5; + pDetail->adFn[52] = pDetail->adAn[52] / tr4p0; + pDetail->adFn[53] = pDetail->adAn[53] * pDetail->dW / tr7p0; + pDetail->adFn[54] = pDetail->adAn[54] / tr3p0; + pDetail->adFn[55] = pDetail->adAn[55] * pDetail->dW; + pDetail->adFn[56] = pDetail->adAn[56] / tr; + pDetail->adFn[57] = pDetail->adAn[57] * pDetail->dQp2; +} + +void Detail_ddetail(Detail *pDetail, NGParSTRUCT *ptNGPar) { + 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; + } + 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; + delx = x1 - x2; + delprv = delx; + } + if (fabs(y3) < fabs(y2)) { + x1 = x2; + x2 = x3; + x3 = x1; + y1 = y2; + y2 = y3; + y3 = y1; + } + 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) { + y2my3 = y2 - y3; + y3my1 = y3 - y1; + y1my2 = y1 - y2; + xdenom = -(y1my2) * (y2my3) * (y3my1); + xnumer = x1 * y2 * y3 * (y2my3) + + x2 * y3 * y1 * (y3my1) + + x3 * y1 * y2 * (y1my2) - x2 * xdenom; + } else { + xnumer = (x2 - x1) * y2; + xdenom = y1 - y2; + } + if (2.0 * fabs(xnumer) < fabs(delprv * xdenom)) { + delprv = delx; + delx = xnumer / xdenom; + } else { + delx = delbis; + delprv = delbis; + } + } + if ((fabs(y2) < epsp * pDetail->dP) && (fabs(delx) < epsr * fabs(x2))) { + pDetail->dRho = x2 + delx; + return; + } + if (fabs(delx) < delmin) { + sgndel = delbis / fabs(delbis); + delx = 1.0000009 * sgndel * delmin; + delprv = delx; + } + boundn = delx * (x2 + delx - x3); + if (boundn > 0.0) { + delx = delbis; + delprv = delbis; + } + x1 = x2; + y1 = y2; + x2 = x2 + delx; + Detail_pdetail(pDetail, x2); + y2 = pDetail->dPCalc - pDetail->dP; + } + ptNGPar->lStatus = MAX_NUM_OF_ITERATIONS_EXCEEDED; + pDetail->dRho = x2; +} + +void Detail_braket(Detail *pDetail, NGParSTRUCT *ptNGPar) { + 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; + videal = RGASKJ * pDetail->dT / pDetail->dP; + if (fabs(pDetail->dB) < (0.167 * videal)) { + rho2 = 0.95 / (videal + pDetail->dB); + } else { + rho2 = 1.15 / videal; + } + 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; + rho2 = rho1 + del; + continue; + } + Detail_pdetail(pDetail, rho2); + p2 = pDetail->dPCalc; + if (p2 > pDetail->dP) { + pDetail->dRhoL = rho1; + pDetail->dPRhoL = p1; + pDetail->dRhoH = rho2; + pDetail->dPRhoH = p2; + ptNGPar->lStatus = NORMAL; + return; + } else if (p2 > p1) { + if (ptNGPar->lStatus == MAX_DENSITY_IN_BRAKET_EXCEEDED) + del *= 2.0; + rho1 = rho2; + p1 = p2; + rho2 = rho1 + del; + continue; + } else { + ptNGPar->lStatus = NEGATIVE_DENSITY_DERIVATIVE; + pDetail->dRho = rho1; + return; + } + } + ptNGPar->lStatus = MAX_NUM_OF_ITERATIONS_EXCEEDED; + pDetail->dRho = rho2; + return; +} + +void Detail_pdetail(Detail *pDetail, double dD) { + pDetail->dPCalc = Detail_zdetail(pDetail, dD) * dD * RGASKJ * pDetail->dT; +} + +double Detail_zdetail(Detail *pDetail, double d) { + 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) + + (pDetail->adFn[16] + pDetail->adFn[17]) * D1 * (exp4 - 1.0 - 4.0 * D4 * exp4) + + (pDetail->adFn[18] + pDetail->adFn[19]) * D2 * 2.0 + + (pDetail->adFn[20] + pDetail->adFn[21] + pDetail->adFn[22]) * D2 * (2.0 - 2.0 * D2) * exp2 + + (pDetail->adFn[23] + pDetail->adFn[24] + pDetail->adFn[25]) * D2 * (2.0 - 4.0 * D4) * exp4 + + pDetail->adFn[26] * D2 * (2.0 - 4.0 * D4) * exp4 + + pDetail->adFn[27] * D3 * 3.0 + + (pDetail->adFn[28] + pDetail->adFn[29]) * D3 * (3.0 - D1) * exp1 + + (pDetail->adFn[30] + pDetail->adFn[31]) * D3 * (3.0 - 2.0 * D2) * exp2 + + (pDetail->adFn[32] + pDetail->adFn[33]) * D3 * (3.0 - 3.0 * D3) * exp3 + + (pDetail->adFn[34] + pDetail->adFn[35] + pDetail->adFn[36]) * D3 * (3.0 - 4.0 * D4) * exp4 + + (pDetail->adFn[37] + pDetail->adFn[38]) * D4 * 4.0 + + (pDetail->adFn[39] + pDetail->adFn[40] + pDetail->adFn[41]) * D4 * (4.0 - 2.0 * D2) * exp2 + + (pDetail->adFn[42] + pDetail->adFn[43]) * D4 * (4.0 - 4.0 * D4) * exp4 + + pDetail->adFn[44] * D5 * 5.0 + + (pDetail->adFn[45] + pDetail->adFn[46]) * D5 * (5.0 - 2.0 * D2) * exp2 + + (pDetail->adFn[47] + pDetail->adFn[48]) * D5 * (5.0 - 4.0 * D4) * exp4 + + pDetail->adFn[49] * D6 * 6.0 + + pDetail->adFn[50] * D6 * (6.0 - 2.0 * D2) * exp2 + + pDetail->adFn[51] * D7 * 7.0 + + pDetail->adFn[52] * D7 * (7.0 - 2.0 * D2) * exp2 + + pDetail->adFn[53] * D8 * (8.0 - D1) * exp1 + + (pDetail->adFn[54] + pDetail->adFn[55]) * D8 * (8.0 - 2.0 * D2) * exp2 + + (pDetail->adFn[56] + pDetail->adFn[57]) * D9 * (9.0 - 2.0 * D2) * exp2; + return pDetail->dZ; +} + +double Detail_dZdT(Detail *pDetail, double d) { + 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 { + pDetail->fx[i] = 0.0; + } + } + pDetail->ddZdT = d * pDetail->ddBdT; + if (pDetail->dF) + pDetail->ddZdT += pDetail->fx[12] - (pDetail->fx[12] * (1.0 - 3.0 * D3) * exp3); + 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); + tmp = (1.0 - 4.0 * D4) * exp4; + pDetail->ddZdT += pDetail->fx[16] - (pDetail->fx[16] * tmp); + pDetail->ddZdT += pDetail->fx[17] - (pDetail->fx[17] * tmp); + pDetail->ddZdT = pDetail->ddZdT - (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 + - pDetail->fx[26] * D1 * (2.0 - 4.0 * D4) * exp4 + - pDetail->fx[27] * D2 * 3.0 + - (pDetail->fx[28] + pDetail->fx[29]) * D2 * (3.0 - D1) * exp1 + - (pDetail->fx[30] + pDetail->fx[31]) * D2 * (3.0 - 2.0 * D2) * exp2 + - (pDetail->fx[32] + pDetail->fx[33]) * D2 * (3.0 - 3.0 * D3) * exp3 + - (pDetail->fx[34] + pDetail->fx[35] + pDetail->fx[36]) * D2 * (3.0 - 4.0 * D4) * exp4 + - pDetail->fx[37] * D3 * 4.0 + - (pDetail->fx[39] + pDetail->fx[40] + pDetail->fx[41]) * D3 * (4.0 - 2.0 * D2) * exp2 + - (pDetail->fx[42] + pDetail->fx[43]) * D3 * (4.0 - 4.0 * D4) * exp4 + - pDetail->fx[44] * D4 * 5.0 + - (pDetail->fx[45] + pDetail->fx[46]) * D4 * (5.0 - 2.0 * D2) * exp2 + - (pDetail->fx[47] + pDetail->fx[48]) * D4 * (5.0 - 4.0 * D4) * exp4 + - pDetail->fx[49] * D5 * 6.0 + - pDetail->fx[50] * D5 * (6.0 - 2.0 * D2) * exp2 + - pDetail->fx[51] * D6 * 7.0 + - pDetail->fx[52] * D6 * (7.0 - 2.0 * D2) * exp2 + - pDetail->fx[53] * D7 * (8.0 - D1) * exp1 + - pDetail->fx[54] * D7 * (8.0 - 2.0 * D2) * exp2 + - pDetail->fx[56] * D8 * (9.0 - 2.0 * D2) * exp2; + return pDetail->ddZdT; +} + +double Detail_d2ZdT2(Detail *pDetail, double d) { + 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); + } else { + pDetail->fx[i] = 0.0; + } + } + pDetail->dd2ZdT2 = d * pDetail->dd2BdT2; + if (pDetail->dF) + pDetail->dd2ZdT2 += pDetail->fx[12] - (pDetail->fx[12] * (1.0 - 3.0 * D3) * exp3); + 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->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 + + pDetail->fx[26] * D1 * (2.0 - 4.0 * D4) * exp4 + + pDetail->fx[27] * D2 * 3.0 + + (pDetail->fx[28] + pDetail->fx[29]) * D2 * (3.0 - D1) * exp1 + + (pDetail->fx[30] + pDetail->fx[31]) * D2 * (3.0 - 2.0 * D2) * exp2 + + (pDetail->fx[32] + pDetail->fx[33]) * D2 * (3.0 - 3.0 * D3) * exp3 + + (pDetail->fx[34] + pDetail->fx[35] + pDetail->fx[36]) * D2 * (3.0 - 4.0 * D4) * exp4 + + pDetail->fx[37] * D3 * 4.0 + + (pDetail->fx[39] + pDetail->fx[40] + pDetail->fx[41]) * D3 * (4.0 - 2.0 * D2) * exp2 + + (pDetail->fx[42] + pDetail->fx[43]) * D3 * (4.0 - 4.0 * D4) * exp4 + + pDetail->fx[44] * D4 * 5.0 + + (pDetail->fx[45] + pDetail->fx[46]) * D4 * (5.0 - 2.0 * D2) * exp2 + + (pDetail->fx[47] + pDetail->fx[48]) * D4 * (5.0 - 4.0 * D4) * exp4 + + pDetail->fx[49] * D5 * 6.0 + + pDetail->fx[50] * D5 * (6.0 - 2.0 * D2) * exp2 + + pDetail->fx[51] * D6 * 7.0 + + pDetail->fx[52] * D6 * (7.0 - 2.0 * D2) * exp2 + + pDetail->fx[53] * D7 * (8.0 - D1) * exp1 + + pDetail->fx[54] * D7 * (8.0 - 2.0 * D2) * exp2 + + pDetail->fx[56] * D8 * (9.0 - 2.0 * D2) * exp2; + return pDetail->dd2ZdT2; +} + +double Detail_dZdD(Detail *pDetail, double d) { + double temp, temp1, temp2, temp3; + 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; + if (pDetail->dF) { + temp1 = -9.0 * D3 * exp3; + temp2 = (1.0 - 3.0 * D3) * exp3; + temp3 = -temp2 * 3.0 * D6; + temp = temp1 + temp2 + temp3; + pDetail->ddZdD += -pDetail->fx[12] + pDetail->fx[12] * temp; + } + temp1 = -4.0 * D2 * exp2; + temp2 = (1.0 - 2.0 * D2) * exp2; + temp3 = -temp2 * 2.0 * D2; + temp = temp1 + temp2 + temp3; + pDetail->ddZdD += -pDetail->fx[13] + pDetail->fx[13] * temp; + pDetail->ddZdD += -pDetail->fx[14] + pDetail->fx[14] * temp; + pDetail->ddZdD += -pDetail->fx[15] + pDetail->fx[15] * temp; + temp1 = -16.0 * D4 * exp4; + temp2 = (1.0 - 4.0 * D4) * exp4; + temp3 = -temp2 * 4.0 * D4; + temp = temp1 + temp2 + temp3; + pDetail->ddZdD += -pDetail->fx[16] + pDetail->fx[16] * temp; + pDetail->ddZdD += -pDetail->fx[17] + pDetail->fx[17] * temp; + temp = 4.0 * D1; + pDetail->ddZdD += pDetail->fx[18] * temp; + pDetail->ddZdD += pDetail->fx[19] * temp; + temp1 = -4.0 * D3 * exp2; + temp2 = (2.0 - 2.0 * D2) * 2.0 * D1 * exp2; + temp3 = -temp2 * D2; + temp = temp1 + temp2 + temp3; + pDetail->ddZdD += pDetail->fx[20] * temp; + pDetail->ddZdD += pDetail->fx[21] * temp; + pDetail->ddZdD += pDetail->fx[22] * temp; + temp1 = -16.0 * D5 * exp4; + temp2 = (2.0 - 4.0 * D4) * 2.0 * D1 * exp4; + temp3 = -temp2 * 2.0 * D4; + temp = temp1 + temp2 + temp3; + pDetail->ddZdD += pDetail->fx[23] * temp; + pDetail->ddZdD += pDetail->fx[24] * temp; + pDetail->ddZdD += pDetail->fx[25] * temp; + pDetail->ddZdD += pDetail->fx[26] * temp; + temp = 9.0 * D2; + pDetail->ddZdD += pDetail->fx[27] * temp; + temp = -D3 * exp1 + (3.0 - D1) * 3.0 * D2 * exp1; + temp -= (3.0 - D1) * D3 * exp1; + pDetail->ddZdD += pDetail->fx[28] * temp; + pDetail->ddZdD += pDetail->fx[29] * temp; + temp1 = -4.0 * D4 * exp2; + temp2 = (3.0 - 2.0 * D2) * 3.0 * D2 * exp2; + temp3 = -(3.0 - 2.0 * D2) * 2.0 * D4 * exp2; + temp = temp1 + temp2 + temp3; + pDetail->ddZdD += pDetail->fx[30] * temp; + pDetail->ddZdD += pDetail->fx[31] * temp; + temp1 = -9.0 * D5 * exp3; + temp2 = (3.0 - 3.0 * D3) * 3.0 * D2 * exp3; + temp3 = -(3.0 - 3.0 * D3) * 3.0 * D5 * exp3; + temp = temp1 + temp2 + temp3; + pDetail->ddZdD += pDetail->fx[32] * temp; + pDetail->ddZdD += pDetail->fx[33] * temp; + temp1 = -16.0 * D6 * exp4; + temp2 = (3.0 - 4.0 * D4) * 3.0 * D2 * exp4; + temp3 = -(3.0 - 4.0 * D4) * D6 * 4.0 * exp4; + temp = temp1 + temp2 + temp3; + pDetail->ddZdD += pDetail->fx[34] * temp; + pDetail->ddZdD += pDetail->fx[35] * temp; + pDetail->ddZdD += pDetail->fx[36] * temp; + temp = 16.0 * D3; + pDetail->ddZdD += pDetail->fx[37] * temp; + pDetail->ddZdD += pDetail->fx[38] * temp; + temp1 = -4.0 * D5 * exp2; + temp2 = (4.0 - 2.0 * D2) * 4.0 * D3 * exp2; + temp3 = -(4.0 - 2.0 * D2) * 2.0 * D5 * exp2; + temp = temp1 + temp2 + temp3; + pDetail->ddZdD += pDetail->fx[39] * temp; + pDetail->ddZdD += pDetail->fx[40] * temp; + pDetail->ddZdD += pDetail->fx[41] * temp; + temp = -16.0 * D7 * exp4 + (4.0 - 4.0 * D4) * 4.0 * D3 * exp4; + temp -= (4.0 - 4.0 * D4) * D7 * 4.0 * exp4; + pDetail->ddZdD += pDetail->fx[42] * temp; + pDetail->ddZdD += pDetail->fx[43] * temp; + temp = 25.0 * D4; + pDetail->ddZdD += pDetail->fx[44] * temp; + temp = -4.0 * D6 * exp2 + (5.0 - 2.0 * D2) * 5.0 * D4 * exp2; + temp -= (5.0 - 2.0 * D2) * D6 * 2.0 * exp2; + pDetail->ddZdD += pDetail->fx[45] * temp; + pDetail->ddZdD += pDetail->fx[46] * temp; + temp = -16.0 * D8 * exp4 + (5.0 - 4.0 * D4) * 5.0 * D4 * exp4; + temp -= (5.0 - 4.0 * D4) * D8 * 4.0 * exp4; + pDetail->ddZdD += pDetail->fx[47] * temp; + pDetail->ddZdD += pDetail->fx[48] * temp; + temp = 36.0 * D5; + pDetail->ddZdD += pDetail->fx[49] * temp; + temp = -4.0 * D7 * exp2 + (6.0 - 2.0 * D2) * 6.0 * D5 * exp2; + temp -= (6.0 - 2.0 * D2) * D7 * 2.0 * exp2; + pDetail->ddZdD += pDetail->fx[50] * temp; + temp = 49.0 * D6; + pDetail->ddZdD += pDetail->fx[51] * temp; + temp = -4.0 * D8 * exp2 + (7.0 - 2.0 * D2) * 7.0 * D6 * exp2; + temp -= (7.0 - 2.0 * D2) * D8 * 2.0 * exp2; + pDetail->ddZdD += pDetail->fx[52] * temp; + temp = -1.0 * D8 * exp1 + (8.0 - D1) * 8.0 * D7 * exp1; + temp -= (8.0 - D1) * D8 * exp1; + pDetail->ddZdD += pDetail->fx[53] * temp; + temp = -4.0 * D1 * D8 * exp2 + (8.0 - 2.0 * D2) * 8.0 * D7 * exp2; + temp -= (8.0 - 2.0 * D2) * D8 * 2.0 * D1 * exp2; + pDetail->ddZdD += pDetail->fx[54] * temp; + pDetail->ddZdD += pDetail->fx[55] * temp; + temp = -4.0 * D2 * D8 * exp2 + (9.0 - 2.0 * D2) * 9.0 * D8 * exp2; + temp -= (9.0 - 2.0 * D2) * D2 * D8 * 2.0 * exp2; + pDetail->ddZdD += pDetail->fx[56] * temp; + pDetail->ddZdD += pDetail->fx[57] * temp; + pDetail->ddZdD *= pDetail->dKp3; + return pDetail->ddZdD; +} + +void Detail_relativedensity(Detail *pDetail, NGParSTRUCT *ptNGPar) { + double dBX, dZa; + const double dMWair = 28.96256; + 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); +} diff --git a/NG/Detail.h b/NG/Detail.h new file mode 100644 index 0000000..bd0cc27 --- /dev/null +++ b/NG/Detail.h @@ -0,0 +1,122 @@ +/************************************************************************* +* �ļ�: detail.h +**************************************************************************/ + +#ifndef _DETAIL_H +#define _DETAIL_H + +#include "NGCal.h" +#include + + +typedef struct Detail { + + int iNCC; + + int aiCID[21]; + + double dOldMixID; + double dOldPb; + double dOldTb; + double dOldPf; + double dOldTf; + + + double adAn[58]; + double adUn[58]; + + + double dMri[21]; + double dEi[21]; + double dKi[21]; + double dGi[21]; + double dQi[21]; + double dFi[21]; + double dSi[21]; + double dWi[21]; + + double dEij[21][21]; + double dUij[21][21]; + double dKij[21][21]; + double dGij[21][21]; + + double adTable6Eij[21][21]; + double adTable6Uij[21][21]; + double adTable6Kij[21][21]; + double adTable6Gij[21][21]; + + double adTable5Qi[21]; + double adTable5Fi[21]; + double adTable5Si[21]; + double adTable5Wi[21]; + + double adTableHhvMol[4][21]; + double adTableLhvMol[4][21]; + + double dXi[21]; + double dPCalc; + double dT; + double dP; + double dRhoTP; + double dB; + double adBcoef[18]; + double adFn[58]; + double fx[58]; + double dU; + double dKp3; + double dW; + double dQp2; + double dF; + double dRho; + double dRhoL; + double dRhoH; + double dPRhoL; + double dPRhoH; + + + double dZ; + double ddZdT; + double dd2ZdT2; + double ddZdD; + double ddBdT; + double dd2BdT2; +} Detail; + + +Detail *Detail_Construct(void); + +void Detail_Destroy(Detail *pDetail); + + +int Detail_compositionchange(Detail *pDetail, NGParSTRUCT *pAGA10); + +int Detail_table(Detail *pDetail); + +void Detail_paramdl(Detail *pDetail); + +void Detail_chardl(Detail *pDetail, NGParSTRUCT *pAGA10); + +void Detail_bvir(Detail *pDetail); + +void Detail_temp(Detail *pDetail); + +void Detail_braket(Detail *pDetail, NGParSTRUCT *pAGA10); + +void Detail_pdetail(Detail *pDetail, double dRho); + +void Detail_ddetail(Detail *pDetail, NGParSTRUCT *pAGA10); + +void Detail_relativedensity(Detail *pDetail, NGParSTRUCT *pAGA10); + + +double Detail_zdetail(Detail *pDetail, double dRho); + +double Detail_dZdT(Detail *pDetail, double dRho); + +double Detail_d2ZdT2(Detail *pDetail, double dRho); + +double Detail_dZdD(Detail *pDetail, double dRho); + +void Detail_Run(Detail *pDetail, NGParSTRUCT *ptNGPar); + +#endif diff --git a/NG/FlowCal.c b/NG/FlowCal.c new file mode 100644 index 0000000..2bfebbd --- /dev/null +++ b/NG/FlowCal.c @@ -0,0 +1,427 @@ +// +// Created by ldeyu on 2025/7/7. +// +#include "NGCal.h" +#include "FlowCal.h" + +void OFlowCal(FlowParSTRUCT *ptFlowPar, NGParSTRUCT *ptNGPar) { + + double tempPatm = ptFlowPar->dPatm * 1000000; + + double tempPf = ptFlowPar->dPf * 1000000; + + + double tempDP = ptFlowPar->dDp * 1000; + + double tempTf = ptFlowPar->dTf + 273.15; + + if (ptFlowPar->dPfType == 0) + { + ptFlowPar->dPf = tempPatm + tempPf; + ptNGPar->dPf = tempPatm + tempPf; + } else { + ptFlowPar->dPf = tempPf; + ptNGPar->dPf = tempPf; + } + + ptFlowPar->dDp = tempDP; + ptFlowPar->dTf = tempTf; + ptNGPar->dTf = tempTf; + ptNGPar->dCbtj = ptFlowPar->dCbtj; + + switch (ptNGPar->dCbtj) { + case 2: + ptNGPar->dPb = 101325; + ptNGPar->dTb = 273.15; + ptFlowPar->dPb_M = (101325); + ptFlowPar->dTb_M = (273.15); + break; + case 1: + + ptNGPar->dPb = (101325); + ptNGPar->dTb = (288.15); + ptFlowPar->dPb_M = (101325); + ptFlowPar->dTb_M = (288.15); + break; + case 0: + + ptNGPar->dPb = (101325); + ptNGPar->dTb = (293.15); + ptFlowPar->dPb_M = (101325); + ptFlowPar->dTb_M = (293.15); + break; + } + + double ngArray[NUMBEROFCOMPONENTS]; + + for (int i = 0; i < NUMBEROFCOMPONENTS; i++) { + ngArray[i] = ptFlowPar->dNG_Compents[i] / 100; + ptNGPar->adMixture[i] = ngArray[i]; + } + Crit(ptNGPar, 0); + + + ptFlowPar->dFpv = ptNGPar->dFpv; + + + + 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 = calculateE(ptFlowPar->dBeta); + ptFlowPar->dFG = calculateFG(ptNGPar->dRD_Real); + ptFlowPar->dFT = calculateFT(ptFlowPar->dTb_M, ptFlowPar->dTf); + + 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; + double beta = ptFlowPar->dBeta; + double P1 = ptFlowPar->dPf; + double deltaP = ptFlowPar->dDp; + double Tf = ptFlowPar->dTf; + + + double C_initial = 0.6; + 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 { + 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"); + } + } 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; + + 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->dVelocityFlow = ptFlowPar->dVFlowf / (M_PI * pow((ptFlowPar->dPipeD / 2000), 2)); +} + + + + + + +double CaiLiaoPzxs(int tempCaiLiao) { + double CaiLiaoPzxs = 0; + switch (tempCaiLiao) { + case 0: + CaiLiaoPzxs = 11.75; + break; + + case 1: + CaiLiaoPzxs = 11.6; + break; + + case 2: + CaiLiaoPzxs = 11.16; + break; + + case 3: + CaiLiaoPzxs = 11.59; + break; + + case 4: + CaiLiaoPzxs = 10.5; + break; + + case 5: + CaiLiaoPzxs = 10.0; + break; + + case 6: + CaiLiaoPzxs = 10.2; + break; + + case 7: + CaiLiaoPzxs = 15.5; + break; + + case 8: + CaiLiaoPzxs = 11.5; + break; + + case 9: + CaiLiaoPzxs = 10.8; + break; + + case 10: + CaiLiaoPzxs = 16.6; + break; + + case 11: + CaiLiaoPzxs = 11.4; + break; + + case 12: + CaiLiaoPzxs = 16.55; + break; + + case 13: + CaiLiaoPzxs = 17.8; + break; + + case 14: + CaiLiaoPzxs = 17.2; + break; + } + return CaiLiaoPzxs; +} + + +double calculateK(int dPipeType) { + double Jdccd; + switch (dPipeType) { + case 0: + Jdccd = 0.029; + break; + case 1: + case 2: + case 3: + Jdccd = 0.075; + break; + case 4: + Jdccd = 0.1; + break; + case 5: + Jdccd = 0.15; + break; + case 6: + Jdccd = 1; + break; + case 7: + Jdccd = 2.1; + break; + case 8: + Jdccd = 0.04; + break; + case 9: + Jdccd = 0.15; + break; + case 10: + Jdccd = 0.13; + break; + case 11: + Jdccd = 0.25; + break; + default: + + fprintf(stderr, "δ֪Ĺܵ: %d\n", dPipeType); + return FLOW_CALC_ERROR; + } + return Jdccd; +} + + + +double calculateRoughnessFactor(double D_pipe, double K, double C) { + + double K_over_D = K / D_pipe; + + + if (K_over_D <= 0.0004) { + return 1.0000; + } + + + double term = (K_over_D * 1e6) - 400; + if (term < 0) { + fprintf(stderr, "K/D <20><>ʽ÷Χ\n"); + return FLOW_CALC_ERROR; + } + + double G_me = 1 + (0.011 / C) * sqrt(term); + return G_me; +} + + + +double calculateE(double beta) { + return 1 / sqrt(1 - pow(beta, 4)); +} + + +double calculateFG(double dRD_Real) { + return 1 / sqrt(dRD_Real); +} + + +double calculateFT(double dTb_M, double dTf) { + return sqrt(dTb_M / dTf); +} + + +double calculateEpsilon(double dPf, double dDp, double beta, double dKappa) { + double tau = (dPf - dDp) / dPf; + double epsilon = 1 - (0.351 + 0.256 * pow(beta, 4) + 0.93 * pow(beta, 8)) * (1 - pow(tau, 1 / dKappa)); + return epsilon; +} + + +double calculateKappa(double dZf) { + + double gamma = 1.3; + double Z = dZf; + + + double kappa = gamma / (1 - (gamma - 1) * (1 / Z - 1)); + return kappa; +} + + +double calculateReD(double Qf, double D, double rho, double mu) { + return (4 * Qf * rho) / (M_PI * D * mu); +} + + +double calculateCd(double beta, double ReD, double D_mm, int ptMode) { + double L1, L2; + + switch (ptMode) { + case 1: + L1 = L2 = 0; + break; + case 0: + L1 = L2 = 25.4 / D_mm; + break; + case 2: + L1 = 1.0; + L2 = 0.47; + break; + default: + fprintf(stderr, "ֵ֧ȡѹʽ: %d\n", ptMode); + return FLOW_CALC_ERROR; + } + + double term1 = 0.5961 + 0.0261 * pow(beta, 2) - 0.216 * pow(beta, 8); + double term2 = 0.000521 * pow(1e6 * beta / ReD, 0.7); + double A = pow(19000 * beta / ReD, 0.8); + double term3 = (0.0188 + 0.0063 * A) * pow(beta, 3.5) * pow(1e6 / ReD, 0.3); + double term4 = (0.043 + 0.08 * exp(-10 * L1) - 0.123 * exp(-7 * L1)) + * (1 - 0.11 * A) * pow(beta, 4) / (1 - pow(beta, 4)); + double term5 = -0.031 * (2 * L2 / (1 - beta) - 0.8 * pow(2 * L2 / (1 - beta), 1.1)) + * pow(beta, 1.3); + + double Cd = term1 + term2 + term3 + term4 + term5; + + + if (D_mm < 71.12) { + Cd += 0.011 * (0.75 - beta) * (2.8 - D_mm / 25.4); + } + return Cd; +} + + +double Dlndjs(double tempP_jy, double tempT) { + double Dlndjs_Dlnd_Data[8][11] = { + {976, 991, 1014, 1044, 1073, 1114, 1156, 1207, 1261, 1331, 1405}, + {1027, 1040, 1063, 1091, 1118, 1151, 1185, 1230, 1276, 1331, 1389}, + {1071, 1082, 1106, 1127, 1149, 1180, 1211, 1250, 1289, 1335, 1383}, + {1123, 1135, 1153, 1174, 1195, 1224, 1253, 1289, 1324, 1366, 1409}, + {1167, 1178, 1196, 1216, 1236, 1261, 1287, 1318, 1350, 1385, 1421}, + {1213, 1224, 1239, 1257, 1275, 1297, 1320, 1346, 1373, 1403, 1435}, + {1260, 1270, 1281, 1297, 1313, 1333, 1352, 1374, 1396, 1424, 1451}, + {1303, 1312, 1323, 1338, 1352, 1372, 1391, 1412, 1432, 1456, 1482} + }; + + double Dlndjs_Dlnd_T[8] = { + -15 + 273.15, 0 + 273.15, 15 + 273.15, 30 + 273.15, + 45 + 273.15, 60 + 273.15, 75 + 273.15, 90 + 273.15 + }; + + double Dlndjs_Dlnd_P[11] = {0.1, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10}; + + double s1, s2, ky, kx; + int i, m = 0, n = 0; + + if (tempT < Dlndjs_Dlnd_T[0]) { + tempT = Dlndjs_Dlnd_T[0]; + } + if (tempT > Dlndjs_Dlnd_T[7]) { + tempT = Dlndjs_Dlnd_T[7]; + } + if (tempP_jy < Dlndjs_Dlnd_P[0]) { + tempP_jy = Dlndjs_Dlnd_P[0]; + } + if (tempP_jy > Dlndjs_Dlnd_P[10]) { + tempP_jy = Dlndjs_Dlnd_P[10]; + } + + for (i = 0; i <= 6; i++) { + if (tempT >= Dlndjs_Dlnd_T[i] && tempT <= Dlndjs_Dlnd_T[i + 1]) { + m = i; + break; + } + } + + for (i = 0; i <= 9; i++) { + if (tempP_jy >= Dlndjs_Dlnd_P[i] && tempP_jy <= Dlndjs_Dlnd_P[i + 1]) { + n = i; + break; + } + } + + if (Dlndjs_Dlnd_P[n + 1] - Dlndjs_Dlnd_P[n] != 0) { + ky = (tempP_jy - Dlndjs_Dlnd_P[n]) / (Dlndjs_Dlnd_P[n + 1] - Dlndjs_Dlnd_P[n]); + } else { + ky = 0; + } + + if (Dlndjs_Dlnd_T[m + 1] - Dlndjs_Dlnd_T[m] != 0) { + kx = (tempT - Dlndjs_Dlnd_T[m]) / (Dlndjs_Dlnd_T[m + 1] - Dlndjs_Dlnd_T[m]); + } else { + kx = 0; + } + + 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; +} diff --git a/NG/FlowCal.h b/NG/FlowCal.h new file mode 100644 index 0000000..d47e2e0 --- /dev/null +++ b/NG/FlowCal.h @@ -0,0 +1,87 @@ +#ifndef FLOWCAL_H +#define FLOWCAL_H +#include "NGCal.h" +typedef struct FlowParSTRUCT { + + + int dCbtj; + double dPb_M; + double dTb_M; + double dPb_E; + double dTb_E; + double dPatm; + double dNG_Compents[21]; + + int dMeterType; + int dCoreType; + int dPtmode; + int dPipeType; + double dPipeD; + int dPipeMaterial; + + double dOrificeD; + int dOrificeMaterial; + + double dPf; + int dPfType; + double dTf; + double dDp; + + double dMeterFactor; + double dPulseNum; + + + + double dE; + double dFG; + double dFT; + double dDViscosity; + double dDExpCoefficient; + double dRnPipe; + double dBk; + double dRoughNessPipe; + double dCd; + double dCdCorrect; + double dCdNozell; + double dVFlowb; + double dVFlowf; + double dMFlowb; + double dEFlowb; + double dVelocityFlow; + double dPressLost; + double dBeta; + double dKappa; + double dFpv; +} FlowParSTRUCT; + +double CaiLiaoPzxs(int tempCaiLiao); + +double calculateK(int dPipeType); + +double calculateRoughnessFactor(double D_pipe, double K, double C); + +void thermalExpansionCorrection(double dOrificeMaterial, double dOrificeD, + double dPipeMaterial, double dPipeD, + double dTf, double correctedValues[3]); + +double calculateE(double beta); + +double calculateFG(double dRD_Real); + +double calculateFT(double dTb_M, double dTf); + +double calculateEpsilon(double dPf, double dDp, double beta, double dKappa); + +double calculateKappa(double dZf); + +double calculateReD(double Qf, double D, double rho, double mu); + +double calculateCd(double beta, double ReD, double D_mm, int ptMode); + + +double Dlndjs(double tempP_jy, double tempT); + +void OFlowCal(FlowParSTRUCT *ptFlowPar, NGParSTRUCT *ptNGPar); + +#endif + diff --git a/NG/main.c b/NG/main.c new file mode 100644 index 0000000..086efa6 --- /dev/null +++ b/NG/main.c @@ -0,0 +1,179 @@ +#include +#include "NGCal.h" +#include "FlowCal.h" + +int main() { + // 定义并初始化 FlowParSTRUCT 结构体变量 + FlowParSTRUCT flowParams = {0}; + NGParSTRUCT ngParams = {0}; + + // 设置基本参数 + flowParams.dPatm = 0.0981; // 标准大气压(bar) + flowParams.dPf = 4; // 压力(MPa) + flowParams.dPfType = 1; // 0=表压,1=绝压 + flowParams.dDp = 12.50; // 差压(kPa) + flowParams.dTf = 10; // 温度(°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] = 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) + // flowParams.dNG_Compents[3] =1.80; // 乙烷(C2H6) + // flowParams.dNG_Compents[4] =0.45; // 丙烷(C3H8) + // flowParams.dNG_Compents[5] =0; // 水(H2O) + // flowParams.dNG_Compents[6] =0; // 硫化氢(H2S) + // flowParams.dNG_Compents[7] =0; // 氢气(H2) + // flowParams.dNG_Compents[8] =0; // 一氧化碳(CO) + // flowParams.dNG_Compents[9] =0; // 氧气(O2) + // flowParams.dNG_Compents[10]= 0.1; // 异丁烷(i-C4H10) + // flowParams.dNG_Compents[11]= 0.1; // 正丁烷(n-C4H10) + // flowParams.dNG_Compents[12]= 0.05; // 异戊烷(i-C5H12) + // flowParams.dNG_Compents[13]= 0.03; // 正戊烷(n-C5H12) + // flowParams.dNG_Compents[14]= 0.07; // 己烷(C6H14) + // flowParams.dNG_Compents[15]= 0; // 庚烷(C7H16) + // flowParams.dNG_Compents[16]= 0; // 辛烷(C8H18) + // flowParams.dNG_Compents[17]= 0; // 壬烷(C9H20) + // flowParams.dNG_Compents[18]= 0; // 癸烷(C10H22) + // flowParams.dNG_Compents[19]= 0; // 氦气(He) + // flowParams.dNG_Compents[20]= 0; // 其他组分 + + // // 显式调用 NGCal_Init 初始化模块 + // if (NGCal_NGCal != NGCal_Init()) { + // printf("错误:NGCal 初始化失败!\n"); + // return -1; // 退出程序 + // } + + // 调用流量计算函数 + OFlowCal(&flowParams, &ngParams); + + // 打印计算结果 + printf("工况条件信息:\n"); + printf("标准参比条件: %d\n", flowParams.dCbtj); + printf("计量参比压力: %.2f\n", flowParams.dPb_M); + printf("计量参比温度: %.2f\n", flowParams.dTb_M); + printf("能量参比压力: %.2f\n", flowParams.dPb_E); + printf("能量参比温度: %.2f\n", flowParams.dTb_E); + printf("大气压力: %.2f Pa\n", flowParams.dPatm); + printf("天然气组分:\n"); + for (int i = 0; i < 21; i++) { + printf(" 组分 %d: %.6f\n", i, flowParams.dNG_Compents[i]); + } + + printf("\n仪表参数:\n"); + printf("仪表类型: %d\n", flowParams.dMeterType); + printf("核心类型: %d\n", flowParams.dCoreType); + printf("取压方式: %d\n", flowParams.dPtmode); + printf("管道类型: %d\n", flowParams.dPipeType); + printf("管道内径: %.2f mm\n", flowParams.dPipeD); + printf("管道材质: %d\n", flowParams.dPipeMaterial); + printf("孔板直径: %.2f mm\n", flowParams.dOrificeD); + printf("孔板材质: %d\n", flowParams.dOrificeMaterial); + + printf("\n测量值:\n"); + printf("压力: %.2f Pa\n", flowParams.dPf); + printf("压力类型: %d\n", flowParams.dPfType); + printf("温度: %.2f K\n", flowParams.dTf); + printf("差压: %.2f Pa\n", flowParams.dDp); + printf("仪表系数: %.6f\n", flowParams.dMeterFactor); + printf("脉冲数: %.2f\n", flowParams.dPulseNum); + + printf("\n计算结果:\n"); + printf("膨胀系数: %.6f\n", flowParams.dE); + printf("相对密度系数: %.6f\n", flowParams.dFG); + printf("超压缩系数: %.6f\n", flowParams.dFT); + printf("动力粘度: %.6f\n", flowParams.dDViscosity); + printf("热膨胀系数: %.6f\n", flowParams.dDExpCoefficient); + printf("管道雷诺数: %.2f\n", flowParams.dRnPipe); + printf("孔板弯曲系数: %.6f\n", flowParams.dBk); + printf("管道粗糙度: %.6f\n", flowParams.dRoughNessPipe); + printf("流出系数: %.6f\n", flowParams.dCd); + printf("流出系数修正: %.6f\n", flowParams.dCdCorrect); + printf("喷嘴流出系数: %.6f\n", flowParams.dCdNozell); + printf("标况体积流量: %.6f Nm3/s\n", flowParams.dVFlowb); + printf("工况体积流量: %.6f m3/s\n", flowParams.dVFlowf); + printf("质量流量: %.6f t/s\n", flowParams.dMFlowb); + printf("能量流量: %.6f MJ/s\n", flowParams.dEFlowb); + printf("流速: %.6f m/s\n", flowParams.dVelocityFlow); + printf("压力损失: %.6f\n", flowParams.dPressLost); + printf("直径比: %.6f\n", flowParams.dBeta); + printf("等熵指数: %.6f\n", flowParams.dKappa); + printf("压缩因子: %.6f\n", flowParams.dFpv); + + printf("状态: %ld\n", ngParams.lStatus); + printf("强制更新标志: %d\n", ngParams.bForceUpdate); + printf("混合比:\n"); + for (int i = 0; i < 21; i++) { + printf(" 组分 %d: %.6f\n", i, ngParams.adMixture[i]); + } + printf("参比条件: %d\n", ngParams.dCbtj); + printf("标准压力: %.2f Pa\n", ngParams.dPb); + printf("标准温度: %.2f K\n", ngParams.dTb); + printf("工作压力: %.2f Pa\n", ngParams.dPf); + printf("工作温度: %.2f K\n", ngParams.dTf); + + printf("\nAGA 8 详细计算结果:\n"); + printf("平均分子量: %.6f\n", ngParams.dMrx); + printf("标准条件下压缩因子: %.6f\n", ngParams.dZb); + printf("工作条件下压缩因子: %.6f\n", ngParams.dZf); + printf("超压缩因子: %.6f\n", ngParams.dFpv); + printf("标准条件下摩尔密度: %.6f moles/dm3\n", ngParams.dDb); + printf("工作条件下摩尔密度: %.6f moles/dm3\n", ngParams.dDf); + printf("标准条件下密度: %.6f kg/m3\n", ngParams.dRhob); + printf("工作条件下密度: %.6f kg/m3\n", ngParams.dRhof); + printf("理想相对密度: %.6f\n", ngParams.dRD_Ideal); + printf("实际相对密度: %.6f\n", ngParams.dRD_Real); + + printf("\n热力学性质:\n"); + printf("理想焓: %.6f\n", ngParams.dHo); + printf("实际焓: %.6f J/kg\n", ngParams.dH); + printf("实际熵: %.6f J/kg-mol.K\n", ngParams.dS); + printf("理想定压比热: %.6f J/kg-mol.K\n", ngParams.dCpi); + printf("实际定压比热: %.6f J/kg-mol.K\n", ngParams.dCp); + printf("实际定容比热: %.6f J/kg-mol.K\n", ngParams.dCv); + printf("比热比: %.6f\n", ngParams.dk); + printf("等熵指数: %.6f\n", ngParams.dKappa); + printf("声速: %.6f m/s\n", ngParams.dSOS); + printf("临界流函数: %.6f\n", ngParams.dCstar); + + printf("\n单位摩尔高热值: %.6f\n", ngParams.dHhvMol); + printf("单位摩尔低热值: %.6f\n", ngParams.dLhvMol); +}