|
| 1 | + |
| 2 | +/** The phase of water. */ |
| 3 | +enum class WaterPhase |
| 4 | +{ |
| 5 | + solid, |
| 6 | + liquid, |
| 7 | + gas, |
| 8 | + supercritical |
| 9 | +}; |
| 10 | + |
| 11 | +/** @returns the String representation of the given water's phase. */ |
| 12 | +inline String toString (WaterPhase mp) noexcept |
| 13 | +{ |
| 14 | + switch (mp) |
| 15 | + { |
| 16 | + case WaterPhase::solid: return TRANS ("Solid"); |
| 17 | + case WaterPhase::liquid: return TRANS ("Liquid"); |
| 18 | + case WaterPhase::gas: return TRANS ("Gas"); |
| 19 | + case WaterPhase::supercritical: return TRANS ("Supercritical Fluid"); |
| 20 | + default: return TRANS ("Unknown"); |
| 21 | + } |
| 22 | +} |
| 23 | + |
| 24 | +class WaterPhaseCalculator final |
| 25 | +{ |
| 26 | +public: |
| 27 | + /** Freezing point (°C) correlated with pressure (Pa). */ |
| 28 | + static double getFreezingPointCelsius (double pressurePa) |
| 29 | + { |
| 30 | + constexpr auto slope = -7.4e-8; // K/Pa near 1 atm |
| 31 | + constexpr auto Pref = kPaPerAtm; |
| 32 | + constexpr auto TrefC = 0.0; |
| 33 | + return TrefC + slope * (pressurePa - Pref); |
| 34 | + } |
| 35 | + |
| 36 | + /** Boiling point (°C) vs pressure (Pa) using Wagner correlation. */ |
| 37 | + static std::optional<double> getBoilingPointCelsius (double pressurePa, |
| 38 | + double tolK = 1e-5, |
| 39 | + int maxIters = 100) |
| 40 | + { |
| 41 | + if (pressurePa < Ptriple || pressurePa >= Pc) |
| 42 | + return std::nullopt; |
| 43 | + |
| 44 | + auto lo = Ttriple; |
| 45 | + auto hi = Tc - 1e-6; |
| 46 | + auto PloOpt = saturationPressurePa_Wagner (lo); |
| 47 | + auto PhiOpt = saturationPressurePa_Wagner (hi); |
| 48 | + if (!PloOpt || !PhiOpt) |
| 49 | + return std::nullopt; |
| 50 | + |
| 51 | + auto Plo = *PloOpt; |
| 52 | + auto Phi = *PhiOpt; |
| 53 | + |
| 54 | + if (!(Plo <= pressurePa && pressurePa <= Phi)) |
| 55 | + return std::nullopt; |
| 56 | + |
| 57 | + for (int i = 0; i < maxIters && (hi - lo) > tolK; ++i) |
| 58 | + { |
| 59 | + auto mid = 0.5 * (lo + hi); |
| 60 | + auto PmidOpt = saturationPressurePa_Wagner(mid); |
| 61 | + |
| 62 | + if (!PmidOpt) |
| 63 | + return std::nullopt; |
| 64 | + |
| 65 | + auto Pmid = *PmidOpt; |
| 66 | + if (Pmid < pressurePa) lo = mid; |
| 67 | + else hi = mid; |
| 68 | + } |
| 69 | + |
| 70 | + const auto T_K = 0.5 * (lo + hi); |
| 71 | + return T_K - 273.15; |
| 72 | + } |
| 73 | + |
| 74 | + /** Phase determination. */ |
| 75 | + static WaterPhase getWaterPhase (double tempC, double pressurePa) |
| 76 | + { |
| 77 | + if (pressurePa < Ptriple) |
| 78 | + return (tempC < (Ttriple - 273.15)) |
| 79 | + ? WaterPhase::solid |
| 80 | + : WaterPhase::gas; |
| 81 | + |
| 82 | + if (pressurePa >= Pc && (tempC + 273.15) >= Tc) |
| 83 | + return WaterPhase::supercritical; |
| 84 | + |
| 85 | + const auto TfC = getFreezingPointCelsius (pressurePa); |
| 86 | + if (tempC < TfC) |
| 87 | + return WaterPhase::solid; |
| 88 | + |
| 89 | + const auto TbOpt = getBoilingPointCelsius (pressurePa); |
| 90 | + if (!TbOpt) |
| 91 | + return (tempC >= (Tc - 273.15)) |
| 92 | + ? WaterPhase::supercritical |
| 93 | + : WaterPhase::liquid; |
| 94 | + |
| 95 | + if (tempC < *TbOpt) |
| 96 | + return WaterPhase::liquid; |
| 97 | + |
| 98 | + return WaterPhase::gas; |
| 99 | + } |
| 100 | + |
| 101 | +private: |
| 102 | + static inline constexpr double kPaPerAtm = 101325.0; |
| 103 | + static inline constexpr double Tc = 647.096; // K |
| 104 | + static inline constexpr double Pc = 22.064e6; // Pa |
| 105 | + static inline constexpr double Ttriple = 273.16; // K |
| 106 | + static inline constexpr double Ptriple = 611.657; // Pa |
| 107 | + |
| 108 | + /** Wagner saturation pressure. */ |
| 109 | + static std::optional<double> saturationPressurePa_Wagner (double T_K) |
| 110 | + { |
| 111 | + if (T_K < Ttriple || T_K > Tc) |
| 112 | + return std::nullopt; |
| 113 | + |
| 114 | + constexpr auto a = -7.85951783; |
| 115 | + constexpr auto b = 1.84408259; |
| 116 | + constexpr auto c = -11.7866497; |
| 117 | + constexpr auto d = 22.6807411; |
| 118 | + constexpr auto e = -15.9618719; |
| 119 | + constexpr auto f = 1.80122502; |
| 120 | + |
| 121 | + auto tau = 1.0 - T_K / Tc; |
| 122 | + auto poly = a * tau |
| 123 | + + b * std::pow (tau, 1.5) |
| 124 | + + c * std::pow (tau, 3.0) |
| 125 | + + d * std::pow (tau, 3.5) |
| 126 | + + e * std::pow (tau, 4.0) |
| 127 | + + f * std::pow (tau, 7.5); |
| 128 | + |
| 129 | + const auto lnPratio = (Tc / T_K) * poly; |
| 130 | + return Pc * std::exp (lnPratio); |
| 131 | + } |
| 132 | +}; |
0 commit comments