From 2fe122d9cb4ef74dc4bd4221f6bffc59b879f52f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?=E5=BB=96=E5=BE=B7=E4=BA=91?= Date: Sat, 12 Jul 2025 17:09:07 +0800 Subject: [PATCH] =?UTF-8?q?=E5=A4=A9=E7=84=B6=E6=B0=94=E6=B5=81=E9=87=8F?= =?UTF-8?q?=E8=AE=A1=E7=AE=97=EF=BC=8C=E5=8E=8B=E7=BC=A9=E5=9B=A0=E5=AD=90?= =?UTF-8?q?=EF=BC=8C=E5=A3=B0=E9=80=9F=E8=AE=A1=E7=AE=97=E7=9A=84C?= =?UTF-8?q?=E8=AF=AD=E8=A8=80=E7=89=88=E6=9C=AC?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- .idea/.gitignore | 8 + .idea/NG.iml | 2 + .idea/editor.xml | 344 +++++++++++++ .idea/encodings.xml | 6 + .idea/misc.xml | 7 + .idea/modules.xml | 8 + .idea/vcs.xml | 6 + CMakeLists.txt | 15 + Detail.c | 1133 +++++++++++++++++++++++++++++++++++++++++++ Detail.h | 120 +++++ FlowCal.c | 501 +++++++++++++++++++ FlowCal.h | 87 ++++ NGCal.c | 86 ++++ NGCal.h | 80 +++ Therm.c | 410 ++++++++++++++++ Therm.h | 60 +++ main.c | 227 +++++++++ main.h | 9 + main.obj | Bin 0 -> 5149 bytes 19 files changed, 3109 insertions(+) create mode 100644 .idea/.gitignore create mode 100644 .idea/NG.iml create mode 100644 .idea/editor.xml create mode 100644 .idea/encodings.xml create mode 100644 .idea/misc.xml create mode 100644 .idea/modules.xml create mode 100644 .idea/vcs.xml create mode 100644 CMakeLists.txt create mode 100644 Detail.c create mode 100644 Detail.h create mode 100644 FlowCal.c create mode 100644 FlowCal.h create mode 100644 NGCal.c create mode 100644 NGCal.h create mode 100644 Therm.c create mode 100644 Therm.h create mode 100644 main.c create mode 100644 main.h create mode 100644 main.obj 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 0000000000000000000000000000000000000000..afda34947bb421933a131569c242919201d938c1 GIT binary patch literal 5149 zcma)A4QvzF9e)>x5GZTRN4t=&)iqpM=rm(+>cj~pmpBPW>yU^9K|4C*_$vn6=Q?(p zP8thJAf#@zP^ylp(wbIEso2`?qp*U;Y85tU$5>mosT)=4CxneLplC%DZ93D+T^EJkBG4Aiz(w z-l%f}e+4x4Gf+iH5&R12uz&5Kc`;wtAMmAo=wtz_JD_QFPPW=*Tf(e^AtP^s?zpDK zA)>kh<(->id6!W)lht~&^5#MIU7%U2>2O?RQn-D~K%?;Gp{`Gjrn;l75_nqS7t%C2 zvd9ka*EB3-D`fWN>0}~jZAr5b$M0Ic&MM1=be3CkHBunoM0d*+b#tevyJL#F1yj_0 zYl=F@6m|85bpD=Az0T^~fMo%>^K>|OL{VJa+|t^(#?$3l*;H2v4|iKvrI~E^G_PuG zYw4^Mam8MQ-=pwz5Mee%q@`|s+vCluKOGCkQ@!gAlt`)@VRKl&YQ_4PFA}$@-5W_A zSvH$=6Y>t3OYS7c(j*tRMf^#%R}H1as@03c_CzWYi+nE<4~uR!mWV`yNvum=t86(p z29sUMOvAL|IIN5eNdv$fX_iE!%GPEtmCbIJs&ru2L9GttOJ!DSr$WV)>M*6&nNk}} zsUA}*VM_ILDQI4 zVqDb_OqG@b2{LP;ymFQ&7Rsw)dG9rt@$C^fPuSbTmVdEF8GV{@ zqaI0|)ZWpP5;h^Elx==n8T}fPMh-T#a1J)La1J)Sa1M6kh8#y$$vjDyjN&*(YgV=# z#x>Fx&P zq};+BR}ZL$a&nGa4e0BXQ%335a0m!ur;KL=UY=*9St?gXx+JRt$x>zTTrC^Bj7M6B zl;{D};mI^hijwJ)WF=F|1cQfb56y>y%4lIYsJLww6IKqENXw?RV6n<576H{&?T5i& za|6q+YQZ;0OFl|V28pP^wpRrz0`-atR)r}%rHkp<&Y?&s7f5+GKnQLY(Pue4EubWX zp;AD1QjT}%-+GWiup`XD5IPo6ZvGt>GDAp`q68!%;GT^0>2~kw`Ogzv8v9G@NBd1m*kr}u6-JN!ej>hoVue70JTJb!k}PSg)H&{~kK4L3*B3;P?OO%OiP;tN_# zWJ(iBi>tMGC5y9bwe&;sg6W~#zFF}`#l9c?YFDyxRq*ot&%Arjsgom3Z!Dhtag*b{ z&wfAeS@}p}y7ij7siOX3;?;YHt;;i+eS^<#yhzR`hV}GMEnO35H-ED5;HMA%^XXQ{ z_C@67A8f1dELzxfb=D5IwSPE#!t1X*_sW^%!=3YfQhVNW^*c>(hjvy)m$~A9-Fqn3 z@aW9!_LB~$^oz^?4(;o>KJB~T|I2UEJKp;5nn#rA(VC~m_dIjgzEc|}NcsWSmmk`3 zZPFn7;V`yE0>{$D!u9L;3guAHm$KRIV}dZCG4RP~iv@5iB#hxN2{@%K0$vl$fxjT2 z(;^I*1(4!D`f)B<1e{Vhfj)kn;Qj-F*Wkpq2spK{))O7`8`%dxuffG`5!~>^S_iGl zw{O&xPz^2v)Exi#fL0D9=;C3qNOK9f8}`D2N<{%+20!}Sy}7^ucn*9>6Ib(Q4mnxGy+sK!(m z-X4s)PF?-fYq(z1TvTJK3%`Xe0v<2|{8}D4{x?ku(VH5i8dF`gN4b7|dF}i`!}TW( zvR?BSCp_aV0&dC}oX=l+iy<$p6j&%UPz^;V;fZ5|p~!Txx)m6Uu$AV(mYv+XoI$B` z#gZj8PCVxb)ew!T5%ZZ^IiGRicP8(uvEvTT82%cD=NI>O{;;$S+sDEj97G7Hfi`oO zJnXgMdiTTfNyBxo=0c6xg(sxp8aVd%7YtXE=0c6xg@>Kt8f*NZ&v4<#1_soaU3kzK zF6Vglu;EH-F4UM^ICBhF_?Q3rqv3i&bD_rU!h_9leXutSn~FY$FK8~*m|eGACM}t}P|}YDA&XK| z_zKGKdZVh}7xnh00uj~Ab~djHw>dT53znO2pi z&{J9_-EFljiW|ImBZWS6rqrY_9F!ZDG&KCbzRI#oc0pT{L0?pMgt}cJyEEi;RaZM} x9JK}?wDZB>Lp*KQ`Dznu=lG?jEJN?O3198vi!${0)%dv#j?dQ{?|_Z|{{S;cm^c6c literal 0 HcmV?d00001