@param T Temperature, in K @param p Pressure (not used) @param psi_w Water mole fraction (mol_w/mol_ha) @param v_bar Mixture molar volume in m^3/mol_ha @returns s_ha Mixture molar entropy on a humid air basis in J/mol_ha/K */
| 1185 | @returns s_ha Mixture molar entropy on a humid air basis in J/mol_ha/K |
| 1186 | */ |
| 1187 | double MolarEntropy(double T, double p, double psi_w, double v_bar) { |
| 1188 | |
| 1189 | // vbar (molar volume) in m^3/mol |
| 1190 | double x1 = 0, x2 = 0, x3 = 0, y1 = 0, y2 = 0, eps = 1e-8, f = 999, R_bar_Lem = 8.314510; |
| 1191 | int iter = 1; |
| 1192 | double sbar_0 = NAN, sbar_a = 0, sbar_w = 0, sbar = NAN, R_bar = 8.314472, vbar_a_guess = NAN, Baa = NAN, Caaa = NAN, vbar_a = 0; |
| 1193 | double B = NAN, dBdT = NAN, C = NAN, dCdT = NAN; |
| 1194 | // Constant for entropy |
| 1195 | sbar_0 = 0.02366427495; //[J/mol/K] |
| 1196 | |
| 1197 | // Calculate vbar_a, the molar volume of dry air |
| 1198 | // B_m, C_m, etc. functions take care of the units |
| 1199 | Baa = B_m(T, 0); |
| 1200 | B = B_m(T, psi_w); |
| 1201 | dBdT = dB_m_dT(T, psi_w); |
| 1202 | Caaa = C_m(T, 0); |
| 1203 | C = C_m(T, psi_w); |
| 1204 | dCdT = dC_m_dT(T, psi_w); |
| 1205 | |
| 1206 | vbar_a_guess = R_bar_Lem * T / p; //[m^3/mol] since p in [Pa] |
| 1207 | |
| 1208 | while ((iter <= 3 || std::abs(f) > eps) && iter < 100) { |
| 1209 | if (iter == 1) { |
| 1210 | x1 = vbar_a_guess; |
| 1211 | vbar_a = x1; |
| 1212 | } |
| 1213 | if (iter == 2) { |
| 1214 | x2 = vbar_a_guess + 0.001; |
| 1215 | vbar_a = x2; |
| 1216 | } |
| 1217 | if (iter > 2) { |
| 1218 | vbar_a = x2; |
| 1219 | } |
| 1220 | f = R_bar_Lem * T / vbar_a * (1 + Baa / vbar_a + Caaa / pow(vbar_a, 2)) - p; |
| 1221 | if (iter == 1) { |
| 1222 | y1 = f; |
| 1223 | } |
| 1224 | if (iter > 1) { |
| 1225 | y2 = f; |
| 1226 | x3 = x2 - y2 / (y2 - y1) * (x2 - x1); |
| 1227 | y1 = y2; |
| 1228 | x1 = x2; |
| 1229 | x2 = x3; |
| 1230 | } |
| 1231 | iter = iter + 1; |
| 1232 | } |
| 1233 | |
| 1234 | if (FlagUseIdealGasEnthalpyCorrelations) { |
| 1235 | std::cout << "Not implemented" << '\n'; |
| 1236 | } else { |
| 1237 | sbar_w = IdealGasMolarEntropy_Water(T, p); |
| 1238 | sbar_a = IdealGasMolarEntropy_Air(T, vbar_a); |
| 1239 | } |
| 1240 | if (psi_w != 0) { |
| 1241 | sbar = sbar_0 + (1 - psi_w) * sbar_a + psi_w * sbar_w |
| 1242 | - R_bar * ((B + T * dBdT) / v_bar + (C + T * dCdT) / (2 * pow(v_bar, 2)) + (1 - psi_w) * log(1 - psi_w) + psi_w * log(psi_w)); |
| 1243 | } else { |
| 1244 | sbar = sbar_0 + sbar_a - R_bar * ((B + T * dBdT) / v_bar + (C + T * dCdT) / (2 * pow(v_bar, 2))); |
no test coverage detected