From 0710fef5f8c7fad8fae0fba666554c3afe399a1a Mon Sep 17 00:00:00 2001 From: jeremy karpowski Date: Mon, 3 Feb 2025 09:06:30 +0100 Subject: [PATCH 1/6] added thermo update --- include/cantera/thermo/IdealGasPhase.h | 2 +- include/cantera/thermo/PengRobinson.h | 55 +++++++++++++++++++++++++- include/cantera/thermo/ThermoPhase.h | 2 +- src/thermo/MixtureFugacityTP.cpp | 17 ++++++-- src/thermo/PengRobinson.cpp | 46 +++++++++++++++++++++ 5 files changed, 115 insertions(+), 7 deletions(-) diff --git a/include/cantera/thermo/IdealGasPhase.h b/include/cantera/thermo/IdealGasPhase.h index dcd551ccea1..674cf9a320d 100644 --- a/include/cantera/thermo/IdealGasPhase.h +++ b/include/cantera/thermo/IdealGasPhase.h @@ -366,7 +366,7 @@ class IdealGasPhase: public ThermoPhase * @param p Pressure (Pa) * @since New in %Cantera 3.0. */ - void setState_DP(double rho, double p) override { + void setState_DP(double rho, double p, double Tguess=300) override { if (p <= 0) { throw CanteraError("IdealGasPhase::setState_DP", "pressure must be positive"); diff --git a/include/cantera/thermo/PengRobinson.h b/include/cantera/thermo/PengRobinson.h index 8fda6dd19d1..03b877136ad 100644 --- a/include/cantera/thermo/PengRobinson.h +++ b/include/cantera/thermo/PengRobinson.h @@ -203,6 +203,7 @@ class PengRobinson : public MixtureFugacityTP * These are stored internally. */ double daAlpha_dT() const; + double daAlpha_dT(double& temp) const; //! Calculate second derivative @f$ d^2(a \alpha)/dT^2 @f$ /*! @@ -211,6 +212,58 @@ class PengRobinson : public MixtureFugacityTP double d2aAlpha_dT2() const; public: + void setState_DP(double rho, double p, double Tguess=300) override { + if (p <= 0) { + throw CanteraError("IdealGasPhase::setState_DP", + "pressure must be positive"); + } + //"Yeah its this one" + setDensity(rho); + setForcedSolutionBranch(-2); + double mV = meanMolecularWeight() / rho; + + // Update individual alpha + double aCalc, bCalc, aAlpha; + + + int iteration = 0; + double T_new = Tguess; //(p + term2) * denum / GasConstant; + double p_iter=p; + + while (iteration < 300) { + // Update a, b, and alpha based on the current guess for T + updateMixingExpressions(T_new, aCalc, bCalc, aAlpha); + // Calculate p using PR EoS + double denom = mV * mV + 2 * mV * bCalc - bCalc * bCalc; + double mV_b = mV - bCalc; + double term2 = (aAlpha) / denom; + p_iter = (GasConstant*T_new)/mV_b - term2; + //std::cout << "T_ : "< 50) && (rho<1000)) + // { + // std::cout <<"T_new: "< 1e-10) { + std::cout<<"yn: "< 1.0E-14) && (fabs(res) > 1.0E-14 * fabs(dresdV) * fabs(Vroot[i]))) { writelog("MixtureFugacityTP::solveCubic(T = {}, p = {}): " - "WARNING root didn't converge V = {}", T, pres, Vroot[i]); + "WARNING root didn't converge V = {}\n" + "X1: {}, X2: {}\n", + T, pres, Vroot[i],moleFractions_[0],moleFractions_[1]); writelogendl(); } } diff --git a/src/thermo/PengRobinson.cpp b/src/thermo/PengRobinson.cpp index 4bb7dd2139f..8fc4d95a054 100644 --- a/src/thermo/PengRobinson.cpp +++ b/src/thermo/PengRobinson.cpp @@ -673,6 +673,28 @@ void PengRobinson::updateMixingExpressions() calculateAB(m_a,m_b,m_aAlpha_mix); } + +void PengRobinson::updateMixingExpressions(double& temp,double& aCalc, double& bCalc, double& aAlphaCalc) +{ + //double temp = temperature(); + + // Update individual alpha + for (size_t j = 0; j < m_kk; j++) { + double critTemp_j = speciesCritTemperature(m_a_coeffs(j,j), m_b_coeffs[j]); + double sqt_alpha = 1 + m_kappa[j] * (1 - sqrt(temp / critTemp_j)); + m_alpha[j] = sqt_alpha*sqt_alpha; + } + + // Update aAlpha_i, j + for (size_t i = 0; i < m_kk; i++) { + for (size_t j = 0; j < m_kk; j++) { + m_aAlpha_binary(i, j) = sqrt(m_alpha[i] * m_alpha[j]) * m_a_coeffs(i,j); + } + } + calculateAB(aCalc,bCalc,aAlphaCalc); +} + + void PengRobinson::calculateAB(double& aCalc, double& bCalc, double& aAlphaCalc) const { bCalc = 0.0; @@ -710,6 +732,30 @@ double PengRobinson::daAlpha_dT() const return daAlphadT; } +double PengRobinson::daAlpha_dT(double& temp) const +{ + double daAlphadT = 0.0, k, Tc, sqtTr, coeff1, coeff2; + for (size_t i = 0; i < m_kk; i++) { + // Calculate first derivative of alpha for individual species + Tc = speciesCritTemperature(m_a_coeffs(i,i), m_b_coeffs[i]); + sqtTr = sqrt(temp / Tc); + coeff1 = 1 / (Tc*sqtTr); + coeff2 = sqtTr - 1; + k = m_kappa[i]; + m_dalphadT[i] = coeff1 * (k*k*coeff2 - k); + } + // Calculate mixture derivative + for (size_t i = 0; i < m_kk; i++) { + for (size_t j = 0; j < m_kk; j++) { + daAlphadT += moleFractions_[i] * moleFractions_[j] * 0.5 + * m_aAlpha_binary(i, j) + * (m_dalphadT[i] / m_alpha[i] + m_dalphadT[j] / m_alpha[j]); + } + } + return daAlphadT; +} + + double PengRobinson::d2aAlpha_dT2() const { for (size_t i = 0; i < m_kk; i++) { From fa6d4d9734954cfb28138b6e089255a3340f4912 Mon Sep 17 00:00:00 2001 From: jeremy karpowski Date: Mon, 3 Feb 2025 10:56:23 +0100 Subject: [PATCH 2/6] added ForceState option to thermo-phase --- include/cantera/thermo/MixtureFugacityTP.h | 2 +- include/cantera/thermo/PengRobinson.h | 53 +--------------------- include/cantera/thermo/ThermoPhase.h | 8 ++++ src/thermo/PengRobinson.cpp | 53 ++++++++++++++++++++++ 4 files changed, 63 insertions(+), 53 deletions(-) diff --git a/include/cantera/thermo/MixtureFugacityTP.h b/include/cantera/thermo/MixtureFugacityTP.h index 76431d7d959..0965533354d 100644 --- a/include/cantera/thermo/MixtureFugacityTP.h +++ b/include/cantera/thermo/MixtureFugacityTP.h @@ -88,7 +88,7 @@ class MixtureFugacityTP : public ThermoPhase * -1 means gas. The value -2 means unrestricted. Values of zero or * greater refer to species dominated condensed phases. */ - void setForcedSolutionBranch(int solnBranch); + void setForcedSolutionBranch(int solnBranch) override; //! Report the solution branch which the solution is restricted to /*! diff --git a/include/cantera/thermo/PengRobinson.h b/include/cantera/thermo/PengRobinson.h index 03b877136ad..7c7d0820918 100644 --- a/include/cantera/thermo/PengRobinson.h +++ b/include/cantera/thermo/PengRobinson.h @@ -212,58 +212,7 @@ class PengRobinson : public MixtureFugacityTP double d2aAlpha_dT2() const; public: - void setState_DP(double rho, double p, double Tguess=300) override { - if (p <= 0) { - throw CanteraError("IdealGasPhase::setState_DP", - "pressure must be positive"); - } - //"Yeah its this one" - setDensity(rho); - setForcedSolutionBranch(-2); - double mV = meanMolecularWeight() / rho; - - // Update individual alpha - double aCalc, bCalc, aAlpha; - - - int iteration = 0; - double T_new = Tguess; //(p + term2) * denum / GasConstant; - double p_iter=p; - - while (iteration < 300) { - // Update a, b, and alpha based on the current guess for T - updateMixingExpressions(T_new, aCalc, bCalc, aAlpha); - // Calculate p using PR EoS - double denom = mV * mV + 2 * mV * bCalc - bCalc * bCalc; - double mV_b = mV - bCalc; - double term2 = (aAlpha) / denom; - p_iter = (GasConstant*T_new)/mV_b - term2; - //std::cout << "T_ : "< 50) && (rho<1000)) - // { - // std::cout <<"T_new: "< 50) && (rho<1000)) + // { + // std::cout <<"T_new: "< Date: Mon, 3 Feb 2025 12:42:24 +0100 Subject: [PATCH 3/6] removed debuging message --- src/thermo/MixtureFugacityTP.cpp | 5 ----- 1 file changed, 5 deletions(-) diff --git a/src/thermo/MixtureFugacityTP.cpp b/src/thermo/MixtureFugacityTP.cpp index bd3e5e0635f..0319933de8c 100644 --- a/src/thermo/MixtureFugacityTP.cpp +++ b/src/thermo/MixtureFugacityTP.cpp @@ -834,11 +834,6 @@ int MixtureFugacityTP::solveCubic(double T, double pres, double a, double b, //check if y = h if (fabs(fabs(h) - fabs(yN)) < 1.0E-10) { if (disc > 1e-10) { - std::cout<<"yn: "< Date: Wed, 5 Feb 2025 14:58:50 +0100 Subject: [PATCH 4/6] setState_DP is now also tested in the PengRobinson test suite (CoolPropValidateDP) --- test/thermo/PengRobinson_Test.cpp | 30 ++++++++++++++++++++++++++++++ 1 file changed, 30 insertions(+) diff --git a/test/thermo/PengRobinson_Test.cpp b/test/thermo/PengRobinson_Test.cpp index e1a19251769..6c29e81bd13 100644 --- a/test/thermo/PengRobinson_Test.cpp +++ b/test/thermo/PengRobinson_Test.cpp @@ -179,6 +179,36 @@ TEST_F(PengRobinson_Test, CoolPropValidate) } } +TEST_F(PengRobinson_Test, CoolPropValidateDP) +{ + // Validate the P-R EoS in Cantera with P-R EoS from CoolProp + + const double rhoCoolProp[10] = { + 9.067928191884574, + 8.318322900591179, + 7.6883521740498155, + 7.150504298001246, + 6.685330199667018, + 6.278630757480957, + 5.919763108091383, + 5.600572727499541, + 5.314694056926007, + 5.057077678380463 + }; + + double p = 5e5; + + // Calculate Temperature using Peng-Robinson EoS from Cantera + for(int i=0; i<10; i++) + { + const double tempReference = 300 + i*25; + set_r(1.0); + test_phase->setState_DP(rhoCoolProp[i], p,300); + + EXPECT_NEAR(test_phase->temperature(),tempReference,1.e-3); + } +} + TEST(PengRobinson, lookupSpeciesProperties) { AnyMap phase_def = AnyMap::fromYamlString( From 5249dcef215fa88615473b24e42487d0a6e11c88 Mon Sep 17 00:00:00 2001 From: "T. Jeremy P. Karpowski" Date: Fri, 14 Feb 2025 14:30:31 +0100 Subject: [PATCH 5/6] removed trailing whitespace - Lint checks complained about trailing whitespaces. These should be solve now. --- src/thermo/PengRobinson.cpp | 18 +++--------------- 1 file changed, 3 insertions(+), 15 deletions(-) diff --git a/src/thermo/PengRobinson.cpp b/src/thermo/PengRobinson.cpp index ad4b5a5d7ec..101e022f3d1 100644 --- a/src/thermo/PengRobinson.cpp +++ b/src/thermo/PengRobinson.cpp @@ -838,53 +838,41 @@ AnyMap PengRobinson::getAuxiliaryData() return parameters; } -void PengRobinson::setState_DP(double rho, double p, double Tguess) +void PengRobinson::setState_DP(double rho, double p, double Tguess) { if (p <= 0 || rho<=0 ) { throw CanteraError("PengRobinson::setState_DP", "pressure and density must be positive"); } - //"Yeah its this one" setDensity(rho); double mV = meanMolecularWeight() / rho; - // Update individual alpha double aCalc, bCalc, aAlpha; - - int iteration = 0; double T_new = Tguess; //(p + term2) * denum / GasConstant; double p_iter=p; - while (iteration < 300) { - // Update a, b, and alpha based on the current guess for T + // Update a, b, and alpha based on the current guess for T updateMixingExpressions(T_new, aCalc, bCalc, aAlpha); // Calculate p using PR EoS double denom = mV * mV + 2 * mV * bCalc - bCalc * bCalc; double mV_b = mV - bCalc; double term2 = (aAlpha) / denom; p_iter = (GasConstant*T_new)/mV_b - term2; - //std::cout << "T_ : "< 50) && (rho<1000)) - // { - // std::cout <<"T_new: "< Date: Fri, 28 Feb 2025 11:23:29 +0100 Subject: [PATCH 6/6] API correction: Added tolerance as third parameter to be consistent with other setState functions. Enabled tolerance and TGuess in the python API .DP .DPX .DPY --- include/cantera/thermo/IdealGasPhase.h | 2 +- include/cantera/thermo/PengRobinson.h | 13 +++++++++++- include/cantera/thermo/ThermoPhase.h | 4 +++- interfaces/cython/cantera/thermo.pyx | 28 ++++++++++++++++++++------ src/thermo/PengRobinson.cpp | 4 ++-- test/thermo/PengRobinson_Test.cpp | 2 +- 6 files changed, 41 insertions(+), 12 deletions(-) diff --git a/include/cantera/thermo/IdealGasPhase.h b/include/cantera/thermo/IdealGasPhase.h index d9d53a3e2cb..7c0fbaee5cd 100644 --- a/include/cantera/thermo/IdealGasPhase.h +++ b/include/cantera/thermo/IdealGasPhase.h @@ -366,7 +366,7 @@ class IdealGasPhase: public ThermoPhase * @param p Pressure (Pa) * @since New in %Cantera 3.0. */ - void setState_DP(double rho, double p, double Tguess=300) override { + void setState_DP(double rho, double p,double tol=1e-9, double Tguess=300) override { if (p <= 0) { throw CanteraError("IdealGasPhase::setState_DP", "pressure must be positive"); diff --git a/include/cantera/thermo/PengRobinson.h b/include/cantera/thermo/PengRobinson.h index 7c7d0820918..cf8c58c8c07 100644 --- a/include/cantera/thermo/PengRobinson.h +++ b/include/cantera/thermo/PengRobinson.h @@ -212,7 +212,18 @@ class PengRobinson : public MixtureFugacityTP double d2aAlpha_dT2() const; public: - void setState_DP(double rho, double p, double Tguess=300) override; + //! Set the density (kg/m**3) and pressure (Pa) at constant composition + /*! + * Within %Cantera, the independent variable of the PR-EoS is the density and Temperature. Therefore, this function solves for + * the temperature that will yield the desired input pressure at the provided density. + * The composition is held constant during this process. + * + * @param rho Density (kg/m^3) + * @param p Pressure (Pa) + * @param tol tolerance for the iterative solver + * @param TGuess guess for the temperature (K) used as start for the iterative solver + */ + void setState_DP(double rho, double p, double tol=1e-9, double Tguess=300) override; double isothermalCompressibility() const override; double thermalExpansionCoeff() const override; diff --git a/include/cantera/thermo/ThermoPhase.h b/include/cantera/thermo/ThermoPhase.h index fa1f84a92b5..b2eed07ecfa 100644 --- a/include/cantera/thermo/ThermoPhase.h +++ b/include/cantera/thermo/ThermoPhase.h @@ -1399,9 +1399,11 @@ class ThermoPhase : public Phase * * @param rho Density (kg/m^3) * @param p Pressure (Pa) + * @param tol tolerance for the iterative solver for cubic EoS + * @param TGuess guess for the temperature (K) used as start for the iterative solver for cubic EoS * @since New in %Cantera 3.0. */ - virtual void setState_DP(double rho, double p, double Tguess=300) { + virtual void setState_DP(double rho, double p, double tol=1e-9, double Tguess=300) { throw NotImplementedError("ThermoPhase::setState_DP", "Not implemented for phase type '{}'", type()); } diff --git a/interfaces/cython/cantera/thermo.pyx b/interfaces/cython/cantera/thermo.pyx index b091f708829..72f21d278ac 100644 --- a/interfaces/cython/cantera/thermo.pyx +++ b/interfaces/cython/cantera/thermo.pyx @@ -1486,32 +1486,48 @@ cdef class ThermoPhase(_SolutionBase): def __get__(self): return self.density, self.P def __set__(self, values): - assert len(values) == 2, 'incorrect number of values' + assert ((len(values) == 2) || (len(values) == 4)), 'incorrect number of values either 2 (DP) or 4 (D, P, tolerance, TGuess)' D = values[0] if values[0] is not None else self.density P = values[1] if values[1] is not None else self.P - self.thermo.setState_DP(D*self._mass_factor(), P) + if (len(values) == 2): + self.thermo.setState_DP(D*self._mass_factor(), P) + else: + tolerance = values[2] + TGuess = values[3] + self.thermo.setState_DP(D*self._mass_factor(), P,tolerance,TGuess) property DPX: """Get/Set density [kg/m^3], pressure [Pa], and mole fractions.""" def __get__(self): return self.density, self.P, self.X def __set__(self, values): - assert len(values) == 3, 'incorrect number of values' + assert ((len(values) == 3) || (len(values) == 5)), 'incorrect number of values either 3 (DPX) or 5 (D, P, X, tolerance, TGuess)' + if D = values[0] if values[0] is not None else self.density P = values[1] if values[1] is not None else self.P self.X = values[2] - self.thermo.setState_DP(D*self._mass_factor(), P) + if (len(values) == 3): + self.thermo.setState_DP(D*self._mass_factor(), P) + else: + tolerance = values[3] + TGuess = values[4] + self.thermo.setState_DP(D*self._mass_factor(), P,tolerance,TGuess) property DPY: """Get/Set density [kg/m^3], pressure [Pa], and mass fractions.""" def __get__(self): return self.density, self.P, self.Y def __set__(self, values): - assert len(values) == 3, 'incorrect number of values' + assert ((len(values) == 3) || (len(values) == 5)), 'incorrect number of values either 3 (DPY) or 5 (D, P, Y, tolerance, TGuess)' D = values[0] if values[0] is not None else self.density P = values[1] if values[1] is not None else self.P self.Y = values[2] - self.thermo.setState_DP(D*self._mass_factor(), P) + if (len(values) == 3): + self.thermo.setState_DP(D*self._mass_factor(), P) + else: + tolerance = values[3] + TGuess = values[4] + self.thermo.setState_DP(D*self._mass_factor(), P,tolerance,TGuess) property HP: """Get/Set enthalpy [J/kg or J/kmol] and pressure [Pa].""" diff --git a/src/thermo/PengRobinson.cpp b/src/thermo/PengRobinson.cpp index 101e022f3d1..03c87cc95a3 100644 --- a/src/thermo/PengRobinson.cpp +++ b/src/thermo/PengRobinson.cpp @@ -838,7 +838,7 @@ AnyMap PengRobinson::getAuxiliaryData() return parameters; } -void PengRobinson::setState_DP(double rho, double p, double Tguess) +void PengRobinson::setState_DP(double rho, double p,double tol, double Tguess) { if (p <= 0 || rho<=0 ) { throw CanteraError("PengRobinson::setState_DP", @@ -861,7 +861,7 @@ void PengRobinson::setState_DP(double rho, double p, double Tguess) p_iter = (GasConstant*T_new)/mV_b - term2; // Check if the result is within the specified tolerance double error =(p_iter-p)/p; - if (fabs(error)<1e-10) { + if (fabs(error)setState_DP(rhoCoolProp[i], p,300); + test_phase->setState_DP(rhoCoolProp[i], p,1e-10,300); EXPECT_NEAR(test_phase->temperature(),tempReference,1.e-3); }