diff --git a/include/cantera/thermo/SoaveRedlichKwong.h b/include/cantera/thermo/SoaveRedlichKwong.h new file mode 100644 index 0000000000..e1d56a8524 --- /dev/null +++ b/include/cantera/thermo/SoaveRedlichKwong.h @@ -0,0 +1,329 @@ +//! @file SoaveRedlichKwong.h + +// This file is part of Cantera. See License.txt in the top-level directory or +// at https://cantera.org/license.txt for license and copyright information. + +#ifndef CT_SOAVEREDLICHKWONG_H +#define CT_SOAVEREDLICHKWONG_H + +#include "MixtureFugacityTP.h" +#include "cantera/base/Array.h" + +namespace Cantera +{ +/** + * Implementation of a multi-species Soave-Redlich-Kwong equation of state + * + * @ingroup thermoprops + */ +class SoaveRedlichKwong : public MixtureFugacityTP +{ +public: + + //! Construct and initialize a SoaveRedlichKwong object directly from an + //! input file + /*! + * @param infile Name of the input file containing the phase YAML data. + * If blank, an empty phase will be created. + * @param id ID of the phase in the input file. If empty, the + * first phase definition in the input file will be used. + */ + explicit SoaveRedlichKwong(const std::string& infile="", + const std::string& id=""); + + virtual std::string type() const { + return "Soave-Redlich-Kwong"; + } + + //! @name Molar Thermodynamic properties + //! @{ + + virtual double cp_mole() const; + virtual double cv_mole() const; + + //! @} + //! @name Mechanical Properties + //! @{ + + //! Return the thermodynamic pressure (Pa). + /*! + * Since the mass density, temperature, and mass fractions are stored, + * this method uses these values to implement the + * mechanical equation of state \f$ P(T, \rho, Y_1, \dots, Y_K) \f$. + * + * \f[ + * P = \frac{RT}{v-b_{mix}} + * - \frac{\left(a\alpha\right)_{mix}}{v\left(v + b_{mix}\right)} + * \f] + * + * where: + * + * \f[ + * \alpha = \left[ 1 + \kappa \left(1-T_r^{0.5}\right)\right]^2 + * \f] + * + * and + * + * \f[ + * \kappa = \left(0.48508 + 1.55171\omega - 0.15613\omega^2\right) + * \f] + * + * Coefficients \f$ a_{mix}, b_{mix} \f$ and \f$(a \alpha)_{mix}\f$ are calculated as + * + * \f[ + * a_{mix} = \sum_i \sum_j X_i X_j a_{i, j} = \sum_i \sum_j X_i X_j \sqrt{a_i a_j} + * \f] + * + * \f[ + * b_{mix} = \sum_i X_i b_i + * \f] + * + * \f[ + * {a \alpha}_{mix} = \sum_i \sum_j X_i X_j {a \alpha}_{i, j} + * = \sum_i \sum_j X_i X_j \sqrt{a_i a_j} \sqrt{\alpha_i \alpha_j} + * \f] + */ + virtual double pressure() const; + + //! @} + + //! Returns the standard concentration \f$ C^0_k \f$, which is used to + //! normalize the generalized concentration. + /*! + * This is defined as the concentration by which the generalized + * concentration is normalized to produce the activity. + * The ideal gas mixture is considered as the standard or reference state here. + * Since the activity for an ideal gas mixture is simply the mole fraction, + * for an ideal gas, \f$ C^0_k = P/\hat R T \f$. + * + * @param k Optional parameter indicating the species. The default is to + * assume this refers to species 0. + * @return + * Returns the standard Concentration in units of m^3 / kmol. + */ + virtual double standardConcentration(size_t k=0) const; + + //! Get the array of non-dimensional activity coefficients at the current + //! solution temperature, pressure, and solution concentration. + /*! + * For all objects with the Mixture Fugacity approximation, we define the + * standard state as an ideal gas at the current temperature and pressure of + * the solution. The activities are based on this standard state. + * + * @param ac Output vector of activity coefficients. Length: m_kk. + */ + virtual void getActivityCoefficients(double* ac) const; + + //! @name Partial Molar Properties of the Solution + //! @{ + + virtual void getChemPotentials(double* mu) const; + virtual void getPartialMolarEnthalpies(double* hbar) const; + virtual void getPartialMolarEntropies(double* sbar) const; + virtual void getPartialMolarIntEnergies(double* ubar) const; + //! Calculate species-specific molar specific heats + /*! + * This function is currently not implemented for Soave-Redlich-Kwong phase. + */ + virtual void getPartialMolarCp(double* cpbar) const; + virtual void getPartialMolarVolumes(double* vbar) const; + //! @} + + //! Calculate species-specific critical temperature + /*! + * The temperature dependent parameter in SRK EoS is calculated as + * \f[ T_{crit} = (0.0778 a)/(0.4572 b R) \f] TODO + * Units: Kelvin + * + * @param a species-specific coefficients used in SRK EoS + * @param b species-specific coefficients used in SRK EoS + */ + virtual double speciesCritTemperature(double a, double b) const; + + //! @name Initialization Methods - For Internal use + //! + //! The following methods are used in the process of constructing + //! the phase and setting its parameters from a specification in an + //! input file. They are not normally used in application programs. + //! To see how they are used, see importPhase(). + //! @{ + + virtual bool addSpecies(shared_ptr spec); + virtual void initThermo(); + virtual void getSpeciesParameters(const std::string& name, + AnyMap& speciesNode) const; + + //! Set the pure fluid interaction parameters for a species + /*! + * @param species Name of the species + * @param a \f$a\f$ parameter in the Soave-Redlich-Kwong model [Pa-m^6/kmol^2] + * @param b \f$a\f$ parameter in the Soave-Redlich-Kwong model [m^3/kmol] + * @param w acentric factor + */ + void setSpeciesCoeffs(const std::string& species, double a, double b, + double w); + + //! Set values for the interaction parameter between two species + /*! + * @param species_i Name of one species + * @param species_j Name of the other species + * @param a \f$a\f$ parameter in the Soave-Redlich-Kwong model [Pa-m^6/kmol^2] + */ + void setBinaryCoeffs(const std::string& species_i, + const std::string& species_j, double a); + //! @} + +protected: + // Special functions inherited from MixtureFugacityTP + virtual double sresid() const; + virtual double hresid() const; + + virtual double liquidVolEst(double T, double& pres) const; + virtual double densityCalc(double T, double pressure, int phase, double rhoguess); + + virtual double densSpinodalLiquid() const; + virtual double densSpinodalGas() const; + virtual double dpdVCalc(double T, double molarVol, double& presCalc) const; + + // Special functions not inherited from MixtureFugacityTP + + //! Calculate temperature derivative \f$d(a \alpha)/dT\f$ + /*! + * These are stored internally. + */ + double daAlpha_dT() const; + + //! Calculate second derivative \f$d^2(a \alpha)/dT^2\f$ + /*! + * These are stored internally. + */ + double d2aAlpha_dT2() const; + +public: + //! Returns the isothermal compressibility. Units: 1/Pa. + double isothermalCompressibility() const; + + //! Return the volumetric thermal expansion coefficient. Units: 1/K. + double thermalExpansionCoeff() const; + + //! Calculate \f$dp/dV\f$ and \f$dp/dT\f$ at the current conditions + /*! + * These are stored internally. + */ + void calculatePressureDerivatives() const; + + //! Update the \f$a\f$, \f$b\f$, and \f$\alpha\f$ parameters + /*! + * The \f$a\f$ and the \f$b\f$ parameters depend on the mole fraction and the + * parameter \f$\alpha\f$ depends on the temperature. This function updates + * the internal numbers based on the state of the object. + */ + virtual void updateMixingExpressions(); + + //! Calculate the \f$a\f$, \f$b\f$, and \f$a \alpha\f$ parameters given the temperature + /*! + * This function doesn't change the internal state of the object, so it is a + * const function. It does use the stored mole fractions in the object. + * + * @param aCalc (output) Returns the a value + * @param bCalc (output) Returns the b value. + * @param aAlpha (output) Returns the (a*alpha) value. + */ + void calculateAB(double& aCalc, double& bCalc, double& aAlphaCalc) const; + + void calcCriticalConditions(double& pc, double& tc, double& vc) const; + + //! Prepare variables and call the function to solve the cubic equation of state + int solveCubic(double T, double pres, double a, double b, double aAlpha, + double Vroot[3]) const; + +protected: + //! Value of \f$b\f$ in the equation of state + /*! + * `m_b` is a function of the mole fractions and species-specific b values. + */ + double m_b; + + //! Value of \f$a\f$ in the equation of state + /*! + * `m_a` depends only on the mole fractions. + */ + double m_a; + + //! Value of \f$\alpha\f$ in the equation of state + /*! + * `m_aAlpha_mix` is a function of the temperature and the mole fractions. + */ + double m_aAlpha_mix; + + // Vectors required to store species-specific a_coeff, b_coeff, alpha, kappa + // and other derivatives. Length = m_kk. + vector_fp m_b_coeffs; + vector_fp m_kappa; + vector_fp m_acentric; //!< acentric factor for each species, length #m_kk + mutable vector_fp m_dalphadT; + mutable vector_fp m_d2alphadT2; + vector_fp m_alpha; + + // Matrices for Binary coefficients a_{i,j} and {a*alpha}_{i.j} are saved in an + // array form. Size = (m_kk, m_kk). + Array2D m_a_coeffs; + Array2D m_aAlpha_binary; + + //! Explicitly-specified binary interaction parameters, to enable serialization + std::map> m_binaryParameters; + + int m_NSolns; + + double m_Vroot[3]; + + //! Temporary storage - length = m_kk. + mutable vector_fp m_pp; + + // Partial molar volumes of the species + mutable vector_fp m_partialMolarVolumes; + + //! The derivative of the pressure with respect to the volume + /*! + * Calculated at the current conditions. temperature and mole number kept + * constant + */ + mutable double m_dpdV; + + //! The derivative of the pressure with respect to the temperature + /*! + * Calculated at the current conditions. Total volume and mole number kept + * constant + */ + mutable double m_dpdT; + + //! Vector of derivatives of pressure with respect to mole number + /*! + * Calculated at the current conditions. Total volume, temperature and + * other mole number kept constant + */ + mutable vector_fp m_dpdni; + + enum class CoeffSource { EoS, CritProps, Database }; + //! For each species, specifies the source of the a, b, and omega coefficients + std::vector m_coeffSource; + +private: + //! Omega constant: omega_a used in Soave-Redlich-Kwong equation of state + /*! + * This value is calculated by solving SRK cubic equation at the critical point. + */ + static const double omega_a; + + //! Omega constant: omega_b used in Soave-Redlich-Kwong equation of state + /*! + * This value is calculated by solving SRK cubic equation at the critical point. + */ + static const double omega_b; + + //! Omega constant for the critical molar volume + static const double omega_vc; +}; +} + +#endif diff --git a/src/thermo/SoaveRedlichKwong.cpp b/src/thermo/SoaveRedlichKwong.cpp new file mode 100644 index 0000000000..62b2c43e6b --- /dev/null +++ b/src/thermo/SoaveRedlichKwong.cpp @@ -0,0 +1,633 @@ +//! @file SoaveRedlichKwong.cpp + +// This file is part of Cantera. See License.txt in the top-level directory or +// at https://cantera.org/license.txt for license and copyright information. + +#include "cantera/thermo/SoaveRedlichKwong.h" +#include "cantera/thermo/ThermoFactory.h" +#include "cantera/thermo/Species.h" +#include "cantera/base/stringUtils.h" + +#include + +namespace bmt = boost::math::tools; + +namespace Cantera +{ + +const double SoaveRedlichKwong::omega_a = 4.27480233540E-01; +const double SoaveRedlichKwong::omega_b = 8.66403499650E-02; +const double SoaveRedlichKwong::omega_vc = 3.33333333333333E-01; + +SoaveRedlichKwong::SoaveRedlichKwong(const std::string& infile, const std::string& id_) : + m_b(0.0), + m_a(0.0), + m_aAlpha_mix(0.0), + m_NSolns(0), + m_dpdV(0.0), + m_dpdT(0.0) +{ + std::fill_n(m_Vroot, 3, 0.0); + initThermoFile(infile, id_); +} + +void SoaveRedlichKwong::setSpeciesCoeffs(const std::string& species, double a, double b, + double w) +{ + size_t k = speciesIndex(species); + if (k == npos) { + throw CanteraError("SoaveRedlichKwong::setSpeciesCoeffs", + "Unknown species '{}'.", species); + } + + // Calculate value of kappa (independent of temperature) + // w is an acentric factor of species + m_kappa[k] = 0.480 + 1.574*w - 0.176*w*w; + m_acentric[k] = w; // store the original acentric factor to enable serialization + + // Calculate alpha (temperature dependent interaction parameter) + double critTemp = speciesCritTemperature(a, b); + double sqt_T_r = sqrt(temperature() / critTemp); + double sqt_alpha = 1 + m_kappa[k] * (1 - sqt_T_r); + m_alpha[k] = sqt_alpha*sqt_alpha; + + m_a_coeffs(k,k) = a; + double aAlpha_k = a*m_alpha[k]; + m_aAlpha_binary(k,k) = aAlpha_k; + + // standard mixing rule for cross-species interaction term + for (size_t j = 0; j < m_kk; j++) { + if (k == j) { + continue; + } + double a0kj = sqrt(m_a_coeffs(j,j) * a); + double aAlpha_j = a*m_alpha[j]; + double a_Alpha = sqrt(aAlpha_j*aAlpha_k); + if (m_a_coeffs(j, k) == 0) { + m_a_coeffs(j, k) = a0kj; + m_aAlpha_binary(j, k) = a_Alpha; + m_a_coeffs(k, j) = a0kj; + m_aAlpha_binary(k, j) = a_Alpha; + } + } + m_b_coeffs[k] = b; +} + +void SoaveRedlichKwong::setBinaryCoeffs(const std::string& species_i, + const std::string& species_j, double a0) +{ + size_t ki = speciesIndex(species_i); + if (ki == npos) { + throw CanteraError("SoaveRedlichKwong::setBinaryCoeffs", + "Unknown species '{}'.", species_i); + } + size_t kj = speciesIndex(species_j); + if (kj == npos) { + throw CanteraError("SoaveRedlichKwong::setBinaryCoeffs", + "Unknown species '{}'.", species_j); + } + + m_a_coeffs(ki, kj) = m_a_coeffs(kj, ki) = a0; + m_binaryParameters[species_i][species_j] = a0; + m_binaryParameters[species_j][species_i] = a0; + // Calculate alpha_ij + double alpha_ij = m_alpha[ki] * m_alpha[kj]; + m_aAlpha_binary(ki, kj) = m_aAlpha_binary(kj, ki) = a0*alpha_ij; +} + +// ------------Molar Thermodynamic Properties ------------------------- + +double SoaveRedlichKwong::cp_mole() const +{ + return 1.0; +} + +double SoaveRedlichKwong::cv_mole() const +{ + _updateReferenceStateThermo(); + double T = temperature(); + calculatePressureDerivatives(); + return (cp_mole() + T * m_dpdT * m_dpdT / m_dpdV); +} + +double SoaveRedlichKwong::pressure() const +{ + _updateReferenceStateThermo(); + // Get a copy of the private variables stored in the State object + double T = temperature(); + double mv = molarVolume(); + return GasConstant * T / (mv - m_b) - m_aAlpha_mix / (mv * (mv + m_b)); +} + +double SoaveRedlichKwong::standardConcentration(size_t k) const +{ + getStandardVolumes(m_tmpV.data()); + return 1.0 / m_tmpV[k]; +} + +void SoaveRedlichKwong::getActivityCoefficients(double* ac) const +{ + +} + +// ---- Partial Molar Properties of the Solution ----------------- + +void SoaveRedlichKwong::getChemPotentials(double* mu) const +{ + +} + +void SoaveRedlichKwong::getPartialMolarEnthalpies(double* hbar) const +{ + +} + +void SoaveRedlichKwong::getPartialMolarEntropies(double* sbar) const +{ + +} + +void SoaveRedlichKwong::getPartialMolarIntEnergies(double* ubar) const +{ + +} + +void SoaveRedlichKwong::getPartialMolarCp(double* cpbar) const +{ + +} + +void SoaveRedlichKwong::getPartialMolarVolumes(double* vbar) const +{ + +} + +double SoaveRedlichKwong::speciesCritTemperature(double a, double b) const +{ + if (b <= 0.0) { + return 1000000.; + } else if (a <= 0.0) { + return 0.0; + } else { + return a * omega_b / (b * omega_a * GasConstant); + } +} + +bool SoaveRedlichKwong::addSpecies(shared_ptr spec) +{ + bool added = MixtureFugacityTP::addSpecies(spec); + if (added) { + m_a_coeffs.resize(m_kk, m_kk, 0.0); + m_b_coeffs.push_back(0.0); + m_aAlpha_binary.resize(m_kk, m_kk, 0.0); + m_kappa.push_back(0.0); + m_acentric.push_back(0.0); + m_alpha.push_back(0.0); + m_dalphadT.push_back(0.0); + m_d2alphadT2.push_back(0.0); + m_pp.push_back(0.0); + m_partialMolarVolumes.push_back(0.0); + m_dpdni.push_back(0.0); + m_coeffSource.push_back(CoeffSource::EoS); + } + return added; +} + +void SoaveRedlichKwong::initThermo() +{ + // Contents of 'critical-properties.yaml', loaded later if needed + AnyMap critPropsDb; + std::unordered_map dbSpecies; + + for (auto& item : m_species) { + auto& data = item.second->input; + size_t k = speciesIndex(item.first); + if (m_a_coeffs(k, k) != 0.0) { + continue; + } + bool foundCoeffs = false; + if (data.hasKey("equation-of-state") && + data["equation-of-state"].hasMapWhere("model", "Soave-Redlich-Kwong")) + { + // Read a and b coefficients and acentric factor w_ac from species input + // information, specified in a YAML input file. + auto eos = data["equation-of-state"].getMapWhere( + "model", "Soave-Redlich-Kwong"); + if (eos.hasKey("a") && eos.hasKey("b") && eos.hasKey("acentric-factor")) { + double a0 = eos.convert("a", "Pa*m^6/kmol^2"); + double b = eos.convert("b", "m^3/kmol"); + // unitless acentric factor: + double w = eos["acentric-factor"].asDouble(); + setSpeciesCoeffs(item.first, a0, b, w); + foundCoeffs = true; + } + + if (eos.hasKey("binary-a")) { + AnyMap& binary_a = eos["binary-a"].as(); + const UnitSystem& units = binary_a.units(); + for (auto& item2 : binary_a) { + double a0 = units.convert(item2.second, "Pa*m^6/kmol^2"); + setBinaryCoeffs(item.first, item2.first, a0); + } + } + if (foundCoeffs) { + m_coeffSource[k] = CoeffSource::EoS; + continue; + } + } + + // Coefficients have not been populated from model-specific input + double Tc = NAN, Pc = NAN, omega_ac = NAN; + if (data.hasKey("critical-parameters")) { + // Use critical state information stored in the species entry to + // calculate a, b, and the acentric factor. + auto& critProps = data["critical-parameters"].as(); + Tc = critProps.convert("critical-temperature", "K"); + Pc = critProps.convert("critical-pressure", "Pa"); + omega_ac = critProps["acentric-factor"].asDouble(); + m_coeffSource[k] = CoeffSource::CritProps; + } else { + // Search 'crit-properties.yaml' to find Tc and Pc. Load data if needed. + if (critPropsDb.empty()) { + critPropsDb = AnyMap::fromYamlFile("critical-properties.yaml"); + dbSpecies = critPropsDb["species"].asMap("name"); + } + + // All names in critical-properties.yaml are upper case + auto ucName = boost::algorithm::to_upper_copy(item.first); + if (dbSpecies.count(ucName)) { + auto& spec = *dbSpecies.at(ucName); + auto& critProps = spec["critical-parameters"].as(); + Tc = critProps.convert("critical-temperature", "K"); + Pc = critProps.convert("critical-pressure", "Pa"); + omega_ac = critProps["acentric-factor"].asDouble(); + m_coeffSource[k] = CoeffSource::Database; + } + } + + // Check if critical properties were found in either location + if (!isnan(Tc)) { + double a = omega_a * std::pow(GasConstant * Tc, 2) / Pc; + double b = omega_b * GasConstant * Tc / Pc; + setSpeciesCoeffs(item.first, a, b, omega_ac); + } else { + throw InputFileError("SoaveRedlichKwong::initThermo", data, + "No Soave-Redlich-Kwong model parameters or critical properties found for " + "species '{}'", item.first); + } + } +} + +void SoaveRedlichKwong::getSpeciesParameters(const std::string& name, + AnyMap& speciesNode) const +{ + MixtureFugacityTP::getSpeciesParameters(name, speciesNode); + size_t k = speciesIndex(name); + checkSpeciesIndex(k); + + // Pure species parameters + if (m_coeffSource[k] == CoeffSource::EoS) { + auto& eosNode = speciesNode["equation-of-state"].getMapWhere( + "model", "Soave-Redlich-Kwong", true); + eosNode["a"].setQuantity(m_a_coeffs(k, k), "Pa*m^6/kmol^2"); + eosNode["b"].setQuantity(m_b_coeffs[k], "m^3/kmol"); + eosNode["acentric-factor"] = m_acentric[k]; + } else if (m_coeffSource[k] == CoeffSource::CritProps) { + auto& critProps = speciesNode["critical-parameters"]; + double Tc = speciesCritTemperature(m_a_coeffs(k, k), m_b_coeffs[k]); + double Pc = omega_b * GasConstant * Tc / m_b_coeffs[k]; + critProps["critical-temperature"].setQuantity(Tc, "K"); + critProps["critical-pressure"].setQuantity(Pc, "Pa"); + critProps["acentric-factor"] = m_acentric[k]; + } + // Nothing to do in the case where the parameters are from the database + + if (m_binaryParameters.count(name)) { + // Include binary parameters regardless of where the pure species parameters + // were found + auto& eosNode = speciesNode["equation-of-state"].getMapWhere( + "model", "Soave-Redlich-Kwong", true); + AnyMap bin_a; + for (const auto& item : m_binaryParameters.at(name)) { + bin_a[item.first].setQuantity(item.second, "Pa*m^6/kmol^2"); + } + eosNode["binary-a"] = std::move(bin_a); + } +} + +double SoaveRedlichKwong::sresid() const +{ // Replace + double molarV = molarVolume(); + double hh = m_b / molarV; + double zz = z(); + double alpha_1 = daAlpha_dT(); + double vpb = molarV + (1.0 + Sqrt2) * m_b; + double vmb = molarV + (1.0 - Sqrt2) * m_b; + double fac = alpha_1 / (2.0 * Sqrt2 * m_b); + double sresid_mol_R = log(zz*(1.0 - hh)) + fac * log(vpb / vmb) / GasConstant; + return GasConstant * sresid_mol_R; +} + +double SoaveRedlichKwong::hresid() const +{ // Replace + double molarV = molarVolume(); + double zz = z(); + double aAlpha_1 = daAlpha_dT(); + double T = temperature(); + double vpb = molarV + (1 + Sqrt2) * m_b; + double vmb = molarV + (1 - Sqrt2) * m_b; + double fac = 1 / (2.0 * Sqrt2 * m_b); + return GasConstant * T * (zz - 1.0) + + fac * log(vpb / vmb) * (T * aAlpha_1 - m_aAlpha_mix); +} + +double SoaveRedlichKwong::liquidVolEst(double T, double& presGuess) const +{ + double v = m_b * 1.1; + double atmp, btmp, aAlphatmp; + calculateAB(atmp, btmp, aAlphatmp); + double pres = std::max(psatEst(T), presGuess); + double Vroot[3]; + bool foundLiq = false; + int m = 0; + while (m < 100 && !foundLiq) { + int nsol = solveCubic(T, pres, atmp, btmp, aAlphatmp, Vroot); + if (nsol == 1 || nsol == 2) { + double pc = critPressure(); + if (pres > pc) { + foundLiq = true; + } + pres *= 1.04; + } else { + foundLiq = true; + } + } + + if (foundLiq) { + v = Vroot[0]; + presGuess = pres; + } else { + v = -1.0; + } + return v; +} + +double SoaveRedlichKwong::densityCalc(double T, double presPa, int phaseRequested, + double rhoGuess) +{ + // It's necessary to set the temperature so that m_aAlpha_mix is set correctly. + setTemperature(T); + double tcrit = critTemperature(); + double mmw = meanMolecularWeight(); + if (rhoGuess == -1.0) { + if (phaseRequested >= FLUID_LIQUID_0) { + double lqvol = liquidVolEst(T, presPa); + rhoGuess = mmw / lqvol; + } + } else { + // Assume the Gas phase initial guess, if nothing is specified to the routine + rhoGuess = presPa * mmw / (GasConstant * T); + } + + double volGuess = mmw / rhoGuess; + m_NSolns = solveCubic(T, presPa, m_a, m_b, m_aAlpha_mix, m_Vroot); + + double molarVolLast = m_Vroot[0]; + if (m_NSolns >= 2) { + if (phaseRequested >= FLUID_LIQUID_0) { + molarVolLast = m_Vroot[0]; + } else if (phaseRequested == FLUID_GAS || phaseRequested == FLUID_SUPERCRIT) { + molarVolLast = m_Vroot[2]; + } else { + if (volGuess > m_Vroot[1]) { + molarVolLast = m_Vroot[2]; + } else { + molarVolLast = m_Vroot[0]; + } + } + } else if (m_NSolns == 1) { + if (phaseRequested == FLUID_GAS || phaseRequested == FLUID_SUPERCRIT + || phaseRequested == FLUID_UNDEFINED) + { + molarVolLast = m_Vroot[0]; + } else { + return -2.0; + } + } else if (m_NSolns == -1) { + if (phaseRequested >= FLUID_LIQUID_0 || phaseRequested == FLUID_UNDEFINED + || phaseRequested == FLUID_SUPERCRIT) + { + molarVolLast = m_Vroot[0]; + } else if (T > tcrit) { + molarVolLast = m_Vroot[0]; + } else { + return -2.0; + } + } else { + molarVolLast = m_Vroot[0]; + return -1.0; + } + return mmw / molarVolLast; +} + +double SoaveRedlichKwong::densSpinodalLiquid() const +{ + double Vroot[3]; + double T = temperature(); + int nsol = solveCubic(T, pressure(), m_a, m_b, m_aAlpha_mix, Vroot); + if (nsol != 3) { + return critDensity(); + } + + auto resid = [this, T](double v) { + double pp; + return dpdVCalc(T, v, pp); + }; + + boost::uintmax_t maxiter = 100; + std::pair vv = bmt::toms748_solve( + resid, Vroot[0], Vroot[1], bmt::eps_tolerance(48), maxiter); + + double mmw = meanMolecularWeight(); + return mmw / (0.5 * (vv.first + vv.second)); +} + +double SoaveRedlichKwong::densSpinodalGas() const +{ + double Vroot[3]; + double T = temperature(); + int nsol = solveCubic(T, pressure(), m_a, m_b, m_aAlpha_mix, Vroot); + if (nsol != 3) { + return critDensity(); + } + + auto resid = [this, T](double v) { + double pp; + return dpdVCalc(T, v, pp); + }; + + boost::uintmax_t maxiter = 100; + std::pair vv = bmt::toms748_solve( + resid, Vroot[1], Vroot[2], bmt::eps_tolerance(48), maxiter); + + double mmw = meanMolecularWeight(); + return mmw / (0.5 * (vv.first + vv.second)); +} + +double SoaveRedlichKwong::dpdVCalc(double T, double molarVol, double& presCalc) const +{ + double denom = molarVol * (molarVol + m_b); + double vmb = molarVol - m_b; + return -GasConstant * T / (vmb * vmb) + m_aAlpha_mix * (2 * molarVol + m_b) / (denom * denom); +} + +double SoaveRedlichKwong::isothermalCompressibility() const +{ + calculatePressureDerivatives(); + return -1 / (molarVolume() * m_dpdV); +} + +double SoaveRedlichKwong::thermalExpansionCoeff() const +{ + calculatePressureDerivatives(); + return -m_dpdT / (molarVolume() * m_dpdV); +} + +void SoaveRedlichKwong::calculatePressureDerivatives() const +{ + double T = temperature(); + double mv = molarVolume(); + double pres; + + m_dpdV = dpdVCalc(T, mv, pres); + m_dpdT = GasConstant / (mv - m_b) - daAlpha_dT() / (mv * (mv + m_b)); +} + +void SoaveRedlichKwong::updateMixingExpressions() +{ + // Same as PengRobinson::updateMixingExpressions() + 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(m_a,m_b,m_aAlpha_mix); +} + +void SoaveRedlichKwong::calculateAB(double& aCalc, double& bCalc, double& aAlphaCalc) const +{ + // Same as PengRobinson::calculateAB() + bCalc = 0.0; + aCalc = 0.0; + aAlphaCalc = 0.0; + for (size_t i = 0; i < m_kk; i++) { + bCalc += moleFractions_[i] * m_b_coeffs[i]; + for (size_t j = 0; j < m_kk; j++) { + aCalc += m_a_coeffs(i, j) * moleFractions_[i] * moleFractions_[j]; + aAlphaCalc += m_aAlpha_binary(i, j) * moleFractions_[i] * moleFractions_[j]; + } + } +} + +double SoaveRedlichKwong::daAlpha_dT() const +{ + // Same as PengRobinson::daAlpha_dT() + 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(temperature() / 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 SoaveRedlichKwong::d2aAlpha_dT2() const +{ + // Same as PengRobinson::d2aAlpha_dT2() + for (size_t i = 0; i < m_kk; i++) { + double Tcrit_i = speciesCritTemperature(m_a_coeffs(i, i), m_b_coeffs[i]); + double sqt_Tr = sqrt(temperature() / Tcrit_i); + double coeff1 = 1 / (Tcrit_i*sqt_Tr); + double coeff2 = sqt_Tr - 1; + // Calculate first and second derivatives of alpha for individual species + double k = m_kappa[i]; + m_dalphadT[i] = coeff1 * (k*k*coeff2 - k); + m_d2alphadT2[i] = (k*k + k) * coeff1 / (2*sqt_Tr*sqt_Tr*Tcrit_i); + } + + // Calculate mixture derivative + double d2aAlphadT2 = 0.0; + for (size_t i = 0; i < m_kk; i++) { + double alphai = m_alpha[i]; + for (size_t j = 0; j < m_kk; j++) { + double alphaj = m_alpha[j]; + double alphaij = alphai * alphaj; + double term1 = m_d2alphadT2[i] / alphai + m_d2alphadT2[j] / alphaj; + double term2 = 2 * m_dalphadT[i] * m_dalphadT[j] / alphaij; + double term3 = m_dalphadT[i] / alphai + m_dalphadT[j] / alphaj; + d2aAlphadT2 += 0.5 * moleFractions_[i] * moleFractions_[j] + * m_aAlpha_binary(i, j) + * (term1 + term2 - 0.5 * term3 * term3); + } + } + return d2aAlphadT2; +} + +void SoaveRedlichKwong::calcCriticalConditions(double& pc, double& tc, double& vc) const +{ + if (m_b <= 0.0) { + tc = 1000000.; + pc = 1.0E13; + vc = omega_vc * GasConstant * tc / pc; + return; + } + if (m_a <= 0.0) { + tc = 0.0; + pc = 0.0; + vc = 2.0 * m_b; + return; + } + tc = m_a * omega_b / (m_b * omega_a * GasConstant); + pc = omega_b * GasConstant * tc / m_b; + vc = omega_vc * GasConstant * tc / pc; +} + +int SoaveRedlichKwong::solveCubic(double T, double pres, double a, double b, double aAlpha, + double Vroot[3]) const +{ + double an = 1.0; + double bn = - GasConstant * T / pres; + double cn = (aAlpha - b * GasConstant * T) / pres - b * b; + double dn = -aAlpha * b / pres; + + double tc = a * omega_b / (b * omega_a * GasConstant); + double pc = omega_b * GasConstant * tc / b; + double vc = omega_vc * GasConstant * tc / pc; + + return MixtureFugacityTP::solveCubic(T, pres, a, b, aAlpha, Vroot, + an, bn, cn, dn, tc, vc); +} + +} diff --git a/src/thermo/ThermoFactory.cpp b/src/thermo/ThermoFactory.cpp index fd6edb4c92..7a477d368e 100644 --- a/src/thermo/ThermoFactory.cpp +++ b/src/thermo/ThermoFactory.cpp @@ -25,6 +25,7 @@ #include "cantera/thermo/PureFluidPhase.h" #include "cantera/thermo/RedlichKwongMFTP.h" #include "cantera/thermo/PengRobinson.h" +#include "cantera/thermo/SoaveRedlichKwong.h" #include "cantera/thermo/SurfPhase.h" #include "cantera/thermo/EdgePhase.h" #include "cantera/thermo/MetalPhase.h" @@ -104,6 +105,7 @@ ThermoFactory::ThermoFactory() reg("binary-solution-tabulated", []() { return new BinarySolutionTabulatedThermo(); }); addAlias("binary-solution-tabulated", "BinarySolutionTabulatedThermo"); reg("Peng-Robinson", []() { return new PengRobinson(); }); + reg("Soave-Redlich-Kwong", []() { return new SoaveRedlichKwong(); }); } ThermoPhase* ThermoFactory::newThermoPhase(const std::string& model) diff --git a/test/data/co2_SRK_example.yaml b/test/data/co2_SRK_example.yaml new file mode 100644 index 0000000000..9a1ce6ee34 --- /dev/null +++ b/test/data/co2_SRK_example.yaml @@ -0,0 +1,129 @@ +units: {length: cm, quantity: mol, activation-energy: cal/mol} + +phases: +- name: CO2-PR + species: [CO2, H2O, H2, CO, CH4, O2, N2] + thermo: Soave-Redlich-Kwong + kinetics: gas + reactions: all + state: {T: 300, P: 1 atm, mole-fractions: {CO2: 0.99, H2: 0.01}} + + +species: +- name: CO2 + composition: {C: 1, O: 2} + thermo: + model: NASA7 + temperature-ranges: [200.0, 1000.0, 3500.0] + data: + - [2.35677352, 8.98459677e-03, -7.12356269e-06, 2.45919022e-09, -1.43699548e-13, + -4.83719697e+04, 9.90105222] + - [3.85746029, 4.41437026e-03, -2.21481404e-06, 5.23490188e-10, -4.72084164e-14, + -4.8759166e+04, 2.27163806] + note: L7/88 + equation-of-state: + model: Soave-Redlich-Kwong + a: 3.700552E+11 + b: 29.6547 + acentric-factor: 0.228 +- name: H2O + composition: {H: 2, O: 1} + thermo: + model: NASA7 + temperature-ranges: [200.0, 1000.0, 3500.0] + data: + - [4.19864056, -2.0364341e-03, 6.52040211e-06, -5.48797062e-09, 1.77197817e-12, + -3.02937267e+04, -0.849032208] + - [3.03399249, 2.17691804e-03, -1.64072518e-07, -9.7041987e-11, 1.68200992e-14, + -3.00042971e+04, 4.9667701] + note: L8/89 + equation-of-state: + model: Soave-Redlich-Kwong + a: 5.608487E+11 + b: 21.1282 + acentric-factor: 0.344 +- name: H2 + composition: {H: 2} + thermo: + model: NASA7 + temperature-ranges: [200.0, 1000.0, 3500.0] + data: + - [2.34433112, 7.98052075e-03, -1.9478151e-05, 2.01572094e-08, -7.37611761e-12, + -917.935173, 0.683010238] + - [3.3372792, -4.94024731e-05, 4.99456778e-07, -1.79566394e-10, 2.00255376e-14, + -950.158922, -3.20502331] + note: TPIS78 + equation-of-state: + model: Soave-Redlich-Kwong + a: 2.494771E+10 + b: 18.4290 + acentric-factor: -0.22 +- name: CO + composition: {C: 1, O: 1} + thermo: + model: NASA7 + temperature-ranges: [200.0, 1000.0, 3500.0] + data: + - [3.57953347, -6.1035368e-04, 1.01681433e-06, 9.07005884e-10, -9.04424499e-13, + -1.4344086e+04, 3.50840928] + - [2.71518561, 2.06252743e-03, -9.98825771e-07, 2.30053008e-10, -2.03647716e-14, + -1.41518724e+04, 7.81868772] + note: TPIS79 + equation-of-state: + model: Soave-Redlich-Kwong + a: 1.502575E+11 + b: 27.4578 + acentric-factor: 0.049 +- name: CH4 + composition: {C: 1, H: 4} + thermo: + model: NASA7 + temperature-ranges: [200.0, 1000.0, 3500.0] + data: + - [5.14987613, -0.0136709788, 4.91800599e-05, -4.84743026e-08, 1.66693956e-11, + -1.02466476e+04, -4.64130376] + - [0.074851495, 0.0133909467, -5.73285809e-06, 1.22292535e-09, -1.0181523e-13, + -9468.34459, 18.437318] + note: L8/88 + equation-of-state: + model: Soave-Redlich-Kwong + a: 2.333891E+11 + b: 29.8499 + acentric-factor: 0.01 +- name: O2 + composition: {O: 2} + thermo: + model: NASA7 + temperature-ranges: [200.0, 1000.0, 3500.0] + data: + - [3.78245636, -2.99673416e-03, 9.84730201e-06, -9.68129509e-09, 3.24372837e-12, + -1063.94356, 3.65767573] + - [3.28253784, 1.48308754e-03, -7.57966669e-07, 2.09470555e-10, -2.16717794e-14, + -1088.45772, 5.45323129] + note: TPIS89 + equation-of-state: + model: Soave-Redlich-Kwong + a: 1.400265E+11 + b: 22.0823 + acentric-factor: 0.022 +- name: N2 + composition: {N: 2} + thermo: + model: NASA7 + temperature-ranges: [300.0, 1000.0, 5000.0] + data: + - [3.298677, 1.4082404e-03, -3.963222e-06, 5.641515e-09, -2.444854e-12, + -1020.8999, 3.950372] + - [2.92664, 1.4879768e-03, -5.68476e-07, 1.0097038e-10, -6.753351e-15, + -922.7977, 5.980528] + note: '121286' + equation-of-state: + model: Soave-Redlich-Kwong + a: 1.388390E+11 + b: 31.2734 + acentric-factor: 0.04 + + +reactions: +- equation: CO2 + H2 <=> CO + H2O # Reaction 1 + rate-constant: {A: 1.2E+3, b: 0, Ea: 0} diff --git a/test/data/consistency-cases.yaml b/test/data/consistency-cases.yaml index dd06e62482..2e2761a16c 100644 --- a/test/data/consistency-cases.yaml +++ b/test/data/consistency-cases.yaml @@ -55,6 +55,15 @@ peng-robinson: - {T: 320, P: 200 bar, X: {CO2: 0.3, CH4: 0.5, H2O: 0.2}} - {T: 600, P: 200 bar, X: {CO2: 0.4, CH4: 0.2, H2O: 0.4}} +soave-redlich-kwong: + setup: + file: co2_SRK_example.yaml + rtol_fd: 1e-5 + states: + - {T: 300, P: 101325, X: {CO2: 0.7, CH4: 0.2, H2O: 0.1}} + - {T: 320, P: 200 bar, X: {CO2: 0.3, CH4: 0.5, H2O: 0.2}} + - {T: 600, P: 200 bar, X: {CO2: 0.4, CH4: 0.2, H2O: 0.4}} + ideal-molal-solution: setup: file: thermo-models.yaml diff --git a/test/thermo/consistency.cpp b/test/thermo/consistency.cpp index 0004d24212..59e0c0fc93 100644 --- a/test/thermo/consistency.cpp +++ b/test/thermo/consistency.cpp @@ -345,6 +345,56 @@ TEST_P(TestConsistency, dSdv_const_T_eq_dPdT_const_V) { } } +TEST_P(TestConsistency, betaT_eq_minus_dmv_dP_const_T_div_mv) +{ + double betaT1; + try { + betaT1 = phase->isothermalCompressibility(); + } catch (NotImplementedError& err) { + GTEST_SKIP() << err.getMethod() << " threw NotImplementedError"; + } + + double T = phase->temperature(); + double P1 = phase->pressure(); + double mv1 = phase->molarVolume(); + + double P2 = P1 * (1 + 1e-7); + phase->setState_TP(T, P2); + double betaT2 = phase->isothermalCompressibility(); + double mv2 = phase->molarVolume(); + + double betaT_mid = 0.5 * (betaT1 + betaT2); + double mv_mid = 0.5 * (mv1 + mv2); + double betaT_fd = -1 / mv_mid * (mv2 - mv1) / (P2 - P1); + + EXPECT_NEAR(betaT_fd, betaT_mid, max({rtol_fd * betaT_mid, rtol_fd * betaT_fd, atol})); +} + +TEST_P(TestConsistency, alphaV_eq_dmv_dT_const_P_div_mv) +{ + double alphaV1; + try { + alphaV1 = phase->thermalExpansionCoeff(); + } catch (NotImplementedError& err) { + GTEST_SKIP() << err.getMethod() << " threw NotImplementedError"; + } + + double P = phase->pressure(); + double T1 = phase->temperature(); + double mv1 = phase->molarVolume(); + + double T2 = T1 * (1 + 1e-7); + phase->setState_TP(T2, P); + double alphaV2 = phase->thermalExpansionCoeff(); + double mv2 = phase->molarVolume(); + + double alphaV_mid = 0.5 * (alphaV1 + alphaV2); + double mv_mid = 0.5 * (mv1 + mv2); + double alphaV_fd = 1 / mv_mid * (mv2 - mv1) / (T2 - T1); + + EXPECT_NEAR(alphaV_fd, alphaV_mid, max({rtol_fd * alphaV_mid, rtol_fd * alphaV_fd, atol})); +} + // ---------- Tests for consistency of standard state properties --------------- TEST_P(TestConsistency, hk0_eq_uk0_plus_p_vk0) @@ -556,6 +606,12 @@ INSTANTIATE_TEST_SUITE_P(PengRobinson, TestConsistency, testing::ValuesIn(getStates("peng-robinson"))) ); +INSTANTIATE_TEST_SUITE_P(SoaveRedlichKwong, TestConsistency, + testing::Combine( + testing::Values(getSetup("soave-redlich-kwong")), + testing::ValuesIn(getStates("soave-redlich-kwong"))) +); + INSTANTIATE_TEST_SUITE_P(IdealMolalSolution, TestConsistency, testing::Combine( testing::Values(getSetup("ideal-molal-solution")),