diff --git a/NG/NGCal.c b/NG/NGCal.c new file mode 100644 index 0000000..a74767f --- /dev/null +++ b/NG/NGCal.c @@ -0,0 +1,89 @@ +#include "NGCal.h" +#include "Therm.h" +#include "Detail.h" +#include "String.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) { + double DH, DDH, S, H; + double tolerance = 1.0; + double R, P, T, Z; + int i; + if (NULL == ptDetail || NULL == ptTherm) { + NGCal_UnInit(); + if (NGCal_NGCal != NGCal_Init(ptNGPar)) { + ptNGPar->lStatus = MEMORY_ALLOCATION_ERROR; + return 0.0; + } + } + + + Therm_Run(ptTherm, ptNGPar, ptDetail); + + DH = (ptNGPar->dSOS * ptNGPar->dSOS - dPlenumVelocity * dPlenumVelocity) / 2.0; + S = ptNGPar->dS; + H = ptNGPar->dH; + R = ptNGPar->dRhof; + P = ptNGPar->dPf; + Z = ptNGPar->dZf; + T = ptNGPar->dTf; + DDH = 10.0; + for (i = 1; i < MAX_NUM_OF_ITERATIONS; i++) { + Therm_HS_Mode(ptTherm, ptNGPar, ptDetail, H - DH, S, true); + Therm_Run(ptTherm, ptNGPar, ptDetail); + 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); + return ptNGPar->dCstar; +} + +double Cperf(NGParSTRUCT *ptNGPar) { + double k, root, exponent; + k = ptNGPar->dKappa; + root = 2.0 / (k + 1.0); + exponent = (k + 1.0) / (k - 1.0); + return (sqrt(k * pow(root, exponent))); +} + +double CRi(NGParSTRUCT *ptNGPar) { + return (Cperf(ptNGPar) / sqrt(ptNGPar->dZf)); +} diff --git a/NG/NGCal.h b/NG/NGCal.h new file mode 100644 index 0000000..3f9099d --- /dev/null +++ b/NG/NGCal.h @@ -0,0 +1,112 @@ +/************************************************************************* +* �ļ�: NGCal.h + +**************************************************************************/ + +#ifndef _NGCal_H +#define _NGCal_H + + + + +#include +#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 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/NG/Therm.c b/NG/Therm.c new file mode 100644 index 0000000..99087c1 --- /dev/null +++ b/NG/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, bool 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/NG/Therm.h b/NG/Therm.h new file mode 100644 index 0000000..37824c0 --- /dev/null +++ b/NG/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, bool bGuess); + +double Therm_H(Therm *therm, NGParSTRUCT *ptNGPar, Detail *detail); + +double Therm_S(Therm *therm, NGParSTRUCT *ptNGPar, Detail *detail); + +#endif