1 //#include "ThermoProperties.h"
2 #include "WaterHGKreaktoro.h"
3 #include "Reaktoro/WaterConstants.hpp"
4 #include "Common/OutputWaterSteamConventionProp.h"
5 #include "ThermoProperties.h"
6
7 namespace ThermoFun {
8
thermoPropertiesWaterHGKreaktoro(Reaktoro_::Temperature T,const WaterThermoState & wt)9 auto thermoPropertiesWaterHGKreaktoro(Reaktoro_::Temperature T, const WaterThermoState& wt) -> ThermoPropertiesSubstance
10 {
11 // Auxiliary data from Helgeson and Kirkham (1974), on page 1098
12 const auto Ttr = 273.15; // unit: K
13 const auto Str = 15.1320 * cal_to_J; // unit: J/(mol*K)
14 const auto Gtr = -56290.0 * cal_to_J; // unit: J/mol
15 const auto Htr = -68767.0 * cal_to_J; // unit: J/mol
16 const auto Utr = -67887.0 * cal_to_J; // unit: J/mol
17 const auto Atr = -55415.0 * cal_to_J; // unit: J/mol
18
19 const auto Sw = waterMolarMass * wt.entropy; // unit: J/(mol*K)
20 const auto Hw = waterMolarMass * wt.enthalpy; // unit: J/mol
21 const auto Uw = waterMolarMass * wt.internal_energy; // unit: J/mol
22
23 // Calculate the standard molal thermodynamic properties of the aqueous species
24 const auto S = Sw + Str;
25 const auto H = Hw + Htr;
26 const auto U = Uw + Utr;
27 const auto G = Hw - T * (Sw + Str) + Ttr * Str + Gtr;
28 const auto A = Uw - T * (Sw + Str) + Ttr * Str + Atr;
29 const auto V = wt.volume * waterMolarMass * 100000; // unit J/bar
30 const auto Cp = wt.cp * waterMolarMass;
31 const auto Cv = wt.cv * waterMolarMass;
32
33 ThermoPropertiesSubstance state;
34 state.entropy = S;
35 state.enthalpy = H;
36 state.internal_energy = U;
37 state.gibbs_energy = G;
38 state.helmholtz_energy = A;
39 state.volume = V;
40 state.heat_capacity_cp = Cp;
41 state.heat_capacity_cv = Cv;
42
43 #ifdef OutputSTEAM_CONVENTION
44 OutputSteamConventionH2OProp("H2OHGKreaktoro.csv", wt);
45 #endif
46
47 return state;
48 }
49
propertiesWaterHGKreaktoro(const WaterThermoState & wt)50 auto propertiesWaterHGKreaktoro(const WaterThermoState& wt) -> PropertiesSolvent
51 {
52 PropertiesSolvent state;
53
54 state.density = wt.density;
55 state.densityP = wt.densityP;
56 state.densityPP = wt.densityPP;
57 state.densityT = wt.densityT;
58 state.densityTP = wt.densityTP;
59 state.densityTT = wt.densityTT;
60 state.pressure = wt.pressure;
61 state.pressureD = wt.pressureD;
62 state.pressureDD = wt.pressureDD;
63 state.pressureT = wt.pressureT;
64 state.pressureTD = wt.pressureTD;
65 state.pressureTT = wt.pressureTT;
66
67 state.Alpha = -wt.densityT/wt.density;
68 state.Beta = wt.densityP/wt.density;
69 state.dAldT = -wt.densityTT/wt.density + state.Alpha*state.Alpha;
70 // const auto alphaP = -wt.densityTP/wt.density - alpha*beta;
71 // const auto betaP = wt.densityPP/wt.density - beta*beta;
72
73 // auto t = Reaktoro::Temperature( wt.temperature.val );
74 // waterIdealGas(t, state);
75
76 return state;
77 }
78
saturatedWaterVaporPressureHGK(Reaktoro_::Temperature TK)79 auto saturatedWaterVaporPressureHGK(Reaktoro_::Temperature TK) -> Reaktoro_::ThermoScalar
80 {
81 Reaktoro_::ThermoScalar pl, psHGK, v, w, b, q, z;
82 int i=-1;
83 double a[8] ={ -.78889166e1, .25514255e1, -.6716169e1, .33239495e2,
84 -.10538479e3, .17435319e3, -.14839348e3, .48631602e2};
85 if (TK.val <= 314.0e0)
86 {
87 pl = 6.3573118e0 - 8858.843e0 / TK + 607.56335e0 * pow(TK,-0.6e0);
88 psHGK = 0.1e0 * exp(pl);
89 }
90 else
91 {
92 v = TK / 647.25e0;
93 w = 1.0e0 - v;
94 if (w.val<0)
95 w=w*-1;
96
97 b = 0.0e0;
98 while (++i <= 7)
99 {
100 z = i + 1;
101 b += a[i] * pow(w,((z + 1.0e0) / 2.0e0));
102 }
103 q = b / v;
104 psHGK = 22.093e0 * exp(q);
105 }
106 return(psHGK*1e1);
107 }
108
109 } // end namespace ThermoFun
110