diff --git a/include/cantera/thermo/IdealGasPhase.h b/include/cantera/thermo/IdealGasPhase.h index 22178cc706f..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) 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/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 8fda6dd19d1..cf8c58c8c07 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,18 @@ class PengRobinson : public MixtureFugacityTP double d2aAlpha_dT2() const; public: + //! 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; @@ -229,7 +242,7 @@ class PengRobinson : public MixtureFugacityTP * the internal numbers based on the state of the object. */ void updateMixingExpressions() override; - + void updateMixingExpressions(double& temp,double& aCalc, double& bCalc, double& aAlphaCalc); //! Calculate the @f$ a @f$, @f$ b @f$, and @f$ \alpha @f$ parameters given the temperature /*! * This function doesn't change the internal state of the object, so it is a diff --git a/include/cantera/thermo/ThermoPhase.h b/include/cantera/thermo/ThermoPhase.h index ec3952983dc..b2eed07ecfa 100644 --- a/include/cantera/thermo/ThermoPhase.h +++ b/include/cantera/thermo/ThermoPhase.h @@ -1399,12 +1399,22 @@ 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) { + 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()); } + //! Set the solution branch to force the ThermoPhase to exist on one branch + //! or another + /*! + * @param solnBranch Branch that the solution is restricted to. the value + * -1 means gas. The value -2 means unrestricted. Values of zero or + * greater refer to species dominated condensed phases. + */ + virtual void setForcedSolutionBranch(int solnBranch){}; //! Set the state using an AnyMap containing any combination of properties //! supported by the thermodynamic model 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/MixtureFugacityTP.cpp b/src/thermo/MixtureFugacityTP.cpp index 39c1697d95d..0319933de8c 100644 --- a/src/thermo/MixtureFugacityTP.cpp +++ b/src/thermo/MixtureFugacityTP.cpp @@ -835,7 +835,8 @@ int MixtureFugacityTP::solveCubic(double T, double pres, double a, double b, if (fabs(fabs(h) - fabs(yN)) < 1.0E-10) { if (disc > 1e-10) { throw CanteraError("MixtureFugacityTP::solveCubic", - "value of yN and h are too high, unrealistic roots may be obtained"); + "value of yN and h are too high, unrealistic roots may be obtained" + ); } disc = 0.0; } diff --git a/src/thermo/PengRobinson.cpp b/src/thermo/PengRobinson.cpp index 82d0db4be66..03c87cc95a3 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++) { @@ -792,4 +838,45 @@ AnyMap PengRobinson::getAuxiliaryData() return parameters; } +void PengRobinson::setState_DP(double rho, double p,double tol, double Tguess) +{ + if (p <= 0 || rho<=0 ) { + throw CanteraError("PengRobinson::setState_DP", + "pressure and density must be positive"); + } + 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 + 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; + // Check if the result is within the specified tolerance + double error =(p_iter-p)/p; + if (fabs(error)setState_DP(rhoCoolProp[i], p,1e-10,300); + + EXPECT_NEAR(test_phase->temperature(),tempReference,1.e-3); + } +} + TEST(PengRobinson, lookupSpeciesProperties) { AnyMap phase_def = AnyMap::fromYamlString(