#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); }