From 6bf362f2122a5cf8b9e192136c1d6c93fa623fc5 Mon Sep 17 00:00:00 2001 From: gdmiron Date: Tue, 23 Jul 2024 11:42:42 +0200 Subject: [PATCH] explicitly set the state of the solvent (0 liquid, 1 vapor), default -1: from database record --- .gitignore | 2 +- Resources/databases/misc/water-thermofun.json | 447 ++++++++++++++++++ ThermoFun/ElectroModelsSolvent.cpp | 16 +- ThermoFun/ElectroModelsSolvent.h | 12 +- .../WaterElectroStateJohnsonNorton.cpp | 2 +- .../WaterElectroStateJohnsonNorton.hpp | 2 +- .../Solvent/WaterElectroFernandez1997.cpp | 12 +- .../Solvent/WaterElectroFernandez1997.h | 2 +- .../Solvent/WaterElectroSverjensky2014.cpp | 12 +- .../Solvent/WaterElectroSverjensky2014.h | 2 +- .../Substances/Solvent/WaterHGK-JNgems.cpp | 9 +- ThermoFun/ThermoEngine.cpp | 71 +-- ThermoFun/ThermoEngine.h | 12 +- ThermoFun/ThermoModelsSolvent.h | 4 +- pytests/Substances/Solvent/test_water_EoS.py | 60 +++ .../Substances/Solvent/water-thermofun.json | 447 ++++++++++++++++++ python/pyThermoFun/pyThermoEngine.cpp | 53 ++- tests/interfaceTest/src/main.cpp | 32 +- 18 files changed, 1117 insertions(+), 80 deletions(-) create mode 100644 Resources/databases/misc/water-thermofun.json create mode 100644 pytests/Substances/Solvent/test_water_EoS.py create mode 100644 pytests/Substances/Solvent/water-thermofun.json diff --git a/.gitignore b/.gitignore index 17e7204f..a45c5a2b 100644 --- a/.gitignore +++ b/.gitignore @@ -63,7 +63,7 @@ python/tests/results_dbc.csv python/pyThermoFun/CMakeLists_old.txt *.log - +tests/build* diff --git a/Resources/databases/misc/water-thermofun.json b/Resources/databases/misc/water-thermofun.json new file mode 100644 index 00000000..54cf8468 --- /dev/null +++ b/Resources/databases/misc/water-thermofun.json @@ -0,0 +1,447 @@ +{ + "thermodataset": "aq17-gem", + "datasources": [ + "db.thermohub.org", + "[2017MIR/WAG]" + ], + "date": "10.12.2019 12:58:35", + "substances": [ + { + "name": "Al(OH)2+", + "symbol": "my_Al(OH)2+", + "formula": "Al(OH)2+", + "formula_charge": 1, + "aggregate_state": { + "4": "AS_AQUEOUS" + }, + "class_": { + "2": "SC_AQSOLUTE" + }, + "Tst": 298.15, + "Pst": 100000, + "TPMethods": [ + { + "method": { + "3": "solute_hkf88_reaktoro" + }, + "eos_hkf_coeffs": { + "values": [ + 0.24940000474453, + -169.08999633789, + 6.4145998954773, + -27091, + 16.743900299072, + -10465, + 53240, + 0 + ], + "units": [ + "cal/(mol*bar)", + "cal/mol", + "(cal*K)/mol", + "cal/(mol*K)", + "(cal*K)/mol", + "cal/mol" + ], + "names": [ + "a1", + "a2", + "a3", + "a4", + "c1", + "c2", + "wref" + ] + } + } + ], + "sm_heat_capacity_p": { + "values": [ + 40.865230560303 + ], + "errors": [ + 0 + ], + "units": [ + "J/(mol*K)" + ] + }, + "sm_gibbs_energy": { + "values": [ + -898292 + ], + "units": [ + "J/mol" + ] + }, + "sm_enthalpy": { + "values": [ + -995581 + ], + "units": [ + "J/mol" + ] + }, + "sm_entropy_abs": { + "values": [ + -50 + ], + "units": [ + "J/(mol*K)" + ] + }, + "sm_volume": { + "values": [ + 0.38507527112961 + ], + "errors": [ + 0 + ], + "units": [ + "J/bar" + ] + }, + "datasources": [ + "[2001TAG/SCH]", + "[2017MIR/WAG]" + ] + }, + { + "name": "Water HGK", + "symbol": "H2O@", + "formula": "H2O@", + "formula_charge": 0, + "aggregate_state": { + "4": "AS_AQUEOUS" + }, + "class_": { + "3": "SC_AQSOLVENT" + }, + "Tst": 298.15, + "Pst": 100000, + "TPMethods": [ + { + "method": { + "29": "water_eos_hgk84_lvs83_gems" + } + }, + { + "method": { + "25": "water_diel_jnort91_reaktoro" + } + } + ], + "sm_heat_capacity_p": { + "values": [ + 75.360527038574 + ], + "errors": [ + 0 + ], + "units": [ + "J/(mol*K)" + ] + }, + "sm_gibbs_energy": { + "values": [ + -237183 + ], + "errors": [ + 100 + ], + "units": [ + "J/mol" + ] + }, + "sm_enthalpy": { + "values": [ + -285881 + ], + "errors": [ + 200 + ], + "units": [ + "J/mol" + ] + }, + "sm_entropy_abs": { + "values": [ + 69.922996520996 + ], + "errors": [ + 0.10000000149012 + ], + "units": [ + "J/(mol*K)" + ] + }, + "sm_volume": { + "values": [ + 1.8068397045136 + ], + "errors": [ + 0 + ], + "units": [ + "J/bar" + ] + }, + "datasources": [ + "[1992JOH/OEL]" + ] + }, + { + "name": "Water HGK", + "symbol": "H2Og", + "formula": "H2O@", + "formula_charge": 0, + "aggregate_state": { + "0": "AS_GAS" + }, + "class_": { + "3": "SC_AQSOLVENT" + }, + "Tst": 298.15, + "Pst": 100000, + "TPMethods": [ + { + "method": { + "29": "water_eos_hgk84_lvs83_gems" + } + }, + { + "method": { + "26": "water_diel_jnort91_gems" + } + } + ], + "sm_heat_capacity_p": { + "values": [ + 75.360527038574 + ], + "errors": [ + 0 + ], + "units": [ + "J/(mol*K)" + ] + }, + "sm_gibbs_energy": { + "values": [ + -237183 + ], + "errors": [ + 100 + ], + "units": [ + "J/mol" + ] + }, + "sm_enthalpy": { + "values": [ + -285881 + ], + "errors": [ + 200 + ], + "units": [ + "J/mol" + ] + }, + "sm_entropy_abs": { + "values": [ + 69.922996520996 + ], + "errors": [ + 0.10000000149012 + ], + "units": [ + "J/(mol*K)" + ] + }, + "sm_volume": { + "values": [ + 1.8068397045136 + ], + "errors": [ + 0 + ], + "units": [ + "J/bar" + ] + }, + "datasources": [ + "[1992JOH/OEL]" + ] + }, + { + "name": "Water HGK", + "symbol": "H2O@", + "formula": "H2O@", + "formula_charge": 0, + "aggregate_state": { + "4": "AS_AQUEOUS" + }, + "class_": { + "3": "SC_AQSOLVENT" + }, + "Tst": 298.15, + "Pst": 100000, + "TPMethods": [ + { + "method": { + "29": "water_eos_hgk84_lvs83_gems" + } + }, + { + "method": { + "25": "water_diel_jnort91_reaktoro" + } + } + ], + "sm_heat_capacity_p": { + "values": [ + 75.360527038574 + ], + "errors": [ + 0 + ], + "units": [ + "J/(mol*K)" + ] + }, + "sm_gibbs_energy": { + "values": [ + -237183 + ], + "errors": [ + 100 + ], + "units": [ + "J/mol" + ] + }, + "sm_enthalpy": { + "values": [ + -285881 + ], + "errors": [ + 200 + ], + "units": [ + "J/mol" + ] + }, + "sm_entropy_abs": { + "values": [ + 69.922996520996 + ], + "errors": [ + 0.10000000149012 + ], + "units": [ + "J/(mol*K)" + ] + }, + "sm_volume": { + "values": [ + 1.8068397045136 + ], + "errors": [ + 0 + ], + "units": [ + "J/bar" + ] + }, + "datasources": [ + "[1992JOH/OEL]" + ] + }, + { + "name": "Water HGK", + "symbol": "H2O@reak", + "formula": "H2O@", + "formula_charge": 0, + "aggregate_state": { + "4": "AS_AQUEOUS" + }, + "class_": { + "3": "SC_AQSOLVENT" + }, + "Tst": 298.15, + "Pst": 100000, + "TPMethods": [ + { + "method": { + "32": "water_eos_iapws95_reaktoro" + } + }, + { + "method": { + "25": "water_diel_jnort91_reaktoro" + } + } + ], + "sm_heat_capacity_p": { + "values": [ + 75.360527038574 + ], + "errors": [ + 0 + ], + "units": [ + "J/(mol*K)" + ] + }, + "sm_gibbs_energy": { + "values": [ + -237183 + ], + "errors": [ + 100 + ], + "units": [ + "J/mol" + ] + }, + "sm_enthalpy": { + "values": [ + -285881 + ], + "errors": [ + 200 + ], + "units": [ + "J/mol" + ] + }, + "sm_entropy_abs": { + "values": [ + 69.922996520996 + ], + "errors": [ + 0.10000000149012 + ], + "units": [ + "J/(mol*K)" + ] + }, + "sm_volume": { + "values": [ + 1.8068397045136 + ], + "errors": [ + 0 + ], + "units": [ + "J/bar" + ] + }, + "datasources": [ + "[1992JOH/OEL]" + ] + } + ], + "reactions": [] +} diff --git a/ThermoFun/ElectroModelsSolvent.cpp b/ThermoFun/ElectroModelsSolvent.cpp index 74a89a4f..56dcb8a6 100644 --- a/ThermoFun/ElectroModelsSolvent.cpp +++ b/ThermoFun/ElectroModelsSolvent.cpp @@ -39,7 +39,7 @@ WaterJNreaktoro::WaterJNreaktoro(const Substance &substance) {} // calculation -auto WaterJNreaktoro::electroPropertiesSolvent(double T, double P, PropertiesSolvent ps) -> ElectroPropertiesSolvent +auto WaterJNreaktoro::electroPropertiesSolvent(double T, double P, PropertiesSolvent ps, int state) -> ElectroPropertiesSolvent { // if (P==0) P = saturatedWaterVaporPressureHGK(T+C_to_K); @@ -57,7 +57,7 @@ auto WaterJNreaktoro::electroPropertiesSolvent(double T, double P, PropertiesSol wts.densityTP = ps.densityTP; wts.densityPP = ps.densityPP; - WaterElectroState wes = waterElectroStateJohnsonNorton(t, /*p,*/ wts); + WaterElectroState wes = waterElectroStateJohnsonNorton(t, /*p,*/ wts, state); return electroPropertiesWaterJNreaktoro(wes); } @@ -88,13 +88,13 @@ WaterJNgems::WaterJNgems(const Substance &substance) {} // calculation -auto WaterJNgems::electroPropertiesSolvent(double T, double P) -> ElectroPropertiesSolvent +auto WaterJNgems::electroPropertiesSolvent(double T, double P, int state) -> ElectroPropertiesSolvent { WaterHGKgems water_hgk; T -= C_to_K; P /= bar_to_Pa; water_hgk.calculateWaterHGKgems(T, P); - return water_hgk.electroPropertiesWaterJNgems(0); // state 0 = liquid + return water_hgk.electroPropertiesWaterJNgems(state); // state 0 = liquid } //======================================================================================================= @@ -126,14 +126,14 @@ WaterElectroSverjensky2014::WaterElectroSverjensky2014(const Substance &substanc {} // calculation -auto WaterElectroSverjensky2014::electroPropertiesSolvent(double T, double P/*, PropertiesSolvent ps*/) -> ElectroPropertiesSolvent +auto WaterElectroSverjensky2014::electroPropertiesSolvent(double T, double P/*, PropertiesSolvent ps*/, int state) -> ElectroPropertiesSolvent { // if (P==0) P = saturatedWaterVaporPressureHGK(T+C_to_K); auto t = Reaktoro_::Temperature(T); t -= C_to_K; auto p = Reaktoro_::Pressure(P); p /= bar_to_Pa; - return electroPropertiesWaterSverjensky2014(/*ps,*/ t, p, pimpl->substance); + return electroPropertiesWaterSverjensky2014(/*ps,*/ t, p, pimpl->substance, state); } //======================================================================================================= @@ -163,14 +163,14 @@ WaterElectroFernandez1997::WaterElectroFernandez1997(const Substance &substance) {} // calculation -auto WaterElectroFernandez1997::electroPropertiesSolvent(double T, double P/*, PropertiesSolvent ps*/) -> ElectroPropertiesSolvent +auto WaterElectroFernandez1997::electroPropertiesSolvent(double T, double P/*, PropertiesSolvent ps*/, int state) -> ElectroPropertiesSolvent { // if (P==0) P = saturatedWaterVaporPressureHGK(T+C_to_K); auto t = Reaktoro_::Temperature(T); t -= C_to_K; auto p = Reaktoro_::Pressure(P); p /= bar_to_Pa; - return electroPropertiesWaterFernandez1997(/*ps,*/ t, p, pimpl->substance); // t (celsius), p (bar) + return electroPropertiesWaterFernandez1997(/*ps,*/ t, p, pimpl->substance, state); // t (celsius), p (bar) } } // End namespace ThermoFun diff --git a/ThermoFun/ElectroModelsSolvent.h b/ThermoFun/ElectroModelsSolvent.h index 3ee14cfb..88c6523d 100644 --- a/ThermoFun/ElectroModelsSolvent.h +++ b/ThermoFun/ElectroModelsSolvent.h @@ -28,7 +28,8 @@ class WaterJNreaktoro /// @param T The temperature value (in units of C) /// @param P The pressure value (in units of bar) /// @param ps structure holding the solvent properties - auto electroPropertiesSolvent(double T, double P, PropertiesSolvent ps) -> ElectroPropertiesSolvent; + /// @param state The solvent state, 0 liquid, 1 gas/vapor, default 0 + auto electroPropertiesSolvent(double T, double P, PropertiesSolvent ps, int state=0) -> ElectroPropertiesSolvent; private: struct Impl; @@ -54,7 +55,8 @@ class WaterJNgems /// Return the electro-chemical properties of the solvent. /// @param T The temperature value (in units of C) /// @param P The pressure value (in units of bar) - auto electroPropertiesSolvent(double T, double P) -> ElectroPropertiesSolvent; + /// @param state The solvent state, 0 liquid, 1 gas/vapor, default 0 + auto electroPropertiesSolvent(double T, double P, int state=0) -> ElectroPropertiesSolvent; private: struct Impl; @@ -84,7 +86,8 @@ class WaterElectroSverjensky2014 /// @param T The temperature value (in units of C) /// @param P The pressure value (in units of bar) /// @param ps structure holding the solvent properties - auto electroPropertiesSolvent(double T, double P) -> ElectroPropertiesSolvent; + /// @param state The solvent state, 0 liquid, 1 gas/vapor, default -1 determined by the record key "aggregate_state" in the database "4": "AS_AQUEOUS" or "O": "AS_GAS" + auto electroPropertiesSolvent(double T, double P, int state = -1) -> ElectroPropertiesSolvent; private: struct Impl; @@ -112,7 +115,8 @@ class WaterElectroFernandez1997 /// @param T The temperature value (in units of C) /// @param P The pressure value (in units of bar) /// @param ps structure holding the solvent properties - auto electroPropertiesSolvent(double T, double P) -> ElectroPropertiesSolvent; + /// @param state The solvent state, 0 liquid, 1 gas/vapor, default -1 determined by the record key "aggregate_state" in the database "4": "AS_AQUEOUS" or "O": "AS_GAS" + auto electroPropertiesSolvent(double T, double P, int state = -1) -> ElectroPropertiesSolvent; private: struct Impl; diff --git a/ThermoFun/Substances/Solvent/Reaktoro/WaterElectroStateJohnsonNorton.cpp b/ThermoFun/Substances/Solvent/Reaktoro/WaterElectroStateJohnsonNorton.cpp index af5b9fa1..d8ce1368 100644 --- a/ThermoFun/Substances/Solvent/Reaktoro/WaterElectroStateJohnsonNorton.cpp +++ b/ThermoFun/Substances/Solvent/Reaktoro/WaterElectroStateJohnsonNorton.cpp @@ -76,7 +76,7 @@ Reaktoro_::ThermoScalar (*k_tt[5])(Reaktoro_::ThermoScalar) = {k0_tt, k1_tt, k2_ } // namespace -auto waterElectroStateJohnsonNorton(Reaktoro_::Temperature T, /*Pressure P,*/ const WaterThermoState& wt) -> WaterElectroState +auto waterElectroStateJohnsonNorton(Reaktoro_::Temperature T, /*Pressure P,*/ const WaterThermoState& wt, int state) -> WaterElectroState { WaterElectroState we; diff --git a/ThermoFun/Substances/Solvent/Reaktoro/WaterElectroStateJohnsonNorton.hpp b/ThermoFun/Substances/Solvent/Reaktoro/WaterElectroStateJohnsonNorton.hpp index e731cf77..d8b25519 100644 --- a/ThermoFun/Substances/Solvent/Reaktoro/WaterElectroStateJohnsonNorton.hpp +++ b/ThermoFun/Substances/Solvent/Reaktoro/WaterElectroStateJohnsonNorton.hpp @@ -27,6 +27,6 @@ struct WaterElectroState; struct WaterThermoState; // Calculate the electrostatic state of water using the model of Johnson and Norton (1991) -auto waterElectroStateJohnsonNorton(Reaktoro_::Temperature T, const WaterThermoState& wts) -> WaterElectroState; +auto waterElectroStateJohnsonNorton(Reaktoro_::Temperature T, const WaterThermoState& wts, int state=-1) -> WaterElectroState; } // namespace Reaktoro diff --git a/ThermoFun/Substances/Solvent/WaterElectroFernandez1997.cpp b/ThermoFun/Substances/Solvent/WaterElectroFernandez1997.cpp index 8e7de63f..a3ede01e 100644 --- a/ThermoFun/Substances/Solvent/WaterElectroFernandez1997.cpp +++ b/ThermoFun/Substances/Solvent/WaterElectroFernandez1997.cpp @@ -42,7 +42,7 @@ auto epsilonF (Reaktoro_::Temperature T, Reaktoro_::ThermoScalar rho) -> Reaktor } -auto electroPropertiesWaterFernandez1997(/*PropertiesSolvent ps,*/ Reaktoro_::Temperature TC, Reaktoro_::Pressure Pbar, Substance substance) -> ElectroPropertiesSolvent +auto electroPropertiesWaterFernandez1997(/*PropertiesSolvent ps,*/ Reaktoro_::Temperature TC, Reaktoro_::Pressure Pbar, Substance substance, int state) -> ElectroPropertiesSolvent { ElectroPropertiesSolvent wep; @@ -52,7 +52,7 @@ auto electroPropertiesWaterFernandez1997(/*PropertiesSolvent ps,*/ Reaktoro_::Te Reaktoro_::Temperature TK(TC.val + 273.15); Reaktoro_::Pressure P(Pbar.val * 1e05); - auto psol = th.propertiesSolvent(TK.val, P.val, substance.symbol()); + auto psol = th.propertiesSolvent(TK.val, P.val, substance.symbol(), state); const auto RHO = psol.density ; const auto eps = epsilonF(TK, RHO); @@ -61,12 +61,12 @@ auto electroPropertiesWaterFernandez1997(/*PropertiesSolvent ps,*/ Reaktoro_::Te // numerical approximation of epsilonT and epsilonTT Reaktoro_::Temperature T_plus (TC.val+TC.val*0.001); Reaktoro_::Temperature TK_plus (T_plus.val + 273.15); - auto RHO_plus = th.propertiesSolvent(TK_plus.val, P.val, substance.symbol()).density; + auto RHO_plus = th.propertiesSolvent(TK_plus.val, P.val, substance.symbol(), state).density; auto eps_plus = epsilonF(TK_plus, RHO_plus); Reaktoro_::Temperature T_minus (TC.val-TC.val*0.001); Reaktoro_::Temperature TK_minus (T_minus.val + 273.15); - auto RHO_minus = th.propertiesSolvent(TK_minus.val, P.val, substance.symbol()).density; + auto RHO_minus = th.propertiesSolvent(TK_minus.val, P.val, substance.symbol(), state).density; auto eps_minus = epsilonF(TK_minus, RHO_minus); const auto epsilonT = (eps_plus - eps_minus) / ((TK_plus-TK_minus)); @@ -75,12 +75,12 @@ auto electroPropertiesWaterFernandez1997(/*PropertiesSolvent ps,*/ Reaktoro_::Te // numerical approximation of epsilonP and epsilonPP Reaktoro_::Pressure P_plus (Pbar.val+Pbar.val*0.001); Reaktoro_::Pressure P_plusPa (P_plus.val*1e05); - RHO_plus = th.propertiesSolvent(TK.val, P_plusPa.val, substance.symbol()).density; + RHO_plus = th.propertiesSolvent(TK.val, P_plusPa.val, substance.symbol(), state).density; eps_plus = epsilonF(TK, RHO_plus); Reaktoro_::Pressure P_minus (Pbar.val-Pbar.val*0.001); Reaktoro_::Pressure P_minusPa (P_minus.val*1e05); - RHO_minus = th.propertiesSolvent(TK.val, P_minusPa.val, substance.symbol()).density; + RHO_minus = th.propertiesSolvent(TK.val, P_minusPa.val, substance.symbol(), state).density; eps_minus = epsilonF(TK, RHO_minus); const auto epsilonP = (eps_plus - eps_minus) / ((P_plus-P_minus)); diff --git a/ThermoFun/Substances/Solvent/WaterElectroFernandez1997.h b/ThermoFun/Substances/Solvent/WaterElectroFernandez1997.h index c63de3fb..1acd8ea2 100644 --- a/ThermoFun/Substances/Solvent/WaterElectroFernandez1997.h +++ b/ThermoFun/Substances/Solvent/WaterElectroFernandez1997.h @@ -8,7 +8,7 @@ namespace ThermoFun { class Substance; struct ElectroPropertiesSolvent; -auto electroPropertiesWaterFernandez1997(Reaktoro_::Temperature TC, Reaktoro_::Pressure Pbar, Substance substance) -> ElectroPropertiesSolvent; +auto electroPropertiesWaterFernandez1997(Reaktoro_::Temperature TC, Reaktoro_::Pressure Pbar, Substance substance, int state =-1) -> ElectroPropertiesSolvent; } diff --git a/ThermoFun/Substances/Solvent/WaterElectroSverjensky2014.cpp b/ThermoFun/Substances/Solvent/WaterElectroSverjensky2014.cpp index 889fbb9a..e5a907b0 100644 --- a/ThermoFun/Substances/Solvent/WaterElectroSverjensky2014.cpp +++ b/ThermoFun/Substances/Solvent/WaterElectroSverjensky2014.cpp @@ -27,7 +27,7 @@ auto epsilonS(Reaktoro_::Temperature T, Reaktoro_::ThermoScalar RHO) -> Reaktoro return exp(b[1]*T + b[2]*pow(T,0.5) + b[3])*pow(RHO,(a[1]*T + a[2]*pow(T,0.5) + a[3])); } -auto electroPropertiesWaterSverjensky2014(/*PropertiesSolvent ps,*/ Reaktoro_::Temperature TC, Reaktoro_::Pressure Pbar, Substance substance) -> ElectroPropertiesSolvent +auto electroPropertiesWaterSverjensky2014(/*PropertiesSolvent ps,*/ Reaktoro_::Temperature TC, Reaktoro_::Pressure Pbar, Substance substance, int state) -> ElectroPropertiesSolvent { ElectroPropertiesSolvent wep; @@ -37,7 +37,7 @@ auto electroPropertiesWaterSverjensky2014(/*PropertiesSolvent ps,*/ Reaktoro_::T Reaktoro_::Temperature TK(TC.val + 273.15); Reaktoro_::Pressure P(Pbar.val * 1e05); - auto psol = th.propertiesSolvent(TK.val, P.val, substance.symbol()); + auto psol = th.propertiesSolvent(TK.val, P.val, substance.symbol(), state); const auto RHO = psol.density /1000; const auto eps = epsilonS(TC, RHO); @@ -45,11 +45,11 @@ auto electroPropertiesWaterSverjensky2014(/*PropertiesSolvent ps,*/ Reaktoro_::T // numerical approximation of epsilonT and epsilonTT Reaktoro_::Temperature T_plus (TC.val+TC.val*0.001); - auto RHO_plus = th.propertiesSolvent(T_plus.val+273.15, P.val, substance.symbol()).density / 1000; + auto RHO_plus = th.propertiesSolvent(T_plus.val+273.15, P.val, substance.symbol(), state).density / 1000; auto eps_plus = epsilonS(T_plus, RHO_plus); Reaktoro_::Temperature T_minus (TC.val-TC.val*0.001); - auto RHO_minus = th.propertiesSolvent(T_minus.val+273.15, P.val, substance.symbol()).density / 1000; + auto RHO_minus = th.propertiesSolvent(T_minus.val+273.15, P.val, substance.symbol(), state).density / 1000; auto eps_minus = epsilonS(T_minus, RHO_minus); const auto epsilonT = (eps_plus - eps_minus) / ((T_plus-T_minus)); @@ -58,12 +58,12 @@ auto electroPropertiesWaterSverjensky2014(/*PropertiesSolvent ps,*/ Reaktoro_::T // numerical approximation of epsilonP and epsilonPP Reaktoro_::Pressure P_plus (Pbar.val+Pbar.val*0.001); Reaktoro_::Pressure P_plusPa (P_plus.val*1e05); - RHO_plus = th.propertiesSolvent(TK.val, P_plusPa.val, substance.symbol()).density / 1000; + RHO_plus = th.propertiesSolvent(TK.val, P_plusPa.val, substance.symbol(), state).density / 1000; eps_plus = epsilonS(TC, RHO_plus); Reaktoro_::Pressure P_minus (Pbar.val-Pbar.val*0.001); Reaktoro_::Pressure P_minusPa (P_minus.val*1e05); - RHO_minus = th.propertiesSolvent(TK.val, P_minusPa.val, substance.symbol()).density / 1000; + RHO_minus = th.propertiesSolvent(TK.val, P_minusPa.val, substance.symbol(), state).density / 1000; eps_minus = epsilonS(TC, RHO_minus); const auto epsilonP = (eps_plus - eps_minus) / ((P_plus-P_minus)); diff --git a/ThermoFun/Substances/Solvent/WaterElectroSverjensky2014.h b/ThermoFun/Substances/Solvent/WaterElectroSverjensky2014.h index e061f554..258e6d29 100644 --- a/ThermoFun/Substances/Solvent/WaterElectroSverjensky2014.h +++ b/ThermoFun/Substances/Solvent/WaterElectroSverjensky2014.h @@ -8,7 +8,7 @@ namespace ThermoFun { class Substance; struct ElectroPropertiesSolvent; -auto electroPropertiesWaterSverjensky2014(Reaktoro_::Temperature TC, Reaktoro_::Pressure Pbar, Substance substance) -> ElectroPropertiesSolvent; +auto electroPropertiesWaterSverjensky2014(Reaktoro_::Temperature TC, Reaktoro_::Pressure Pbar, Substance substance, int state = -1) -> ElectroPropertiesSolvent; } diff --git a/ThermoFun/Substances/Solvent/WaterHGK-JNgems.cpp b/ThermoFun/Substances/Solvent/WaterHGK-JNgems.cpp index e28f2f25..80ace0e5 100644 --- a/ThermoFun/Substances/Solvent/WaterHGK-JNgems.cpp +++ b/ThermoFun/Substances/Solvent/WaterHGK-JNgems.cpp @@ -390,7 +390,7 @@ auto WaterHGKgems::calculateWaterHGKgems(double T, double &P) -> void if( T < 0.01 && T >= 0.0 ) // Deg. C! T = 0.01; - if( P < 6.11732e-3 ) // 6.11732e-3 is P_sat at triple point of H2O + if( P < 6.11732e-3 ) // 6.11732e-3 is P_sat at triple point of H2O in bar // At lower pressures HKF/HGK runs unstable or crashes P = 0.0; // 06.12.2006 DK @@ -409,7 +409,7 @@ auto WaterHGKgems::calculateWaterHGKgems(double T, double &P) -> void auto diff = fabs(fabs(P) - fabs(aSta.Psat)); if( (fabs( P ) == 0) || (diff < 1e-14) ) - { // set only T + { // set only T, for Psat aSpc.isat=1; aSpc.iopt=1; aSpc.metastable = 0; @@ -483,7 +483,10 @@ auto WaterHGKgems::calculateWaterHGKgems(double T, double &P) -> void WPROPS tw = wl; memcpy(&wl, &wr, sizeof(WPROPS)); - memcpy(&wr, &tw, sizeof(WPROPS)); + if (aSta.Temp <= cr->Tc) // only below critical T we have two phases + memcpy(&wr, &tw, sizeof(WPROPS)); + else + aSta.Dens[1] = aSta.Dens[0]; //double temp = aSta.Dens[0]; //aSta.Dens[0] = aSta.Dens[1]; diff --git a/ThermoFun/ThermoEngine.cpp b/ThermoFun/ThermoEngine.cpp index d18d5932..4963a56c 100644 --- a/ThermoFun/ThermoEngine.cpp +++ b/ThermoFun/ThermoEngine.cpp @@ -14,6 +14,7 @@ #include "OptimizationUtils.h" #include +#include namespace ThermoFun { @@ -58,10 +59,10 @@ using ThermoPropertiesSubstanceFunction = std::function; using ElectroPropertiesSolventFunction = - std::function; + std::function; using PropertiesSolventFunction = - std::function; + std::function; using ThermoPropertiesReactionFunction = std::function; @@ -105,15 +106,15 @@ struct ThermoEngine::Impl }; thermo_properties_substance_fn = memoize(thermo_properties_substance_fn); - electro_properties_solvent_fn = [=](double T, double P_, double &P, std::string symbol) { + electro_properties_solvent_fn = [=](double T, double P_, double &P, std::string symbol, int state) { auto x = P_; - return electroPropertiesSolvent(T, P, symbol); + return electroPropertiesSolvent(T, P, symbol, state); }; electro_properties_solvent_fn = memoize(electro_properties_solvent_fn); - properties_solvent_fn = [=](double T, double P_, double &P, std::string symbol) { + properties_solvent_fn = [=](double T, double P_, double &P, std::string symbol, int state) { auto x = P_; - return propertiesSolvent(T, P, symbol); + return propertiesSolvent(T, P, symbol, state); }; properties_solvent_fn = memoize(properties_solvent_fn); @@ -248,26 +249,26 @@ struct ThermoEngine::Impl } case MethodGenEoS_Thrift::type::CTPM_HKF: { - tps = SoluteHKFgems(pref.workSubstance).thermoProperties(T, P, properties_solvent_fn(T, P, P, solventSymbol), electro_properties_solvent_fn(T, P, P, solventSymbol)); + tps = SoluteHKFgems(pref.workSubstance).thermoProperties(T, P, properties_solvent_fn(T, P, P, solventSymbol, -1), electro_properties_solvent_fn(T, P, P, solventSymbol, -1)); break; } case MethodGenEoS_Thrift::type::CTPM_HKFR: { - tps = SoluteHKFreaktoro(pref.workSubstance).thermoProperties(T, P, properties_solvent_fn(T, P, P, solventSymbol), electro_properties_solvent_fn(T, P, P, solventSymbol)); + tps = SoluteHKFreaktoro(pref.workSubstance).thermoProperties(T, P, properties_solvent_fn(T, P, P, solventSymbol, -1), electro_properties_solvent_fn(T, P, P, solventSymbol, -1)); break; } case MethodGenEoS_Thrift::type::CTPM_HP98: { double Pr = database.getSubstance(solventSymbol).referenceP(); double Tr = database.getSubstance(solventSymbol).referenceT(); - tps = SoluteHollandPowell98(pref.workSubstance).thermoProperties(T, P, properties_solvent_fn(Tr, Pr, Pr, solventSymbol), properties_solvent_fn(T, P, P, solventSymbol)); + tps = SoluteHollandPowell98(pref.workSubstance).thermoProperties(T, P, properties_solvent_fn(Tr, Pr, Pr, solventSymbol, -1), properties_solvent_fn(T, P, P, solventSymbol, -1)); break; } case MethodGenEoS_Thrift::type::CTPM_AN91: { double Pr = database.getSubstance(solventSymbol).referenceP(); double Tr = database.getSubstance(solventSymbol).referenceT(); - tps = SoluteAnderson91(pref.workSubstance).thermoProperties(T, P, properties_solvent_fn(Tr, Pr, Pr, solventSymbol), properties_solvent_fn(T, P, P, solventSymbol)); + tps = SoluteAnderson91(pref.workSubstance).thermoProperties(T, P, properties_solvent_fn(Tr, Pr, Pr, solventSymbol, -1), properties_solvent_fn(T, P, P, solventSymbol, -1)); break; } // default: @@ -302,7 +303,7 @@ struct ThermoEngine::Impl { double Pr = database.getSubstance(solventSymbol).referenceP(); double Tr = database.getSubstance(solventSymbol).referenceT(); - tps = SoluteAkinfievDiamondEOS(pref.workSubstance).thermoProperties(T, P, tps, thermo_properties_substance_fn(T, P, P, solventSymbol), WaterIdealGasWoolley(database.getSubstance(solventSymbol)).thermoProperties(T, P), properties_solvent_fn(T, P, P, solventSymbol), thermo_properties_substance_fn(Tr, Pr, Pr, solventSymbol), WaterIdealGasWoolley(database.getSubstance(solventSymbol)).thermoProperties(Tr, Pr), properties_solvent_fn(Tr, Pr, Pr, solventSymbol)); + tps = SoluteAkinfievDiamondEOS(pref.workSubstance).thermoProperties(T, P, tps, thermo_properties_substance_fn(T, P, P, solventSymbol), WaterIdealGasWoolley(database.getSubstance(solventSymbol)).thermoProperties(T, P), properties_solvent_fn(T, P, P, solventSymbol, -1), thermo_properties_substance_fn(Tr, Pr, Pr, solventSymbol), WaterIdealGasWoolley(database.getSubstance(solventSymbol)).thermoProperties(Tr, Pr), properties_solvent_fn(Tr, Pr, Pr, solventSymbol, -1)); break; } case MethodCorrP_Thrift::type::CPM_CEH: @@ -433,15 +434,17 @@ struct ThermoEngine::Impl return tps; } - auto electroPropertiesSolvent(double T, double &P, std::string solvent) const -> ElectroPropertiesSolvent + auto electroPropertiesSolvent(double T, double &P, std::string solvent, int state) const -> ElectroPropertiesSolvent { - return electroPropertiesSolvent(T, P, database.getSubstance(solvent)); + return electroPropertiesSolvent(T, P, database.getSubstance(solvent), state); } - auto electroPropertiesSolvent(double T, double &P, const Substance &solvent) const -> ElectroPropertiesSolvent + auto electroPropertiesSolvent(double T, double &P, const Substance &solvent, int state) const -> ElectroPropertiesSolvent { ThermoPreferences pref = getThermoPreferencesSubstance(solvent); - PropertiesSolvent ps = propertiesSolvent(T, P, solvent); + if (state>-1 && state<2) + pref.solventState = state; + PropertiesSolvent ps = propertiesSolvent(T, P, solvent, pref.solventState); ElectroPropertiesSolvent eps; if (pref.isH2OSolvent) @@ -450,22 +453,22 @@ struct ThermoEngine::Impl { case MethodGenEoS_Thrift::type::CTPM_WJNR: { - eps = WaterJNreaktoro(pref.workSubstance).electroPropertiesSolvent(T, P, ps); + eps = WaterJNreaktoro(pref.workSubstance).electroPropertiesSolvent(T, P, ps, pref.solventState); break; } case MethodGenEoS_Thrift::type::CTPM_WJNG: { - eps = WaterJNgems(pref.workSubstance).electroPropertiesSolvent(T, P /*, ps*/); + eps = WaterJNgems(pref.workSubstance).electroPropertiesSolvent(T, P /*, ps*/, pref.solventState); break; } case MethodGenEoS_Thrift::type::CTPM_WSV14: { - eps = WaterElectroSverjensky2014(pref.workSubstance).electroPropertiesSolvent(T, P /*, ps*/); + eps = WaterElectroSverjensky2014(pref.workSubstance).electroPropertiesSolvent(T, P /*, ps*/, pref.solventState); break; } case MethodGenEoS_Thrift::type::CTPM_WF97: { - eps = WaterElectroFernandez1997(pref.workSubstance).electroPropertiesSolvent(T, P /*, ps*/); + eps = WaterElectroFernandez1997(pref.workSubstance).electroPropertiesSolvent(T, P /*, ps*/, pref.solventState); break; } // default: @@ -476,14 +479,16 @@ struct ThermoEngine::Impl return eps; } - auto propertiesSolvent(double T, double &P, std::string solvent) const -> PropertiesSolvent + auto propertiesSolvent(double T, double &P, std::string solvent, int state) const -> PropertiesSolvent { - return propertiesSolvent(T, P, database.getSubstance(solvent)); + return propertiesSolvent(T, P, database.getSubstance(solvent), state); } - auto propertiesSolvent(double T, double &P, const Substance &solvent) const -> PropertiesSolvent + auto propertiesSolvent(double T, double &P, const Substance &solvent, int state) const -> PropertiesSolvent { ThermoPreferences pref = getThermoPreferencesSubstance(solvent); + if (state > -1 && state <2) + pref.solventState = state; PropertiesSolvent ps; if (pref.isH2OSolvent) @@ -625,17 +630,17 @@ struct ThermoEngine::Impl } case MethodCorrT_Thrift::type::CTM_DMD: // Dolejs-Maning 2010 density model { - tpr = ReactionDolejsManning10(pref.workReaction).thermoProperties(T, P, properties_solvent_fn(T, P, P, solventSymbol)); + tpr = ReactionDolejsManning10(pref.workReaction).thermoProperties(T, P, properties_solvent_fn(T, P, P, solventSymbol, -1)); break; } case MethodCorrT_Thrift::type::CTM_DKR: // Marshall-Franck density model { - tpr = ReactionFrantzMarshall(pref.workReaction).thermoProperties(T, P, properties_solvent_fn(T, P, P, solventSymbol)); + tpr = ReactionFrantzMarshall(pref.workReaction).thermoProperties(T, P, properties_solvent_fn(T, P, P, solventSymbol, -1)); break; } case MethodCorrT_Thrift::type::CTM_MRB: // Calling modified Ryzhenko-Bryzgalin model TW KD 08.2007 { - return ReactionRyzhenkoBryzgalin(pref.workReaction).thermoProperties(T, P, properties_solvent_fn(T, P, P, solventSymbol)); // NOT TESTED!!! + return ReactionRyzhenkoBryzgalin(pref.workReaction).thermoProperties(T, P, properties_solvent_fn(T, P, P, solventSymbol, -1)); // NOT TESTED!!! } case MethodCorrT_Thrift::type::CTM_IKZ: { @@ -776,14 +781,14 @@ auto ThermoEngine::thermoPropertiesSubstance(double T, double &P, std::string su return pimpl->thermo_properties_substance_fn(T, P, P, substance); } -auto ThermoEngine::electroPropertiesSolvent(double T, double &P, std::string solvent) const -> ElectroPropertiesSolvent +auto ThermoEngine::electroPropertiesSolvent(double T, double &P, std::string solvent, int state) const -> ElectroPropertiesSolvent { - return pimpl->electro_properties_solvent_fn(T, P, P, solvent); + return pimpl->electro_properties_solvent_fn(T, P, P, solvent, state); } -auto ThermoEngine::propertiesSolvent(double T, double &P, std::string solvent) const -> PropertiesSolvent +auto ThermoEngine::propertiesSolvent(double T, double &P, std::string solvent, int state) const -> PropertiesSolvent { - return pimpl->properties_solvent_fn(T, P, P, solvent); + return pimpl->properties_solvent_fn(T, P, P, solvent, state); } auto ThermoEngine::thermoPropertiesSubstance(double T, double &P, const Substance& substance) const -> ThermoPropertiesSubstance @@ -791,14 +796,14 @@ auto ThermoEngine::thermoPropertiesSubstance(double T, double &P, const Substanc return pimpl->thermoPropertiesSubstance(T, P, substance); } -auto ThermoEngine::electroPropertiesSolvent(double T, double &P, const Substance& solvent) const -> ElectroPropertiesSolvent +auto ThermoEngine::electroPropertiesSolvent(double T, double &P, const Substance& solvent, int state) const -> ElectroPropertiesSolvent { - return pimpl->electroPropertiesSolvent(T, P, solvent); + return pimpl->electroPropertiesSolvent(T, P, solvent, state); } -auto ThermoEngine::propertiesSolvent(double T, double &P, const Substance& solvent) const -> PropertiesSolvent +auto ThermoEngine::propertiesSolvent(double T, double &P, const Substance& solvent, int state) const -> PropertiesSolvent { - return pimpl->propertiesSolvent(T, P, solvent); + return pimpl->propertiesSolvent(T, P, solvent, state); } // Reaction diff --git a/ThermoFun/ThermoEngine.h b/ThermoFun/ThermoEngine.h index 6f644df8..e3e01827 100644 --- a/ThermoFun/ThermoEngine.h +++ b/ThermoFun/ThermoEngine.h @@ -85,25 +85,29 @@ class ThermoEngine /// @param T The temperature value (in units of K) /// @param P The pressure value (in units of Pa) /// @param solvent The symbol of the substance solvent - auto electroPropertiesSolvent(double T, double &P, std::string solvent) const -> ElectroPropertiesSolvent; + /// @param state The solvent state, 0 liquid, 1 gas/vapor, default -1 determined by the record key "aggregate_state" in the database "4": "AS_AQUEOUS" or "O": "AS_GAS" + auto electroPropertiesSolvent(double T, double &P, std::string solvent, int state=-1) const -> ElectroPropertiesSolvent; /// Calculate the electro-chemical properties of a substance. /// @param T The temperature value (in units of K) /// @param P The pressure value (in units of Pa) /// @param solvent substance solvent object - auto electroPropertiesSolvent(double T, double &P, const Substance& solvent) const -> ElectroPropertiesSolvent; + /// @param state The solvent state, 0 liquid, 1 gas/vapor, default -1 determined by the record key "aggregate_state" in the database "4": "AS_AQUEOUS" or "O": "AS_GAS" + auto electroPropertiesSolvent(double T, double &P, const Substance& solvent, int state=-1) const -> ElectroPropertiesSolvent; /// Calculate the physical properties of a substance. /// @param T The temperature value (in units of K) /// @param P The pressure value (in units of Pa) /// @param solvent The symbol of the substance solvent - auto propertiesSolvent(double T, double &P, std::string solvent) const -> PropertiesSolvent; + /// @param state The solvent state, 0 liquid, 1 gas/vapor, default -1 determined by the record key "aggregate_state" in the database "4": "AS_AQUEOUS" or "O": "AS_GAS" + auto propertiesSolvent(double T, double &P, std::string solvent, int state=-1) const -> PropertiesSolvent; /// Calculate the physical properties of a substance. /// @param T The temperature value (in units of K) /// @param P The pressure value (in units of Pa) /// @param solvent substance solvent object - auto propertiesSolvent(double T, double &P, const Substance& solvent) const -> PropertiesSolvent; + /// @param state The solvent state, 0 liquid, 1 gas/vapor, default -1 determined by the record key "aggregate_state" in the database "4": "AS_AQUEOUS" or "O": "AS_GAS" + auto propertiesSolvent(double T, double &P, const Substance& solvent, int state=-1) const -> PropertiesSolvent; // Reaction /// Calculate the thermodynamic properties of a reaction. diff --git a/ThermoFun/ThermoModelsSolvent.h b/ThermoFun/ThermoModelsSolvent.h index 7f7f2d84..542982ec 100644 --- a/ThermoFun/ThermoModelsSolvent.h +++ b/ThermoFun/ThermoModelsSolvent.h @@ -96,7 +96,7 @@ class WaterWP95reaktoro /// @param state defines the state of the solvent (liquid 0, vapor 1, below the critical point) auto propertiesSolvent (double T, double &P, int state) -> PropertiesSolvent; - /// Return the themrodynamic properties of the solvent. + /// Return the thermodynamic properties of the solvent. /// @param T The temperature value (in units of K) /// @param P The pressure value (in units of Pa) /// @param state defines the state of the solvent (liquid 0, vapor 1, below the critical point) @@ -128,7 +128,7 @@ class WaterZhangDuan2005 /// @param state defines the state of the solvent (liquid 0, vapor 1, below the critical point) auto propertiesSolvent (double T, double P, int state) -> PropertiesSolvent; - /// Return the themrodynamic properties of the solvent. + /// Return the thermodynamic properties of the solvent. /// @param T The temperature value (in units of K) /// @param P The pressure value (in units of Pa) /// @param state defines the state of the solvent (liquid 0, vapor 1, below the critical point) diff --git a/pytests/Substances/Solvent/test_water_EoS.py b/pytests/Substances/Solvent/test_water_EoS.py new file mode 100644 index 00000000..1f68fcbb --- /dev/null +++ b/pytests/Substances/Solvent/test_water_EoS.py @@ -0,0 +1,60 @@ +import thermofun as thermofun +import pytest as pytest +import unittest + + +class TestWaterEoS(unittest.TestCase): + + def setUp(self): + self.engine = thermofun.ThermoEngine('pytests/Substances/Solvent/water-thermofun.json') + + def test_water_electro_properties(self): + P = 300e5 # region crit + T = 420+273.15 + assert self.engine.propertiesSolvent(T,P, "H2O@", 1).density.val == pytest.approx(203.47268444814958, 1e-5, 1e-14) + assert self.engine.electroPropertiesSolvent(T,P, "H2O@", 1).epsilon.val == pytest.approx(3.1942838457177976, 1e-5, 1e-14) + assert self.engine.propertiesSolvent(T,P, "H2O@", 0).density.val == pytest.approx(203.47268444814958, 1e-5, 1e-14) + assert self.engine.electroPropertiesSolvent(T,P, "H2O@", 0).epsilon.val == pytest.approx(3.1942838457177976, 1e-5, 1e-14) + + assert self.engine.propertiesSolvent(T,P, "H2O@reak", 1).density.val == pytest.approx(203.2389910219359, 1e-5, 1e-14) + assert self.engine.electroPropertiesSolvent(T,P, "H2O@reak", 1).epsilon.val == pytest.approx(3.190715033274149, 1e-5, 1e-14) + assert self.engine.propertiesSolvent(T,P, "H2O@reak", 0).density.val == pytest.approx(203.2389910219359, 1e-5, 1e-14) + assert self.engine.electroPropertiesSolvent(T,P, "H2O@reak", 0).epsilon.val == pytest.approx(3.190715033274149, 1e-5, 1e-14) + + P = 40e5 # vapour + T = 350+273.15 + assert self.engine.propertiesSolvent(T,P, "H2O@", 1).density.val == pytest.approx(15.05176552064709, 1e-5, 1e-14) + assert self.engine.electroPropertiesSolvent(T,P, "H2O@", 1).epsilon.val == pytest.approx(1.1120564097130872, 1e-5, 1e-14) + assert self.engine.propertiesSolvent(T,P, "H2O@", 0).density.val == pytest.approx(452.42695689302792, 1e-5, 1e-14) + assert self.engine.electroPropertiesSolvent(T,P, "H2O@", 0).epsilon.val == pytest.approx(9.4564825264140211, 1e-5, 1e-14) + + assert self.engine.propertiesSolvent(T,P, "H2O@reak", 1).density.val == pytest.approx(15.043737111107781, 1e-5, 1e-14) + assert self.engine.electroPropertiesSolvent(T,P, "H2O@reak", 1).epsilon.val == pytest.approx(1.1119933554710066, 1e-5, 1e-14) + assert self.engine.propertiesSolvent(T,P, "H2O@reak", 0).density.val == pytest.approx(15.043737111107781, 1e-5, 1e-14) + assert self.engine.electroPropertiesSolvent(T,P, "H2O@reak", 0).epsilon.val == pytest.approx(1.1119933554710066, 1e-5, 1e-14) + + # region 3 + P = 300e5 + T = 380+273.15 + assert self.engine.propertiesSolvent(T,P, "H2O@", 1).density.val == pytest.approx(533.9645612493881, 1e-5, 1e-14) + assert self.engine.electroPropertiesSolvent(T,P, "H2O@", 1).epsilon.val == pytest.approx(11.315069043656928, 1e-5, 1e-14) + assert self.engine.propertiesSolvent(T,P, "H2O@", 0).density.val == pytest.approx(533.9645612493881, 1e-5, 1e-14) + assert self.engine.electroPropertiesSolvent(T,P, "H2O@", 0).epsilon.val == pytest.approx(11.315069043656928, 1e-5, 1e-14) + + assert self.engine.propertiesSolvent(T,P, "H2O@reak", 1).density.val == pytest.approx(533.9301819281667, 1e-5, 1e-14) + assert self.engine.electroPropertiesSolvent(T,P, "H2O@reak", 1).epsilon.val == pytest.approx(11.31400838132288, 1e-5, 1e-14) + assert self.engine.propertiesSolvent(T,P, "H2O@reak", 0).density.val == pytest.approx(533.9301819281667, 1e-5, 1e-14) + assert self.engine.electroPropertiesSolvent(T,P, "H2O@reak", 0).epsilon.val == pytest.approx(11.31400838132288, 1e-5, 1e-14) + + # region 2 + P = 600e5 + T = 450+273.15 + assert self.engine.propertiesSolvent(T,P, "H2O@", 1).density.val == pytest.approx(479.86864831430904, 1e-5, 1e-14) + assert self.engine.electroPropertiesSolvent(T,P, "H2O@", 1).epsilon.val == pytest.approx(8.635464564374157, 1e-5, 1e-14) + assert self.engine.propertiesSolvent(T,P, "H2O@", 0).density.val == pytest.approx(479.86864831430904, 1e-5, 1e-14) + assert self.engine.electroPropertiesSolvent(T,P, "H2O@", 0).epsilon.val == pytest.approx(8.635464564374157, 1e-5, 1e-14) + + assert self.engine.propertiesSolvent(T,P, "H2O@reak", 1).density.val == pytest.approx(479.5064053309655, 1e-5, 1e-14) + assert self.engine.electroPropertiesSolvent(T,P, "H2O@reak", 1).epsilon.val == pytest.approx(8.626211060320887, 1e-5, 1e-14) + assert self.engine.propertiesSolvent(T,P, "H2O@reak", 0).density.val == pytest.approx(479.5064053309655, 1e-5, 1e-14) + assert self.engine.electroPropertiesSolvent(T,P, "H2O@reak", 0).epsilon.val == pytest.approx(8.626211060320887, 1e-5, 1e-14) \ No newline at end of file diff --git a/pytests/Substances/Solvent/water-thermofun.json b/pytests/Substances/Solvent/water-thermofun.json new file mode 100644 index 00000000..54cf8468 --- /dev/null +++ b/pytests/Substances/Solvent/water-thermofun.json @@ -0,0 +1,447 @@ +{ + "thermodataset": "aq17-gem", + "datasources": [ + "db.thermohub.org", + "[2017MIR/WAG]" + ], + "date": "10.12.2019 12:58:35", + "substances": [ + { + "name": "Al(OH)2+", + "symbol": "my_Al(OH)2+", + "formula": "Al(OH)2+", + "formula_charge": 1, + "aggregate_state": { + "4": "AS_AQUEOUS" + }, + "class_": { + "2": "SC_AQSOLUTE" + }, + "Tst": 298.15, + "Pst": 100000, + "TPMethods": [ + { + "method": { + "3": "solute_hkf88_reaktoro" + }, + "eos_hkf_coeffs": { + "values": [ + 0.24940000474453, + -169.08999633789, + 6.4145998954773, + -27091, + 16.743900299072, + -10465, + 53240, + 0 + ], + "units": [ + "cal/(mol*bar)", + "cal/mol", + "(cal*K)/mol", + "cal/(mol*K)", + "(cal*K)/mol", + "cal/mol" + ], + "names": [ + "a1", + "a2", + "a3", + "a4", + "c1", + "c2", + "wref" + ] + } + } + ], + "sm_heat_capacity_p": { + "values": [ + 40.865230560303 + ], + "errors": [ + 0 + ], + "units": [ + "J/(mol*K)" + ] + }, + "sm_gibbs_energy": { + "values": [ + -898292 + ], + "units": [ + "J/mol" + ] + }, + "sm_enthalpy": { + "values": [ + -995581 + ], + "units": [ + "J/mol" + ] + }, + "sm_entropy_abs": { + "values": [ + -50 + ], + "units": [ + "J/(mol*K)" + ] + }, + "sm_volume": { + "values": [ + 0.38507527112961 + ], + "errors": [ + 0 + ], + "units": [ + "J/bar" + ] + }, + "datasources": [ + "[2001TAG/SCH]", + "[2017MIR/WAG]" + ] + }, + { + "name": "Water HGK", + "symbol": "H2O@", + "formula": "H2O@", + "formula_charge": 0, + "aggregate_state": { + "4": "AS_AQUEOUS" + }, + "class_": { + "3": "SC_AQSOLVENT" + }, + "Tst": 298.15, + "Pst": 100000, + "TPMethods": [ + { + "method": { + "29": "water_eos_hgk84_lvs83_gems" + } + }, + { + "method": { + "25": "water_diel_jnort91_reaktoro" + } + } + ], + "sm_heat_capacity_p": { + "values": [ + 75.360527038574 + ], + "errors": [ + 0 + ], + "units": [ + "J/(mol*K)" + ] + }, + "sm_gibbs_energy": { + "values": [ + -237183 + ], + "errors": [ + 100 + ], + "units": [ + "J/mol" + ] + }, + "sm_enthalpy": { + "values": [ + -285881 + ], + "errors": [ + 200 + ], + "units": [ + "J/mol" + ] + }, + "sm_entropy_abs": { + "values": [ + 69.922996520996 + ], + "errors": [ + 0.10000000149012 + ], + "units": [ + "J/(mol*K)" + ] + }, + "sm_volume": { + "values": [ + 1.8068397045136 + ], + "errors": [ + 0 + ], + "units": [ + "J/bar" + ] + }, + "datasources": [ + "[1992JOH/OEL]" + ] + }, + { + "name": "Water HGK", + "symbol": "H2Og", + "formula": "H2O@", + "formula_charge": 0, + "aggregate_state": { + "0": "AS_GAS" + }, + "class_": { + "3": "SC_AQSOLVENT" + }, + "Tst": 298.15, + "Pst": 100000, + "TPMethods": [ + { + "method": { + "29": "water_eos_hgk84_lvs83_gems" + } + }, + { + "method": { + "26": "water_diel_jnort91_gems" + } + } + ], + "sm_heat_capacity_p": { + "values": [ + 75.360527038574 + ], + "errors": [ + 0 + ], + "units": [ + "J/(mol*K)" + ] + }, + "sm_gibbs_energy": { + "values": [ + -237183 + ], + "errors": [ + 100 + ], + "units": [ + "J/mol" + ] + }, + "sm_enthalpy": { + "values": [ + -285881 + ], + "errors": [ + 200 + ], + "units": [ + "J/mol" + ] + }, + "sm_entropy_abs": { + "values": [ + 69.922996520996 + ], + "errors": [ + 0.10000000149012 + ], + "units": [ + "J/(mol*K)" + ] + }, + "sm_volume": { + "values": [ + 1.8068397045136 + ], + "errors": [ + 0 + ], + "units": [ + "J/bar" + ] + }, + "datasources": [ + "[1992JOH/OEL]" + ] + }, + { + "name": "Water HGK", + "symbol": "H2O@", + "formula": "H2O@", + "formula_charge": 0, + "aggregate_state": { + "4": "AS_AQUEOUS" + }, + "class_": { + "3": "SC_AQSOLVENT" + }, + "Tst": 298.15, + "Pst": 100000, + "TPMethods": [ + { + "method": { + "29": "water_eos_hgk84_lvs83_gems" + } + }, + { + "method": { + "25": "water_diel_jnort91_reaktoro" + } + } + ], + "sm_heat_capacity_p": { + "values": [ + 75.360527038574 + ], + "errors": [ + 0 + ], + "units": [ + "J/(mol*K)" + ] + }, + "sm_gibbs_energy": { + "values": [ + -237183 + ], + "errors": [ + 100 + ], + "units": [ + "J/mol" + ] + }, + "sm_enthalpy": { + "values": [ + -285881 + ], + "errors": [ + 200 + ], + "units": [ + "J/mol" + ] + }, + "sm_entropy_abs": { + "values": [ + 69.922996520996 + ], + "errors": [ + 0.10000000149012 + ], + "units": [ + "J/(mol*K)" + ] + }, + "sm_volume": { + "values": [ + 1.8068397045136 + ], + "errors": [ + 0 + ], + "units": [ + "J/bar" + ] + }, + "datasources": [ + "[1992JOH/OEL]" + ] + }, + { + "name": "Water HGK", + "symbol": "H2O@reak", + "formula": "H2O@", + "formula_charge": 0, + "aggregate_state": { + "4": "AS_AQUEOUS" + }, + "class_": { + "3": "SC_AQSOLVENT" + }, + "Tst": 298.15, + "Pst": 100000, + "TPMethods": [ + { + "method": { + "32": "water_eos_iapws95_reaktoro" + } + }, + { + "method": { + "25": "water_diel_jnort91_reaktoro" + } + } + ], + "sm_heat_capacity_p": { + "values": [ + 75.360527038574 + ], + "errors": [ + 0 + ], + "units": [ + "J/(mol*K)" + ] + }, + "sm_gibbs_energy": { + "values": [ + -237183 + ], + "errors": [ + 100 + ], + "units": [ + "J/mol" + ] + }, + "sm_enthalpy": { + "values": [ + -285881 + ], + "errors": [ + 200 + ], + "units": [ + "J/mol" + ] + }, + "sm_entropy_abs": { + "values": [ + 69.922996520996 + ], + "errors": [ + 0.10000000149012 + ], + "units": [ + "J/(mol*K)" + ] + }, + "sm_volume": { + "values": [ + 1.8068397045136 + ], + "errors": [ + 0 + ], + "units": [ + "J/bar" + ] + }, + "datasources": [ + "[1992JOH/OEL]" + ] + } + ], + "reactions": [] +} diff --git a/python/pyThermoFun/pyThermoEngine.cpp b/python/pyThermoFun/pyThermoEngine.cpp index ddd21b30..b8251cae 100644 --- a/python/pyThermoFun/pyThermoEngine.cpp +++ b/python/pyThermoFun/pyThermoEngine.cpp @@ -39,15 +39,48 @@ void exportThermoEngine(py::module& m) auto appendData1 = static_cast(&ThermoEngine::appendData); auto appendData2 = static_cast, std::string)>(&ThermoEngine::appendData); + // Using lambda functions to provide default value for `state` + auto electroPropertiesSolvent1 = [](const ThermoEngine& self, double param1, double& param2, const std::string& param3, int state = -1) { + return self.electroPropertiesSolvent(param1, param2, param3, state); + }; + + auto electroPropertiesSolvent1_default = [](const ThermoEngine& self, double param1, double& param2, const std::string& param3) { + return self.electroPropertiesSolvent(param1, param2, param3, -1); + }; + + auto electroPropertiesSolvent2 = [](const ThermoEngine& self, double param1, double& param2, const Substance& param3, int state = -1) { + return self.electroPropertiesSolvent(param1, param2, param3, state); + }; + + auto electroPropertiesSolvent2_default = [](const ThermoEngine& self, double param1, double& param2, const Substance& param3) { + return self.electroPropertiesSolvent(param1, param2, param3, -1); + }; + + auto propertiesSolvent1 = [](const ThermoEngine& self, double param1, double& param2, const std::string& param3, int state = -1) { + return self.propertiesSolvent(param1, param2, param3, state); + }; + + auto propertiesSolvent1_default = [](const ThermoEngine& self, double param1, double& param2, const std::string& param3) { + return self.propertiesSolvent(param1, param2, param3, -1); + }; + + auto propertiesSolvent2 = [](const ThermoEngine& self, double param1, double& param2, const Substance& param3, int state = -1) { + return self.propertiesSolvent(param1, param2, param3, state); + }; + + auto propertiesSolvent2_default = [](const ThermoEngine& self, double param1, double& param2, const Substance& param3) { + return self.propertiesSolvent(param1, param2, param3, -1); + }; + auto thermoPropertiesSubstance1 = static_cast(&ThermoEngine::thermoPropertiesSubstance); - auto electroPropertiesSolvent1 = static_cast(&ThermoEngine::electroPropertiesSolvent); - auto propertiesSolvent1 = static_cast(&ThermoEngine::propertiesSolvent); + //auto electroPropertiesSolvent1 = static_cast(&ThermoEngine::electroPropertiesSolvent); + //auto propertiesSolvent1 = static_cast(&ThermoEngine::propertiesSolvent); auto thermoPropertiesReaction1 = static_cast(&ThermoEngine::thermoPropertiesReaction); auto thermoPropertiesReactionFromReactants1 = static_cast(&ThermoEngine::thermoPropertiesReactionFromReactants); auto thermoPropertiesSubstance2 = static_cast(&ThermoEngine::thermoPropertiesSubstance); - auto electroPropertiesSolvent2 = static_cast(&ThermoEngine::electroPropertiesSolvent); - auto propertiesSolvent2 = static_cast(&ThermoEngine::propertiesSolvent); + //auto electroPropertiesSolvent2 = static_cast(&ThermoEngine::electroPropertiesSolvent); + //auto propertiesSolvent2 = static_cast(&ThermoEngine::propertiesSolvent); auto thermoPropertiesReaction2 = static_cast(&ThermoEngine::thermoPropertiesReaction); auto thermoPropertiesReactionFromReactants2 = static_cast(&ThermoEngine::thermoPropertiesReactionFromReactants); @@ -61,13 +94,17 @@ void exportThermoEngine(py::module& m) .def("setSolventSymbol", &ThermoEngine::setSolventSymbol, "Sets the symbol of the solvent record present in the thermodynamic dataset. Will be used to calculate the solvent properties ", py::arg("symbol")) .def("solventSymbol", &ThermoEngine::solventSymbol, "Returns the symbol of the solvent record used to calculate the solvent properties") .def("thermoPropertiesSubstance", thermoPropertiesSubstance1, "Calculate the thermodynamic properties of a substance with a given symbol.", py::arg("temperature"), py::arg("pressure"), py::arg("symbol")) - .def("electroPropertiesSolvent", electroPropertiesSolvent1, "Calculate the electro-chemical properties of a substance solvent with a given symbol.", py::arg("temperature"), py::arg("pressure"), py::arg("symbol")) - .def("propertiesSolvent", propertiesSolvent1, "Calculate the properties of a substance solvent with a given symbol.", py::arg("temperature"), py::arg("pressure"), py::arg("symbol")) + .def("electroPropertiesSolvent", electroPropertiesSolvent1, "Calculate the electro-chemical properties of a substance solvent with a given symbol.", py::arg("temperature"), py::arg("pressure"), py::arg("symbol"), py::arg("state")) + .def("propertiesSolvent", propertiesSolvent1, "Calculate the properties of a substance solvent with a given symbol.", py::arg("temperature"), py::arg("pressure"), py::arg("symbol"), py::arg("state")) + .def("electroPropertiesSolvent", electroPropertiesSolvent1_default, "Calculate the electro-chemical properties of a substance solvent with a given symbol.", py::arg("temperature"), py::arg("pressure"), py::arg("symbol")) + .def("propertiesSolvent", propertiesSolvent1_default, "Calculate the properties of a substance solvent with a given symbol.", py::arg("temperature"), py::arg("pressure"), py::arg("symbol")) .def("thermoPropertiesReaction", thermoPropertiesReaction1, "Calculate the thermodynamic properties of a reaction with a given symbol, or for a given reaction equation.", py::arg("temperature"), py::arg("pressure"), py::arg("symbol")) .def("thermoPropertiesReactionFromReactants", thermoPropertiesReactionFromReactants1, "Calculate the thermodynamic properties of a reaction from reactants with a given symbol.", py::arg("temperature"), py::arg("pressure"), py::arg("symbol")) .def("thermoPropertiesSubstance", thermoPropertiesSubstance2, "Calculate the thermodynamic properties of a given substance object.", py::arg("temperature"), py::arg("pressure"), py::arg("substance")) - .def("electroPropertiesSolvent", electroPropertiesSolvent2, "Calculate the electro-chemical properties of a given substance solvent object.", py::arg("temperature"), py::arg("pressure"), py::arg("solvent")) - .def("propertiesSolvent", propertiesSolvent2, "Calculate the properties of a given substance solvent object.", py::arg("temperature"), py::arg("pressure"), py::arg("solvent")) + .def("electroPropertiesSolvent", electroPropertiesSolvent2, "Calculate the electro-chemical properties of a given substance solvent object.", py::arg("temperature"), py::arg("pressure"), py::arg("solvent"), py::arg("state")) + .def("propertiesSolvent", propertiesSolvent2, "Calculate the properties of a given substance solvent object.", py::arg("temperature"), py::arg("pressure"), py::arg("solvent"), py::arg("state")) + .def("electroPropertiesSolvent", electroPropertiesSolvent2_default, "Calculate the electro-chemical properties of a given substance solvent object.", py::arg("temperature"), py::arg("pressure"), py::arg("solvent")) + .def("propertiesSolvent", propertiesSolvent2_default, "Calculate the properties of a given substance solvent object.", py::arg("temperature"), py::arg("pressure"), py::arg("solvent")) .def("thermoPropertiesReaction", thermoPropertiesReaction2, "Calculate the thermodynamic properties of a given reaction object.", py::arg("temperature"), py::arg("pressure"), py::arg("reaction")) .def("thermoPropertiesReactionFromReactants", thermoPropertiesReactionFromReactants2, "Calculate the thermodynamic properties from the reactants of a given reaction object.", py::arg("temperature"), py::arg("pressure"), py::arg("reaction")) ; diff --git a/tests/interfaceTest/src/main.cpp b/tests/interfaceTest/src/main.cpp index 242728c3..26299c72 100644 --- a/tests/interfaceTest/src/main.cpp +++ b/tests/interfaceTest/src/main.cpp @@ -2,6 +2,9 @@ #include "ChemicalFun/FormulaParser.h" #include "GlobalVariables.h" //#include "ThermoFun/Common/ThermoScalar.hpp" +#include + + using namespace std; using namespace ThermoFun; @@ -58,8 +61,33 @@ {"gibbs_energy","entropy", "volume", "enthalpy"} // list of properties ).toCSVPropertyGrid("grid.csv"); */ // output + // Water solvent + + ThermoEngine engine("aq17-thermofun.json"); + + engine.appendData("water-thermofun.json"); + + double P = 300e5; + double T = 380+273.15; + + auto rho_solvent = engine.propertiesSolvent(T,P, "H2O@", 1).density.val; + + auto eplsilon_solvent = engine.electroPropertiesSolvent(T,P, "H2O@", 1).epsilon.val; + + auto rho_solvent2 = engine.propertiesSolvent(T,P, "H2O@", 0).density.val; - // Test + auto eplsilon_solvent2 = engine.electroPropertiesSolvent(T,P, "H2O@", 0).epsilon.val; + + auto rho_solvent3 = engine.propertiesSolvent(T,P, "H2O@").density.val; + + auto eplsilon_solvent3 = engine.electroPropertiesSolvent(T,P, "H2O@").epsilon.val; + + std::cout <<"end"<< std::endl; + + + + // Test mines add string + /* Database db("mines16-thermofun.json"); auto test = db.getSubstance("O2"); @@ -249,6 +277,8 @@ // batch.thermoPropertiesReaction(25, 1, "H2O@ = H+ + OH-", "logKr").toCSV("test_reac.cvs"); // batch.thermoPropertiesReaction(25, 1, "Al+3 + 4H2O@ + 0Ca+2 = 1Al(OH)4- + 4 \nH+", "logKr").toCSV("test_reac2.cvs"); +*/ + return 0; }