commit 2fe122d9cb4ef74dc4bd4221f6bffc59b879f52f Author: 廖德云 Date: Sat Jul 12 17:09:07 2025 +0800 天然气流量计算，压缩因子，声速计算的C语言版本 diff --git a/.idea/.gitignore b/.idea/.gitignore new file mode 100644 index 0000000..35410ca --- /dev/null +++ b/.idea/.gitignore @@ -0,0 +1,8 @@ +# 默认忽略的文件 +/shelf/ +/workspace.xml +# 基于编辑器的 HTTP 客户端请求 +/httpRequests/ +# Datasource local storage ignored files +/dataSources/ +/dataSources.local.xml diff --git a/.idea/NG.iml b/.idea/NG.iml new file mode 100644 index 0000000..f08604b --- /dev/null +++ b/.idea/NG.iml @@ -0,0 +1,2 @@ + + \ No newline at end of file diff --git a/.idea/editor.xml b/.idea/editor.xml new file mode 100644 index 0000000..25c6c37 --- /dev/null +++ b/.idea/editor.xml @@ -0,0 +1,344 @@ + + + + + \ No newline at end of file diff --git a/.idea/encodings.xml b/.idea/encodings.xml new file mode 100644 index 0000000..7605d94 --- /dev/null +++ b/.idea/encodings.xml @@ -0,0 +1,6 @@ + + + + + + \ No newline at end of file diff --git a/.idea/misc.xml b/.idea/misc.xml new file mode 100644 index 0000000..0b76fe5 --- /dev/null +++ b/.idea/misc.xml @@ -0,0 +1,7 @@ + + + + + + \ No newline at end of file diff --git a/.idea/modules.xml b/.idea/modules.xml new file mode 100644 index 0000000..1ad4cbd --- /dev/null +++ b/.idea/modules.xml @@ -0,0 +1,8 @@ + + + + + + + + \ No newline at end of file diff --git a/.idea/vcs.xml b/.idea/vcs.xml new file mode 100644 index 0000000..b2bdec2 --- /dev/null +++ b/.idea/vcs.xml @@ -0,0 +1,6 @@ + + + + + + \ No newline at end of file diff --git a/CMakeLists.txt b/CMakeLists.txt new file mode 100644 index 0000000..5b5bb14 --- /dev/null +++ b/CMakeLists.txt @@ -0,0 +1,15 @@ +cmake_minimum_required(VERSION 3.31) +project(NG C) + +set(CMAKE_C_STANDARD 99) + +add_executable(NG + main.c + Detail.c + FlowCal.h # 头文件明确列出 + main.c + NGCal.c + NGCal.h + # 确保包含实现文件 + Therm.c + FlowCal.c) diff --git a/Detail.c b/Detail.c new file mode 100644 index 0000000..4645113 --- /dev/null +++ b/Detail.c @@ -0,0 +1,1133 @@ +#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, const NGParSTRUCT *ptNGPar) { + double dMixID = 0.0; + for (int i = 0; i < NUMBEROFCOMPONENTS; i++) + dMixID += ((i + 2) * ptNGPar->adMixture[i]); + if (dMixID != pDetail->dOldMixID) { + pDetail->dOldMixID = dMixID; + return 1; + } + return 0; +} + +void Detail_Run(Detail *pDetail, NGParSTRUCT *ptNGPar) { + const int bCompChange = Detail_compositionchange(pDetail, ptNGPar); + ptNGPar->bForceUpdate = ptNGPar->bForceUpdate || bCompChange; + if (ptNGPar->bForceUpdate) { + pDetail->iNCC = 0; + for (int 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 = 1; + } + 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; + } + 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 = 0; + } +} + +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; + + + + return 0; +} + +void Detail_paramdl(Detail *pDetail) { + 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--) { + 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]]; + } + 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_dhvMol(Detail *pDetail, NGParSTRUCT *ptNGPar) { + const double adTableHhvMol[4][NUMBEROFCOMPONENTS] = { + { //0 + 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 + }, + { //15 + 891.56, 0, 0, 1562.14, 2221.1, 44.43, 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 + }, + { //20 + 891.09, 0, 0, 1561.41, 2220.13, 44.22, 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 + }, + { //25 + 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 + } + }; + + const double adTableLhvMol[4][NUMBEROFCOMPONENTS] = { + { //0 + 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 + }, + { //15 + 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 + }, + { //20 + 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 + }, + { //25 + 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 + } + }; + for (int i = 0; i dCbtj) + { + case 0: + ptNGPar->dHhvMol += adTableHhvMol[0][i] * ptNGPar->adMixture[i]; + ptNGPar->dLhvMol +=adTableLhvMol[0][i] * ptNGPar->adMixture[i]; + break; + case 1: + ptNGPar->dHhvMol += adTableHhvMol[1][i] * ptNGPar->adMixture[i]; + ptNGPar->dLhvMol += adTableLhvMol[1][i] * ptNGPar->adMixture[i]; + break; + case 2: + ptNGPar->dHhvMol += adTableHhvMol[2][i] * ptNGPar->adMixture[i]; + ptNGPar->dLhvMol += adTableLhvMol[2][i] * ptNGPar->adMixture[i]; + break; + default: ; + } + } +} + +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]; + + 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++) { + const 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); + 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); + 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; + 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); + 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; + } + } + + + 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; + 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 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); + 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; + + 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) { + double xnumer, xdenom; + const int imax = 150; + 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; + 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; + } + const double delmin = epsmin * fabs(x2); + const double 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; + 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) { + const double sgndel = delbis / fabs(delbis); + delx = 1.0000009 * sgndel * delmin; + delprv = delx; + } + const double 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) { + double rho2; + const int imax = 200; + double rho1 = 0.0; + double p1 = 0.0; + double rhomax = 1.0 / pDetail->dKp3; + if (pDetail->dT > 1.2593 * pDetail->dU) + rhomax = 20.0 * rhomax; + const double 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++) { + 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); + const double 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; +} + +void Detail_pdetail(Detail *pDetail, double dD) { + pDetail->dPCalc = Detail_zdetail(pDetail, dD) * dD * RGASKJ * pDetail->dT; +} + +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); + 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) { + 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++) { + 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); + double 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) { + 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++) { + 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); + double 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; + 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++) { + 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(const Detail *pDetail, NGParSTRUCT *ptNGPar) { + 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); + ptNGPar->dRD_Ideal = ptNGPar->dMrx / dMWair; + ptNGPar->dRD_Real = ptNGPar->dRD_Ideal * (dZa / ptNGPar->dZb); +} diff --git a/Detail.h b/Detail.h new file mode 100644 index 0000000..f350318 --- /dev/null +++ b/Detail.h @@ -0,0 +1,120 @@ +/************************************************************************* +* �ļ�: detail.h +**************************************************************************/ + +#ifndef _DETAIL_H +#define _DETAIL_H + +#include "NGCal.h" + + +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 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, const NGParSTRUCT *pAGA10); + +int Detail_table(Detail *pDetail); + +void Detail_paramdl(Detail *pDetail); + +void Detail_chardl(Detail *pDetail, NGParSTRUCT *pAGA10); + +void Detail_dhvMol(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(const 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/FlowCal.c b/FlowCal.c new file mode 100644 index 0000000..ebed82c --- /dev/null +++ b/FlowCal.c @@ -0,0 +1,501 @@ +// +// Created by ldeyu on 2025/7/7. +// +#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]; + + // 生成动态格式字符串（如 "%.4f"） + snprintf(format_str, sizeof(format_str), "%%.%df", digits); + + // 格式化数值到字符串 + snprintf(buffer, sizeof(buffer), format_str, value); + + // 转换回 double + double result = strtod(buffer, NULL); + return result; +} + +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; + default: ; + } + + 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 = format_double(ptNGPar->dFpv, 4); + + + 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->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->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); + + 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; + 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"); + } + } 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->dPb_M * 1e-6 / RGASKJ / ptFlowPar->dTb_M; + ptFlowPar->dVelocityFlow = ptFlowPar->dVFlowf / (M_PI * pow((ptFlowPar->dPipeD / 2000), 2)); +} + + +double CaiLiaoPzxs(const 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; + default: ; + } + 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 * 1000) / (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 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; + } + + 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; + return (s1 + (s2 - s1) * kx) / 100000.0; +} + +// 压力损失计算 +double YaLiSunShi(double tempLiuChuXiShu, double tempZjb, double tempDp, int JieLiuZhuangZhi) { + double ylss = 0; + switch (JieLiuZhuangZhi) { + case 0: + case 1: + case 2: + ylss = tempDp * (sqrt(1 - tempZjb) - tempLiuChuXiShu * pow(tempZjb, 2)) + / (sqrt(1 - tempZjb) + tempLiuChuXiShu * pow(tempZjb, 2)); + break; + default: ; + } + return ylss; +} + + +// 标况转工况流量转换 +double FlowConvert_BaseToWork(double Pf, double Tf, double Zb, double Zf, double FlowBase, int Cbtj) { + double tempPn = 0; + double tempTn = 0; + + + switch (Cbtj) { + case 2: + tempPn = 101325; + tempTn = 273.15; + break; + case 1: + tempPn = 101325; + tempTn = 288.15; + break; + case 0: + tempPn = 101325; + tempTn = 293.15; + break; + case 3: + tempPn = 10155981; + tempTn = 288.7055555; + break; + default: ; + } + return FlowBase * tempPn * Tf * Zf / Pf / tempTn / Zb; +} + +// 工况转标况流量转换 +double FlowConvert_WorkToBase(double Pf, double Tf, double Zb, double Zf, double FlowWork, int Cbtj) { + double tempPn = 0; + double tempTn = 0; + + switch (Cbtj) { + case 2: + tempPn = 101325; + tempTn = 273.15; + break; + case 1: + tempPn = 101325; + tempTn = 288.15; + break; + case 0: + tempPn = 101325; + tempTn = 293.15; + break; + case 3: + tempPn = 10155981; + tempTn = 288.7055555; + break; + default: ; + } + return FlowWork * Pf * tempTn * Zb / tempPn / Tf / Zf; +} diff --git a/FlowCal.h b/FlowCal.h new file mode 100644 index 0000000..23cf7ea --- /dev/null +++ b/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/NGCal.c b/NGCal.c new file mode 100644 index 0000000..a9d8e1b --- /dev/null +++ b/NGCal.c @@ -0,0 +1,86 @@ +#include "NGCal.h" +#include "Therm.h" +#include "Detail.h" +#include "math.h" +static Therm *ptTherm; +static Detail *ptDetail; + +int NGCal_Init(NGParSTRUCT *ptNGPar) { + + ptDetail = Detail_Construct(); + if (NULL == ptDetail) { + return MEMORY_ALLOCATION_ERROR; + } + + + ptTherm = (Therm *) malloc(sizeof(Therm)); + Therm_Init(ptTherm); + if (NULL == ptTherm) { + return MEMORY_ALLOCATION_ERROR; + } + + return NGCal_NGCal; +} + +int NGCal_UnInit(void) { + if (ptDetail) free(ptDetail); + if (ptTherm) free(ptTherm); + return 0; +} + +double SOS(NGParSTRUCT *ptNGPar) { + if (NULL == ptDetail || NULL == ptTherm) { + NGCal_UnInit(); + NGCal_Init(ptNGPar); + } + Therm_Run(ptTherm, ptNGPar, ptDetail); + ptNGPar->dCstar = 0.0; + return ptNGPar->dSOS; +} + +double Crit(NGParSTRUCT *ptNGPar, double dPlenumVelocity) { + if (NULL == ptDetail || NULL == ptTherm) { + NGCal_UnInit(); + if (NGCal_NGCal != NGCal_Init(ptNGPar)) { + ptNGPar->lStatus = MEMORY_ALLOCATION_ERROR; + return 0.0; + } + } + + + 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); + Therm_Run(ptTherm, ptNGPar, ptDetail); + double DDH = DH; + DH = (ptNGPar->dSOS * ptNGPar->dSOS - dPlenumVelocity * dPlenumVelocity) / 2.0; + if (fabs(DDH - DH) < tolerance) break; + } + ptNGPar->dCstar = (ptNGPar->dRhof * ptNGPar->dSOS) / sqrt(R * P * Z); + ptNGPar->dPf = P; + ptNGPar->dTf = T; + Therm_Run(ptTherm, ptNGPar, ptDetail); + Detail_dhvMol(ptDetail,ptNGPar); + 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); + return (sqrt(k * pow(root, exponent))); +} + +double CRi(const NGParSTRUCT *ptNGPar) { + return (Cperf(ptNGPar) / sqrt(ptNGPar->dZf)); +} diff --git a/NGCal.h b/NGCal.h new file mode 100644 index 0000000..bfa0869 --- /dev/null +++ b/NGCal.h @@ -0,0 +1,80 @@ +/************************************************************************* +* �ļ�: NGCal.h +**************************************************************************/ +#ifndef _NGCal_H +#define _NGCal_H +#include +#include +#define NORMAL 9000 +#define NGCal_NGCal 9001 +#define MEMORY_ALLOCATION_ERROR 9002 +#define GENERAL_CALCULATION_FAILURE 9003 +#define MAX_NUM_OF_ITERATIONS_EXCEEDED 9004 +#define NEGATIVE_DENSITY_DERIVATIVE 9005 +#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; +} NGParSTRUCT; + + +enum gascomp { + XiC1=0, XiN2, XiCO2, XiC2, XiC3, + XiH2O, XiH2S, XiH2, XiCO, XiO2, + XiIC4, XiNC4, XiIC5, XiNC5, XiNC6, + XiNC7, XiNC8, XiNC9, XiNC10, XiHe, XiAr +}; + + + int NGCal_Init(NGParSTRUCT * ptNGPar); + int NGCal_UnInit(void); + + + double SOS(NGParSTRUCT *); + + + double Crit(NGParSTRUCT *, double); + +#endif diff --git a/Therm.c b/Therm.c new file mode 100644 index 0000000..52a0c0d --- /dev/null +++ b/Therm.c @@ -0,0 +1,410 @@ +#include "therm.h" +#include +#include "Detail.h" + + +void Therm_Init(Therm *therm) { + therm->CAL_TH = 4.1840; + therm->coefA = 0; + therm->coefB = 1; + therm->coefC = 2; + therm->coefD = 3; + therm->coefE = 4; + therm->coefF = 5; + therm->coefG = 6; + therm->coefH = 7; + therm->coefI = 8; + therm->coefJ = 9; + therm->coefK = 10; + + therm->dPdD = 0.0; + therm->dPdT = 0.0; + therm->dSi = 0.0; + therm->dTold = 0.0; + therm->dMrxold = 0.0; + + therm->GK_points = 5; + + + therm->GK_root[0] = 0.14887433898163121088; + therm->GK_root[1] = 0.43339539412924719080; + therm->GK_root[2] = 0.67940956829902440263; + therm->GK_root[3] = 0.86506336668898451073; + therm->GK_root[4] = 0.97390652851717172008; + + + therm->GK_weight[0] = 0.29552422471475286217; + therm->GK_weight[1] = 0.26926671930999634918; + therm->GK_weight[2] = 0.21908636251598204295; + therm->GK_weight[3] = 0.14945134915058059038; + therm->GK_weight[4] = 0.066671344308688137179; + + + double thermConstants[21][11] = { + {-29776.4, 7.95454, 43.9417, 1037.09, 1.56373, 813.205, -24.9027, 1019.98, -10.1601, 1070.14, -20.0615}, + {-3495.34, 6.95587, 0.272892, 662.738, -0.291318, -680.562, 1.78980, 1740.06, 0.0, 100.0, 4.49823}, + {20.7307, 6.96237, 2.68645, 500.371, -2.56429, -530.443, 3.91921, 500.198, 2.13290, 2197.22, 5.81381}, + {-37524.4, 7.98139, 24.3668, 752.320, 3.53990, 272.846, 8.44724, 1020.13, -13.2732, 869.510, -22.4010}, + {-56072.1, 8.14319, 37.0629, 735.402, 9.38159, 247.190, 13.4556, 1454.78, -11.7342, 984.518, -24.0426}, + {-13773.1, 7.97183, 6.27078, 2572.63, 2.05010, 1156.72, 0.0, 100.0, 0.0, 100.0, -3.24989}, + {-10085.4, 7.94680, -0.08380, 433.801, 2.85539, 843.792, 6.31595, 1481.43, -2.88457, 1102.23, -0.51551}, + {-5565.60, 6.66789, 2.33458, 2584.98, 0.749019, 559.656, 0.0, 100.0, 0.0, 100.0, -7.94821}, + {-2753.49, 6.95854, 2.02441, 1541.22, 0.096774, 3674.81, 0.0, 100.0, 0.0, 100.0, 6.23387}, + {-3497.45, 6.96302, 2.40013, 2522.05, 2.21752, 1154.15, 0.0, 100.0, 0.0, 100.0, 9.19749}, + {-72387.0, 17.8143, 58.2062, 1787.39, 40.7621, 808.645, 0.0, 100.0, 0.0, 100.0, -44.1341}, + {-72674.8, 18.6383, 57.4178, 1792.73, 38.6599, 814.151, 0.0, 100.0, 0.0, 100.0, -46.1938}, + {-91505.5, 21.3861, 74.3410, 1701.58, 47.0587, 775.899, 0.0, 100.0, 0.0, 100.0, -60.2474}, + {-83845.2, 22.5012, 69.5789, 1719.58, 46.2164, 802.174, 0.0, 100.0, 0.0, 100.0, -62.2197}, + {-94982.5, 26.6225, 80.3819, 1718.49, 55.6598, 802.069, 0.0, 100.0, 0.0, 100.0, -77.5366}, + {-103353.0, 30.4029, 90.6941, 1669.32, 63.2028, 786.001, 0.0, 100.0, 0.0, 100.0, -92.0164}, + {-109674.0, 34.0847, 100.253, 1611.55, 69.7675, 768.847, 0.0, 100.0, 0.0, 100.0, -106.149}, + {-122599.0, 38.5014, 111.446, 1646.48, 80.5015, 781.588, 0.0, 100.0, 0.0, 100.0, -122.444}, + {-133564.0, 42.7143, 122.173, 1654.85, 90.2255, 785.564, 0.0, 100.0, 0.0, 100.0, -138.006}, + {0.0, 4.9680, 0.0, 100.0, 0.0, 100.0, 0.0, 100.0, 0.0, 100.0, 0.0}, + {0.0, 4.9680, 0.0, 100.0, 0.0, 100.0, 0.0, 100.0, 0.0, 100.0, 0.0} + }; + + + for (int i = 0; i < 21; i++) { + for (int j = 0; j < 11; j++) { + therm->ThermConstants[i][j] = thermConstants[i][j]; + } + } +} + + +void Therm_Run(Therm *therm, NGParSTRUCT *ptNGPar, Detail *detail) { + + double c, x, y, z; + Detail_Run(detail, ptNGPar); + Detail_dZdD(detail, ptNGPar->dDf); + Therm_CprCvrHS(therm, ptNGPar, detail); + ptNGPar->dk = ptNGPar->dCp / ptNGPar->dCv; + x = ptNGPar->dk * RGAS * 1000.0 * ptNGPar->dTf; + y = ptNGPar->dMrx; + z = ptNGPar->dZf + ptNGPar->dDf * detail->ddZdD; + c = (x / y) * z; + ptNGPar->dSOS = sqrt(c); + ptNGPar->dKappa = (c * ptNGPar->dRhof) / ptNGPar->dPf; +} + + +double Therm_CpiMolar(Therm *therm, NGParSTRUCT *ptNGPar) { + double cp = 0.0; + double Cpx; + double DT, FT, HT, JT; + double Dx, Fx, Hx, Jx; + double T; + + T = ptNGPar->dTf; + + for (int i = 0; i < NUMBEROFCOMPONENTS; i++) { + if (ptNGPar->adMixture[i] > 0) { + Cpx = 0.0; + DT = therm->ThermConstants[i][therm->coefD] / T; + FT = therm->ThermConstants[i][therm->coefF] / T; + HT = therm->ThermConstants[i][therm->coefH] / T; + JT = therm->ThermConstants[i][therm->coefJ] / T; + Dx = DT / sinh(DT); + Fx = FT / cosh(FT); + Hx = HT / sinh(HT); + Jx = JT / cosh(JT); + Cpx += therm->ThermConstants[i][therm->coefB]; + Cpx += therm->ThermConstants[i][therm->coefC] * Dx * Dx; + Cpx += therm->ThermConstants[i][therm->coefE] * Fx * Fx; + Cpx += therm->ThermConstants[i][therm->coefG] * Hx * Hx; + Cpx += therm->ThermConstants[i][therm->coefI] * Jx * Jx; + Cpx *= ptNGPar->adMixture[i]; + cp += Cpx; + } + } + cp *= therm->CAL_TH; + return cp; +} + + +double Therm_coth(double x) { + return 1.0 / tanh(x); +} + + +double Therm_Ho(Therm *therm, NGParSTRUCT *ptNGPar) { + double H = 0.0; + double Hx; + double DT, FT, HT, JT; + double cothDT, tanhFT, cothHT, tanhJT; + double T = ptNGPar->dTf; + + for (int i = 0; i < NUMBEROFCOMPONENTS; i++) { + if (ptNGPar->adMixture[i] <= 0.0) continue; + Hx = 0.0; + DT = therm->ThermConstants[i][therm->coefD] / T; + FT = therm->ThermConstants[i][therm->coefF] / T; + HT = therm->ThermConstants[i][therm->coefH] / T; + JT = therm->ThermConstants[i][therm->coefJ] / T; + cothDT = Therm_coth(DT); + tanhFT = tanh(FT); + cothHT = Therm_coth(HT); + tanhJT = tanh(JT); + Hx += therm->ThermConstants[i][therm->coefA]; + Hx += therm->ThermConstants[i][therm->coefB] * T; + Hx += therm->ThermConstants[i][therm->coefC] * therm->ThermConstants[i][therm->coefD] * cothDT; + Hx -= therm->ThermConstants[i][therm->coefE] * therm->ThermConstants[i][therm->coefF] * tanhFT; + Hx += therm->ThermConstants[i][therm->coefG] * therm->ThermConstants[i][therm->coefH] * cothHT; + Hx -= therm->ThermConstants[i][therm->coefI] * therm->ThermConstants[i][therm->coefJ] * tanhJT; + Hx *= ptNGPar->adMixture[i]; + H += Hx; + } + + H *= therm->CAL_TH; + H /= ptNGPar->dMrx; + return H * 1000.0; +} + + +double Therm_So(Therm *therm, NGParSTRUCT *ptNGPar) { + double S = 0.0; + double Sx; + double DT, FT, HT, JT; + double cothDT, tanhFT, cothHT, tanhJT; + double sinhDT, coshFT, sinhHT, coshJT; + double T = ptNGPar->dTf; + for (int i = 0; i < NUMBEROFCOMPONENTS; i++) { + if (ptNGPar->adMixture[i] <= 0.0) continue; + Sx = 0.0; + DT = therm->ThermConstants[i][therm->coefD] / T; + FT = therm->ThermConstants[i][therm->coefF] / T; + HT = therm->ThermConstants[i][therm->coefH] / T; + JT = therm->ThermConstants[i][therm->coefJ] / T; + cothDT = Therm_coth(DT); + tanhFT = tanh(FT); + cothHT = Therm_coth(HT); + tanhJT = tanh(JT); + sinhDT = sinh(DT); + coshFT = cosh(FT); + sinhHT = sinh(HT); + coshJT = cosh(JT); + Sx += therm->ThermConstants[i][therm->coefK]; + Sx += therm->ThermConstants[i][therm->coefB] * log(T); + Sx += therm->ThermConstants[i][therm->coefC] * (DT * cothDT - log(sinhDT)); + Sx -= therm->ThermConstants[i][therm->coefE] * (FT * tanhFT - log(coshFT)); + Sx += therm->ThermConstants[i][therm->coefG] * (HT * cothHT - log(sinhHT)); + Sx -= therm->ThermConstants[i][therm->coefI] * (JT * tanhJT - log(coshJT)); + + Sx *= ptNGPar->adMixture[i]; + S += Sx; + } + + S *= therm->CAL_TH; + S /= ptNGPar->dMrx; + return S * 1000.0; +} + + +void Therm_CprCvrHS(Therm *therm, NGParSTRUCT *ptNGPar, Detail *detail) { + double Cvinc = 0.0; + double Cvr, Cpr; + double Hinc = 0.0; + double Sinc = 0.0; + double Smixing = 0.0; + double Si; + + double Cp = Therm_CpiMolar(therm, ptNGPar); + ptNGPar->dHo = Therm_Ho(therm, ptNGPar); + Si = Therm_So(therm, ptNGPar); + ptNGPar->dCpi = (Cp * 1000.0) / ptNGPar->dMrx; + for (int i = 0; i < therm->GK_points; i++) { + double x = ptNGPar->dDf * (1.0 + therm->GK_root[i]) / 2.0; + Detail_zdetail(detail, x); + Detail_dZdT(detail, x); + Detail_d2ZdT2(detail, x); + Hinc += therm->GK_weight[i] * detail->ddZdT / x; + Cvinc += therm->GK_weight[i] * (2.0 * detail->ddZdT + ptNGPar->dTf * detail->dd2ZdT2) / x; + Sinc += therm->GK_weight[i] * (detail->dZ + ptNGPar->dTf * detail->ddZdT - 1.0) / x; + x = ptNGPar->dDf * (1.0 - therm->GK_root[i]) / 2.0; + Detail_zdetail(detail, x); + Detail_dZdT(detail, x); + Detail_d2ZdT2(detail, x); + Hinc += therm->GK_weight[i] * detail->ddZdT / x; + Cvinc += therm->GK_weight[i] * (2.0 * detail->ddZdT + ptNGPar->dTf * detail->dd2ZdT2) / x; + Sinc += therm->GK_weight[i] * (detail->dZ + ptNGPar->dTf * detail->ddZdT - 1.0) / x; + } + + Detail_zdetail(detail, ptNGPar->dDf); + Detail_dZdT(detail, ptNGPar->dDf); + Detail_d2ZdT2(detail, ptNGPar->dDf); + Cvr = Cp - RGAS * (1.0 + ptNGPar->dTf * Cvinc * 0.5 * ptNGPar->dDf); + double a = (ptNGPar->dZf + ptNGPar->dTf * detail->ddZdT); + double b = (ptNGPar->dZf + ptNGPar->dDf * detail->ddZdD); + double dPdT = RGAS * ptNGPar->dDf * a; + double dPdD = RGAS * ptNGPar->dTf * b; + Cpr = Cvr + RGAS * ((a * a) / b); + + Cpr /= ptNGPar->dMrx; + Cvr /= ptNGPar->dMrx; + + ptNGPar->dCv = Cvr * 1000.0; + ptNGPar->dCp = Cpr * 1000.0; + + ptNGPar->dH = ptNGPar->dHo + 1000.0 * RGAS * ptNGPar->dTf * ( + ptNGPar->dZf - 1.0 - ptNGPar->dTf * Hinc * 0.5 * ptNGPar->dDf) / ptNGPar->dMrx; + + for (int i = 0; i < NUMBEROFCOMPONENTS; i++) { + if (ptNGPar->adMixture[i] != 0) Smixing += ptNGPar->adMixture[i] * log(ptNGPar->adMixture[i]); + } + Smixing *= RGAS; + + ptNGPar->dS = Si - Smixing - 1000.0 * RGAS * ( + log(ptNGPar->dPf / 101325.0) - log(ptNGPar->dZf) + Sinc * 0.5 * ptNGPar->dDf) / ptNGPar->dMrx; +} + + +void Therm_HS_Mode(Therm *therm, NGParSTRUCT *ptNGPar, Detail *detail, double H, double S, int bGuess) { + double s0 = S; + double h0 = H; + double t1, t2, tmin, tmax; + double p1, p2, px, pmin, pmax; + double delta1, delta2; + double tolerance = 0.001; + + if (bGuess) { + t1 = ptNGPar->dTf; + px = ptNGPar->dPf; + pmax = px * 2.0; + pmin = px * 0.1; + tmax = t1 * 1.5; + tmin = t1 * 0.67; + } else { + t1 = 273.15; + px = 1013250.0; + pmax = P_MAX; + pmin = 10000.0; + tmax = T_MAX; + tmin = T_MIN; + } + + t2 = t1 + 10.0; + + Detail_Run(detail, ptNGPar); + double h1 = Therm_H(therm, ptNGPar, detail) - h0; + + for (int i = 0; i < MAX_NUM_OF_ITERATIONS; i++) { + ptNGPar->dTf = t2; + p1 = px; + p2 = px * 0.1; + ptNGPar->dPf = p1; + Detail_Run(detail, ptNGPar); + double s1 = Therm_S(therm, ptNGPar, detail) - s0; + + for (int j = 0; j < MAX_NUM_OF_ITERATIONS; j++) { + ptNGPar->dPf = p2; + Detail_Run(detail, ptNGPar); + double s2 = Therm_S(therm, ptNGPar, detail) - s0; + delta2 = fabs(s1 - s2) / s0; + if (delta2 < tolerance) break; + double p0 = p2; + p2 = (p1 * s2 - p2 * s1) / (s2 - s1); + if (p2 <= pmin) p2 = pmin; + if (p2 >= pmax) p2 = pmax; + p1 = p0; + s1 = s2; + } + if (ptNGPar->lStatus == MAX_NUM_OF_ITERATIONS_EXCEEDED) break; + double h2 = Therm_H(therm, ptNGPar, detail) - h0; + delta1 = fabs(h1 - h2) / h0; + if (delta1 < tolerance && i > 0) break; + double t0 = t2; + t2 = (t1 * h2 - t2 * h1) / (h2 - h1); + if (t2 >= tmax) t2 = tmax; + if (t2 <= tmin) { + t2 = t0 + 10.0; + ptNGPar->dTf = t2; + Detail_Run(detail, ptNGPar); + h2 = Therm_H(therm, ptNGPar, detail) - h0; + } + + t1 = t0; + h1 = h2; + } + + if (ptNGPar->lStatus == MAX_NUM_OF_ITERATIONS_EXCEEDED) { + ptNGPar->lStatus = MAX_NUM_OF_ITERATIONS_EXCEEDED; + } +} + + +double Therm_H(Therm *therm, NGParSTRUCT *ptNGPar, Detail *detail) { + double Hinc = 0.0; + ptNGPar->dHo = Therm_Ho(therm, ptNGPar); + + for (int i = 0; i < therm->GK_points; i++) { + double x = ptNGPar->dDf * (1.0 + therm->GK_root[i]) / 2.0; + Detail_zdetail(detail, x); + Detail_dZdT(detail, x); + Detail_d2ZdT2(detail, x); + Hinc += therm->GK_weight[i] * detail->ddZdT / x; + if (i == 10) break; + x = ptNGPar->dDf * (1.0 - therm->GK_root[i]) / 2.0; + Detail_zdetail(detail, x); + Detail_dZdT(detail, x); + Detail_d2ZdT2(detail, x); + + Hinc += therm->GK_weight[i] * detail->ddZdT / x; + } + + Detail_zdetail(detail, ptNGPar->dDf); + Detail_dZdT(detail, ptNGPar->dDf); + Detail_d2ZdT2(detail, ptNGPar->dDf); + + ptNGPar->dH = ptNGPar->dHo + 1000.0 * RGAS * ptNGPar->dTf * ( + ptNGPar->dZf - 1.0 - ptNGPar->dTf * Hinc * 0.5 * ptNGPar->dDf) / ptNGPar->dMrx; + return ptNGPar->dH; +} + + +double Therm_S(Therm *therm, NGParSTRUCT *ptNGPar, Detail *detail) { + double Sinc; + double Smixing; + double x; + int i; + Sinc = 0.0; + Smixing = 0.0; + for (i = 0; i < therm->GK_points; i++) { + + x = ptNGPar->dDf * (1.0 + therm->GK_root[i]) / 2.0; + Detail_zdetail(detail, x); + Detail_dZdT(detail, x); + Detail_d2ZdT2(detail, x); + Sinc += therm->GK_weight[i] * (detail->dZ + ptNGPar->dTf * detail->ddZdT - 1.0) / x; + if (i == 10) break; + + x = ptNGPar->dDf * (1.0 - therm->GK_root[i]) / 2.0; + Detail_zdetail(detail, x); + Detail_dZdT(detail, x); + Detail_d2ZdT2(detail, x); + Sinc += therm->GK_weight[i] * (detail->dZ + ptNGPar->dTf * detail->ddZdT - 1.0) / x; + } + + + Detail_zdetail(detail, ptNGPar->dDf); + Detail_dZdT(detail, ptNGPar->dDf); + Detail_d2ZdT2(detail, ptNGPar->dDf); + + + if (ptNGPar->dTf != therm->dTold || ptNGPar->dMrx != therm->dMrxold) { + therm->dSi = Therm_So(therm, ptNGPar); + therm->dTold = ptNGPar->dTf; + therm->dMrxold = ptNGPar->dMrx; + } + + + for (i = 0; i < NUMBEROFCOMPONENTS; i++) { + if (ptNGPar->adMixture[i] != 0) Smixing += ptNGPar->adMixture[i] * log(ptNGPar->adMixture[i]); + } + Smixing *= RGAS; + + + ptNGPar->dS = therm->dSi - Smixing - 1000.0 * RGAS * ( + log(ptNGPar->dPf / 101325.0) - log(ptNGPar->dZf) + Sinc * 0.5 * ptNGPar->dDf) / ptNGPar->dMrx; + return (ptNGPar->dS); +} diff --git a/Therm.h b/Therm.h new file mode 100644 index 0000000..47962a4 --- /dev/null +++ b/Therm.h @@ -0,0 +1,60 @@ +/************************************************************************* +* �ļ�: therm.h +* ����: Therm���ͷ�ļ� +**************************************************************************/ + +#ifndef _THERM_H +#define _THERM_H + +#include "NGCal.h" +#include "Detail.h" + +typedef struct Therm { + double CAL_TH; + int coefA; + int coefB; + int coefC; + int coefD; + int coefE; + int coefF; + int coefG; + int coefH; + int coefI; + int coefJ; + int coefK; + + double dPdD; + double dPdT; + double dSi; + double dTold; + double dMrxold; + + int GK_points; + double GK_root[5]; + double GK_weight[5]; + + double ThermConstants[21][11]; +} Therm; + + +void Therm_Init(Therm *therm); + +void Therm_Run(Therm *therm, NGParSTRUCT *ptNGPar, Detail *detail); + +double Therm_CpiMolar(Therm *therm, NGParSTRUCT *ptNGPar); + +double Therm_coth(double x); + +double Therm_Ho(Therm *therm, NGParSTRUCT *ptNGPar); + +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); + +double Therm_H(Therm *therm, NGParSTRUCT *ptNGPar, Detail *detail); + +double Therm_S(Therm *therm, NGParSTRUCT *ptNGPar, Detail *detail); + +#endif diff --git a/main.c b/main.c new file mode 100644 index 0000000..e40086a --- /dev/null +++ b/main.c @@ -0,0 +1,227 @@ +#include +#include "NGCal.h" +#include "FlowCal.h" + +int main() { + // 定义并初始化 FlowParSTRUCT 结构体变量 + 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; // 温度(°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 = 9; // 镍铬合金 + + // 设置天然气组分(示例: 95%甲烷，5%其他) + // 初始化天然气组分数组(GB/T 21446-2008 典型示例组成) + for (int i = 0; i < NUMBEROFCOMPONENTS; i++) { + 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) + // 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[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] = 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) + // 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); +} diff --git a/main.h b/main.h new file mode 100644 index 0000000..1476f32 --- /dev/null +++ b/main.h @@ -0,0 +1,9 @@ +#ifndef __MAIN_H +#define __MAIN_H + +#include "stm32h7xx_hal.h" + +static void SystemClock_Config(void); + +#endif /* __MAIN_H */ + diff --git a/main.obj b/main.obj new file mode 100644 index 0000000..afda349 Binary files /dev/null and b/main.obj differ