diff --git a/include/cantera/thermo/PengRobinsonMFTP.h b/include/cantera/thermo/PengRobinsonMFTP.h new file mode 100644 index 00000000000..2d476742649 --- /dev/null +++ b/include/cantera/thermo/PengRobinsonMFTP.h @@ -0,0 +1,430 @@ +//! @file PengRobinsonMFTP.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_PENGROBINSONMFTP_H +#define CT_PENGROBINSONMFTP_H + +#include "MixtureFugacityTP.h" +#include "cantera/base/Array.h" + +namespace Cantera +{ +/** + * Implementation of a multi-species Peng-Robinson equation of state + * + * @ingroup thermoprops + */ +class PengRobinsonMFTP : public MixtureFugacityTP +{ +public: + //! @name Constructors and Duplicators + //! @{ + + //! Base constructor. + PengRobinsonMFTP(); + + //! Construct and initialize a PengRobinsonMFTP object directly from an + //! ASCII input file + /*! + * @param infile Name of the input file containing the phase XML data + * to set up the object + * @param id ID of the phase in the input file. Defaults to the empty + * string. + */ + PengRobinsonMFTP(const std::string& infile, const std::string& id=""); + + //! Construct and initialize a PengRobinsonMFTP object directly from an + //! XML database + /*! + * @param phaseRef XML phase node containing the description of the phase + * @param id id attribute containing the name of the phase. (default + * is the empty string) + */ + PengRobinsonMFTP(XML_Node& phaseRef, const std::string& id = ""); + + virtual std::string type() const { + return "PengRobinson"; + } + + //! @name Molar Thermodynamic properties + //! @{ + + virtual double enthalpy_mole() const; + virtual double entropy_mole() const; + 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(\alpha a\right)_{mix}}{v^2 + 2b_{mix}v - b_{mix}^2 } + * \f] + * + * where: + * + * \f[ + * \alpha = \left[ 1 + \kappa \left(1-T_r^{0.5}\right)\right]^2 + * \f] + * + * and + * + * \f[ + * \kappa = \left(0.37464 + 1.54226\omega - 0.26992\omega^2\right) for omega <= 0.491 + * \kappa = \left(0.379642 + 1.487503\omega - 0.164423\omega^2 + 0.016667\omega^3 \right) for omega > 0.491 + * \f] + * + *Coefficients a_mix, b_mix and (a \alpha)_{mix} are caclulated 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; + + // @} + +protected: + /** + * Calculate the density of the mixture using the partial molar volumes and + * mole fractions as input + * + * The formula for this is + * + * \f[ + * \rho = \frac{\sum_k{X_k W_k}}{\sum_k{X_k V_k}} + * \f] + * + * where \f$X_k\f$ are the mole fractions, \f$W_k\f$ are the molecular + * weights, and \f$V_k\f$ are the pure species molar volumes. + * + * Note, the basis behind this formula is that in an ideal solution the + * partial molar volumes are equal to the species standard state molar + * volumes. The species molar volumes may be functions of temperature and + * pressure. + */ + virtual void calcDensity(); + + virtual void setTemperature(const double temp); + virtual void compositionChanged(); + +public: + virtual void getActivityConcentrations(double* c) 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. In many cases, this + * quantity will be the same for all species in a phase. + * 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 m3 kmol-1. + */ + 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 + //@{ + + //! Get the array of non-dimensional species chemical potentials. + //! These are partial molar Gibbs free energies. + /*! + * \f$ \mu_k / \hat R T \f$. + * Units: unitless + * + * We close the loop on this function here calling getChemPotentials() and + * then dividing by RT. No need for child classes to handle. + * + * @param mu Output vector of non-dimensional species chemical potentials + * Length: m_kk. + */ + virtual void getChemPotentials_RT(double* mu) const; + + virtual void getChemPotentials(double* mu) const; + virtual void getPartialMolarEnthalpies(double* hbar) const; + virtual void getPartialMolarEntropies(double* sbar) const; + virtual void getPartialMolarIntEnergies(double* ubar) const; + virtual void getPartialMolarCp(double* cpbar) const; + virtual void getPartialMolarVolumes(double* vbar) const; + + //! Calculate the temperature dependent interaction parameter alpha needed for P-R EoS + /* + * The temperature dependent parameter in P-R EoS is calculated as + * \alpha = [1 + \kappa(1 - sqrt{T/T_crit}]^2 + * kappa is a function calulated based on the accentric factor. + * Units: unitless + */ + virtual void calculateAlpha(const std::string& species, double a, double b, double w); + //@} + /// @name Critical State Properties. + //@{ + + virtual double critTemperature() const; + virtual double critPressure() const; + virtual double critVolume() const; + virtual double critCompressibility() const; + virtual double critDensity() const; + virtual double speciesCritTemperature(double a, double b) const; + +public: + //@} + //! @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 setParametersFromXML(const XML_Node& thermoNode); + virtual void setToEquilState(const double* lambda_RT); + virtual void initThermoXML(XML_Node& phaseNode, const std::string& id); + + //! Retrieve a and b coefficients by looking up tabulated critical parameters + /*! + * If pureFluidParameters are not provided for any species in the phase, + * consult the critical properties tabulated in build/data/thermo/critProperties.xml. + * If the species is found there, calculate pure fluid parameters a_k and b_k as: + * \f[ a_k = 0.4278*R**2*T_c^2/P_c \f] + * + * and: + * \f[ b_k = 0.08664*R*T_c/P_c \f] + * + * @param iName Name of the species + */ + virtual std::vector getCoeff(const std::string& iName); + + //! Set the pure fluid interaction parameters for a species + /*! + * The "a" parameter for species *i* in the Peng-Robinson model is assumed + * to be a linear function of temperature: + * \f[ a = a_0 + a_1 T \f] + * + * @param species Name of the species + * @param a0 constant term in the expression for the "a" parameter + * of the specified species [Pa-m^6/kmol^2] + * @param a1 temperature-proportional term in the expression for the + * "a" parameter of the specified species [Pa-m^6/kmol^2/K] + * @param b "b" parameter in the Peng-Robinson model [m^3/kmol] + * @param alpha dimensionless function of T_r and \omega + * @param omega acentric factor + */ + void setSpeciesCoeffs(const std::string& species, double a, double b, + double w); + + //! Set values for the interaction parameter between two species + /*! + * The "a" parameter for interactions between species *i* and *j* is + * assumed by default to be computed as: + * \f[ a_{ij} = \sqrt(a_{i, 0} a_{j, 0}) + \sqrt(a_{i, 1} a_{j, 1}) T \f] + * + * This function overrides the defaults with the specified parameters: + * \f[ a_{ij} = a_{ij, 0} + a_{ij, 1} T \f] + * + * @param species_i Name of one species + * @param species_j Name of the other species + * @param a0 constant term in the "a" expression [Pa-m^6/kmol^2] + * @param a1 temperature-proportional term in the "a" expression + * [Pa-m^6/kmol^2/K] + */ + void setBinaryCoeffs(const std::string& species_i, + const std::string& species_j, double a0, double a1); + +private: + //! Read the pure species PengRobinson input parameters + /*! + * @param pureFluidParam XML_Node for the pure fluid parameters + */ + void readXMLPureFluid(XML_Node& pureFluidParam); + + //! Read the cross species PengRobinson input parameters + /*! + * @param crossFluidParam XML_Node for the cross fluid parameters + */ + void readXMLCrossFluid(XML_Node& crossFluidParam); + + // @} + +protected: + // Special functions inherited from MixtureFugacityTP + virtual double sresid() const; + virtual double hresid() const; + +public: + virtual double liquidVolEst(double TKelvin, double& pres) const; + virtual double densityCalc(double TKelvin, double pressure, int phase, double rhoguess); + + virtual double densSpinodalLiquid() const; + virtual double densSpinodalGas() const; + virtual double pressureCalc(double TKelvin, double molarVol) const; + virtual double dpdVCalc(double TKelvin, double molarVol, double& presCalc) const; + + //! Calculate dpdV and dpdT at the current conditions + /*! + * These are stored internally. + */ + void pressureDerivatives() const; + + virtual void updateMixingExpressions(); + + //! Update the a and b parameters + /*! + * The a and the b parameters depend on the mole fraction and the + * temperature. This function updates the internal numbers based on the + * state of the object. + */ + void updateAB(); + + //! Calculate the a and the b 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 temp Temperature (TKelvin) + * @param aCalc (output) Returns the a value + * @param bCalc (output) Returns the b value. + */ + void calculateAB(double temp, double& aCalc, double& bCalc, double& aAlpha) const; + + // Special functions not inherited from MixtureFugacityTP + + double daAlpha_dT() const; + double d2aAlpha_dT2() const; + + void calcCriticalConditions(double a, double b,double& pc, double& tc, double& vc) const; + + //! Solve the cubic equation of state + /*! + * The P-R equation of state may be solved via the following formula: + * + * V**3 - V**2(RT/P - b) - V(2bRT/P - \alpha a/P + 3*b*b) - (a \alpha b/p - b*b RT/P - b*b*b) = 0 + * + * Returns the number of solutions found. If it only finds the liquid + * branch solution, it will return a -1 or a -2 instead of 1 or 2. If it + * returns 0, then there is an error. + * The cubic equation is solved using Nickall's method (Ref: The Mathematical Gazette(1993), 77(November), 354–359, https://www.jstor.org/stable/3619777) + */ + int NicholsSolve(double TKelvin, double pres, double a, double b, double aAlpha, + double Vroot[3]) const; + +protected: + //! Form of the temperature parameterization + /*! + * 0 = There is no temperature parameterization of a or b + * 1 = The a_ij parameter is a linear function of the temperature + */ + int m_formTempParam; + + //! Value of b in the equation of state + /*! + * m_b_current is a function of the temperature and the mole fractions. + */ + double m_b_current; + + //! Value of a and alpha in the equation of state + /*! + * m_aAlpha_current is a function of the temperature and the mole fractions. m_a_current depends only on the mole fractions. + */ + double m_a_current; + double m_aAlpha_current; + + // Vectors required to store a_coeff, b_coeff, alpha, kappa and other values for every species. Length = m_kk + vector_fp a_vec_Curr_; + vector_fp b_vec_Curr_; + vector_fp aAlpha_vec_Curr_; + vector_fp alpha_vec_Curr_; + vector_fp kappa_vec_; + mutable vector_fp dalphadT_vec_Curr_; + mutable vector_fp d2alphadT2_; + + Array2D a_coeff_vec; + Array2D aAlpha_coeff_vec; + + int NSolns_; + + double Vroot_[3]; + + //! Temporary storage - length = m_kk. + mutable vector_fp m_pp; + + //! Temporary storage - length = m_kk. + mutable vector_fp m_tmpV; + + // 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 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 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 dpdni_; + +public: + //! Omega constants: a0 (= omega_a) and b0 (= omega_b) values used in Peng-Robinson equation of state + /*! + * These values are calculated by solving P-R cubic equation at the critical point. + */ + static const double omega_a; + + //! Omega constant for b + static const double omega_b; + + //! Omega constant for the critical molar volume + static const double omega_vc; +}; +} + +#endif diff --git a/interfaces/cython/cantera/ctml_writer.py b/interfaces/cython/cantera/ctml_writer.py index 26cba15b02e..b0fa5abb313 100644 --- a/interfaces/cython/cantera/ctml_writer.py +++ b/interfaces/cython/cantera/ctml_writer.py @@ -835,23 +835,27 @@ class pureFluidParameters(activityCoefficients): """ """ - def __init__(self, species = None, a_coeff = [], b_coeff = 0): + def __init__(self, species = None, a_coeff = [], b_coeff = 0, acentric_factor = None): """ """ self._species = species self._acoeff = a_coeff self._bcoeff = b_coeff + self._w_ac = acentric_factor def build(self,a): f= a.addChild("pureFluidParameters") f['species'] = self._species s = '%.10g, %.10g\n' % (self._acoeff[0], self._acoeff[1]) - ac = f.addChild("a_coeff",s) + ac = f.addChild("a_coeff", s) ac["units"] = _upres+'-'+_ulen+'6/'+_umol+'2' ac["model"] = "linear_a" s = '%.10g\n' % self._bcoeff - bc = f.addChild("b_coeff",s) + bc = f.addChild("b_coeff", s) bc["units"] = _ulen+'3/'+_umol + if self._w_ac: + s = '%.10g\n' % self._w_ac + cc = f.addChild("acentric_factor", s) class crossFluidParameters(activityCoefficients): @@ -865,12 +869,12 @@ def build(self,a): f["species2"] = self._species2 f["species1"] = self._species1 s = '%.10g, %.10g\n' % (self._acoeff[0], self._acoeff[1]) - ac = f.addChild("a_coeff",s) + ac = f.addChild("a_coeff", s) ac["units"] = _upres+'-'+_ulen+'6/'+_umol+'2' ac["model"] = "linear_a" if self._bcoeff: s = '%.10g\n' % self._bcoeff - bc = f.addChild("b_coeff",s) + bc = f.addChild("b_coeff", s) bc["units"] = _ulen+'3/'+_umol @@ -2464,6 +2468,47 @@ def build(self, p): k = ph.addChild("kinetics") k['model'] = self._kin +class PengRobinsonMFTP(phase): + """A multi-component fluid model for non-ideal gas fluids. """ + + def __init__(self, + name = '', + elements = '', + species = '', + note = '', + reactions = 'none', + kinetics = 'GasKinetics', + initial_state = None, + activity_coefficients = None, + transport = 'None', + options = []): + + phase.__init__(self, name, 3, elements, species, note, reactions, + initial_state, options) + self._pure = 0 + self._kin = kinetics + self._tr = transport + self._activityCoefficients = activity_coefficients + + def build(self, p): + ph = phase.build(self, p) + e = ph.child("thermo") + e['model'] = 'PengRobinsonMFTP' + if self._activityCoefficients: + a = e.addChild("activityCoefficients") + if isinstance(self._activityCoefficients, activityCoefficients): + self._activityCoefficients.build(a) + else: + na = len(self._activityCoefficients) + for n in range(na): + self._activityCoefficients[n].build(a) + + if self._tr: + t = ph.addChild('transport') + t['model'] = self._tr + if self._kin: + k = ph.addChild("kinetics") + k['model'] = self._kin class ideal_interface(phase): """A chemically-reacting ideal surface solution of multiple species.""" diff --git a/src/thermo/PengRobinsonMFTP.cpp b/src/thermo/PengRobinsonMFTP.cpp new file mode 100644 index 00000000000..012468e6f5c --- /dev/null +++ b/src/thermo/PengRobinsonMFTP.cpp @@ -0,0 +1,1300 @@ +//! @file PengRobinsonMFTP.cpp + +// This file is part of Cantera. See License.txt in the top-level directory or +// at http://www.cantera.org/license.txt for license and copyright information. + +#include "cantera/thermo/PengRobinsonMFTP.h" +#include "cantera/thermo/ThermoFactory.h" +#include "cantera/base/stringUtils.h" +#include "cantera/base/ctml.h" + +#include + +#define _USE_MATH_DEFINES +#include + +using namespace std; +namespace bmt = boost::math::tools; + +namespace Cantera +{ + +const double PengRobinsonMFTP::omega_a = 4.5723552892138218E-01; +const double PengRobinsonMFTP::omega_b = 7.77960739038885E-02; +const double PengRobinsonMFTP::omega_vc = 3.07401308698703833E-01; + +PengRobinsonMFTP::PengRobinsonMFTP() : + m_formTempParam(0), + m_b_current(0.0), + m_a_current(0.0), + m_aAlpha_current(0.0), + NSolns_(0), + dpdV_(0.0), + dpdT_(0.0) +{ + fill_n(Vroot_, 3, 0.0); +} + +PengRobinsonMFTP::PengRobinsonMFTP(const std::string& infile, const std::string& id_) : + m_formTempParam(0), + m_b_current(0.0), + m_a_current(0.0), + m_aAlpha_current(0.0), + NSolns_(0), + dpdV_(0.0), + dpdT_(0.0) +{ + fill_n(Vroot_, 3, 0.0); + initThermoFile(infile, id_); +} + +PengRobinsonMFTP::PengRobinsonMFTP(XML_Node& phaseRefRoot, const std::string& id_) : + m_formTempParam(0), + m_b_current(0.0), + m_a_current(0.0), + m_aAlpha_current(0.0), + NSolns_(0), + dpdV_(0.0), + dpdT_(0.0) +{ + fill_n(Vroot_, 3, 0.0); + importPhase(phaseRefRoot, this); +} + +void PengRobinsonMFTP::calculateAlpha(const std::string& species, double a, double b, double w) +{ + size_t k = speciesIndex(species); + if (k == npos) { + throw CanteraError("PengRobinsonMFTP::setSpeciesCoeffs", + "Unknown species '{}'.", species); + } + + // Calculate value of kappa (independent of temperature) + // w is an acentric factor of species and must be specified in the CTI file + + if (w <= 0.491) { + kappa_vec_[k] = 0.37464 + 1.54226*w - 0.26992*w*w; + } else { + kappa_vec_[k] = 0.374642 + 1.487503*w - 0.164423*w*w + 0.016666*w*w*w; + } + + //Calculate alpha (temperature dependent interaction parameter) + double criTemp = speciesCritTemperature(a, b); // critical temperature of individual species + double sqt_T_r = sqrt(temperature() / criTemp); + double sqt_alpha = 1 + kappa_vec_[k] * (1 - sqt_T_r); + alpha_vec_Curr_[k] = sqt_alpha*sqt_alpha; +} + +void PengRobinsonMFTP::setSpeciesCoeffs(const std::string& species, + double a, double b, double w) +{ + size_t k = speciesIndex(species); + if (k == npos) { + throw CanteraError("PengRobinsonMFTP::setSpeciesCoeffs", + "Unknown species '{}'.", species); + } + size_t counter = k + m_kk * k; + a_coeff_vec(0, counter) = a; + // we store this locally because it is used below to calculate a_Alpha: + double aAlpha_k = a*alpha_vec_Curr_[k]; + aAlpha_coeff_vec(0, counter) = 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(a_coeff_vec(0, j + m_kk * j) * a); + double aAlpha_j = a*alpha_vec_Curr_[j]; + double a_Alpha = sqrt(aAlpha_j*aAlpha_k); + if (a_coeff_vec(0, j + m_kk * k) == 0) { + a_coeff_vec(0, j + m_kk * k) = a0kj; + aAlpha_coeff_vec(0, j + m_kk * k) = a_Alpha; + a_coeff_vec(0, k + m_kk * j) = a0kj; + aAlpha_coeff_vec(0, k + m_kk * j) = a_Alpha; + } + } + a_coeff_vec.getRow(0, a_vec_Curr_.data()); + aAlpha_coeff_vec.getRow(0, aAlpha_vec_Curr_.data()); + b_vec_Curr_[k] = b; +} + +void PengRobinsonMFTP::setBinaryCoeffs(const std::string& species_i, + const std::string& species_j, double a0, double alpha) +{ + size_t ki = speciesIndex(species_i); + if (ki == npos) { + throw CanteraError("PengRobinsonMFTP::setBinaryCoeffs", + "Unknown species '{}'.", species_i); + } + size_t kj = speciesIndex(species_j); + if (kj == npos) { + throw CanteraError("PengRobinsonMFTP::setBinaryCoeffs", + "Unknown species '{}'.", species_j); + } + + size_t counter1 = ki + m_kk * kj; + size_t counter2 = kj + m_kk * ki; + a_coeff_vec(0, counter1) = a_coeff_vec(0, counter2) = a0; + aAlpha_coeff_vec(0, counter1) = aAlpha_coeff_vec(0, counter2) = a0*alpha; + a_vec_Curr_[counter1] = a_vec_Curr_[counter2] = a0; + aAlpha_vec_Curr_[counter1] = aAlpha_vec_Curr_[counter2] = a0*alpha; +} + +// ------------Molar Thermodynamic Properties ------------------------- + +double PengRobinsonMFTP::enthalpy_mole() const +{ + _updateReferenceStateThermo(); + double h_ideal = RT() * mean_X(m_h0_RT); + double h_nonideal = hresid(); + return h_ideal + h_nonideal; +} + +double PengRobinsonMFTP::entropy_mole() const +{ + _updateReferenceStateThermo(); + double sr_ideal = GasConstant * (mean_X(m_s0_R) - sum_xlogx() + - std::log(pressure()/refPressure())); + double sr_nonideal = sresid(); + return sr_ideal + sr_nonideal; +} + +double PengRobinsonMFTP::cp_mole() const +{ + _updateReferenceStateThermo(); + double TKelvin = temperature(); + double mv = molarVolume(); + double vpb = mv + (1 + M_SQRT2)*m_b_current; + double vmb = mv + (1 - M_SQRT2)*m_b_current; + pressureDerivatives(); + double cpref = GasConstant * mean_X(m_cp0_R); + double dHdT_V = cpref + mv * dpdT_ - GasConstant + + 1.0 / (2.0 * M_SQRT2 *m_b_current) * log(vpb / vmb) * TKelvin *d2aAlpha_dT2(); + return dHdT_V - (mv + TKelvin * dpdT_ / dpdV_) * dpdT_; +} + +double PengRobinsonMFTP::cv_mole() const +{ + _updateReferenceStateThermo(); + double TKelvin = temperature(); + pressureDerivatives(); + return (cp_mole() + TKelvin* dpdT_* dpdT_ / dpdV_); +} + +double PengRobinsonMFTP::pressure() const +{ + _updateReferenceStateThermo(); + // Get a copy of the private variables stored in the State object + double TKelvin = temperature(); + double mv = molarVolume(); + double denom = mv * mv + 2 * mv * m_b_current - m_b_current * m_b_current; + double pp = GasConstant * TKelvin / (mv - m_b_current) - m_aAlpha_current / denom; + return pp; +} + +void PengRobinsonMFTP::calcDensity() +{ + // Calculate the molarVolume of the solution (m**3 kmol-1) + const double* const dtmp = moleFractdivMMW(); + getPartialMolarVolumes(m_tmpV.data()); + double invDens = dot(m_tmpV.begin(), m_tmpV.end(), dtmp); + + // Set the density in the parent State object directly, by calling the + // Phase::setDensity() function. + Phase::setDensity(1.0/invDens); +} + +void PengRobinsonMFTP::setTemperature(const double temp) +{ + Phase::setTemperature(temp); + _updateReferenceStateThermo(); + updateAB(); +} + +void PengRobinsonMFTP::compositionChanged() +{ + MixtureFugacityTP::compositionChanged(); + updateAB(); +} + +void PengRobinsonMFTP::getActivityConcentrations(double* c) const +{ + getActivityCoefficients(c); + double p_RT = pressure() / RT(); + for (size_t k = 0; k < m_kk; k++) { + c[k] *= moleFraction(k)* p_RT; + } +} + +double PengRobinsonMFTP::standardConcentration(size_t k) const +{ + getStandardVolumes(m_tmpV.data()); + return 1.0 / m_tmpV[k]; +} + +void PengRobinsonMFTP::getActivityCoefficients(double* ac) const +{ + double mv = molarVolume(); + double T = temperature(); + double vpb2 = mv + (1 + M_SQRT2)*m_b_current; + double vmb2 = mv + (1 - M_SQRT2)*m_b_current; + double vmb = mv - m_b_current; + double pres = pressure(); + + for (size_t k = 0; k < m_kk; k++) { + m_pp[k] = 0.0; + for (size_t i = 0; i < m_kk; i++) { + size_t counter = k + m_kk*i; + m_pp[k] += moleFractions_[i] * aAlpha_vec_Curr_[counter]; + } + } + double num = 0; + double den = 2 * M_SQRT2 * m_b_current * m_b_current; + double den2 = m_b_current*(mv * mv + 2 * mv * m_b_current - m_b_current * m_b_current); + double RTkelvin = RT(); + for (size_t k = 0; k < m_kk; k++) { + num = 2 * m_b_current * m_pp[k] - m_aAlpha_current* b_vec_Curr_[k]; + ac[k] = (-RTkelvin *log(pres*mv/ RTkelvin) + RTkelvin * log(mv / vmb) + + RTkelvin * b_vec_Curr_[k] / vmb + - (num /den) * log(vpb2/vmb2) + - m_aAlpha_current* b_vec_Curr_[k] * mv/den2 + ); + } + for (size_t k = 0; k < m_kk; k++) { + ac[k] = exp(ac[k]/ RTkelvin); + } +} + +// ---- Partial Molar Properties of the Solution ----------------- + +void PengRobinsonMFTP::getChemPotentials_RT(double* muRT) const +{ + getChemPotentials(muRT); + double RTkelvin = RT(); + for (size_t k = 0; k < m_kk; k++) { + muRT[k] *= 1.0 / RTkelvin; + } +} + +void PengRobinsonMFTP::getChemPotentials(double* mu) const +{ + getGibbs_ref(mu); + double RTkelvin = RT(); + for (size_t k = 0; k < m_kk; k++) { + double xx = std::max(SmallNumber, moleFraction(k)); + mu[k] += RTkelvin *(log(xx)); + } + + double mv = molarVolume(); + double vmb = mv - m_b_current; + double vpb2 = mv + (1 + M_SQRT2)*m_b_current; + double vmb2 = mv + (1 - M_SQRT2)*m_b_current; + + for (size_t k = 0; k < m_kk; k++) { + m_pp[k] = 0.0; + for (size_t i = 0; i < m_kk; i++) { + size_t counter = k + m_kk*i; + m_pp[k] += moleFractions_[i] * aAlpha_vec_Curr_[counter]; + } + } + double pres = pressure(); + double refP = refPressure(); + double num = 0; + double den = 2 * M_SQRT2 * m_b_current * m_b_current; + double den2 = m_b_current*(mv * mv + 2 * mv * m_b_current - m_b_current * m_b_current); + + for (size_t k = 0; k < m_kk; k++) { + num = 2 * m_b_current * m_pp[k] - m_aAlpha_current* b_vec_Curr_[k]; + + mu[k] += (RTkelvin * log(pres/refP) - RTkelvin * log(pres * mv / RTkelvin) + + RTkelvin * log(mv / vmb) + + RTkelvin * b_vec_Curr_[k] / vmb + - (num /den) * log(vpb2/vmb2) + - m_aAlpha_current* b_vec_Curr_[k] * mv/den2 + ); + } +} + +void PengRobinsonMFTP::getPartialMolarEnthalpies(double* hbar) const +{ + // First we get the reference state contributions + getEnthalpy_RT_ref(hbar); + scale(hbar, hbar+m_kk, hbar, RT()); + + // We calculate dpdni_ + double TKelvin = temperature(); + double mv = molarVolume(); + double vmb = mv - m_b_current; + double vpb2 = mv + (1 + M_SQRT2)*m_b_current; + double vmb2 = mv + (1 - M_SQRT2)*m_b_current; + + for (size_t k = 0; k < m_kk; k++) { + m_pp[k] = 0.0; + for (size_t i = 0; i < m_kk; i++) { + size_t counter = k + m_kk*i; + m_pp[k] += moleFractions_[i] * aAlpha_vec_Curr_[counter]; + } + } + + double den = mv * mv + 2 * mv * m_b_current - m_b_current * m_b_current; + double den2 = den*den; + double RTkelvin = RT(); + for (size_t k = 0; k < m_kk; k++) { + dpdni_[k] = RTkelvin /vmb + RTkelvin * b_vec_Curr_[k] / (vmb * vmb) - 2.0 * m_pp[k] / den + + 2 * vmb * m_aAlpha_current * b_vec_Curr_[k] / den2; + } + + double daAlphadT = daAlpha_dT(); + double fac = TKelvin * daAlphadT - m_aAlpha_current; + + pressureDerivatives(); + double fac2 = mv + TKelvin * dpdT_ / dpdV_; + double fac3 = 2 * M_SQRT2 * m_b_current *m_b_current; + for (size_t k = 0; k < m_kk; k++) { + double hE_v = mv * dpdni_[k] - RTkelvin + (2 * m_b_current - b_vec_Curr_[k]) / fac3 * log(vpb2 / vmb2)*fac + + (mv * b_vec_Curr_[k]) /(m_b_current*den) * fac; + hbar[k] = hbar[k] + hE_v; + hbar[k] -= fac2 * dpdni_[k]; + } +} + +void PengRobinsonMFTP::getPartialMolarEntropies(double* sbar) const +{ + getEntropy_R_ref(sbar); + scale(sbar, sbar+m_kk, sbar, GasConstant); + double TKelvin = temperature(); + double mv = molarVolume(); + double vmb = mv - m_b_current; + double vpb2 = mv + (1 + M_SQRT2)*m_b_current; + double vmb2 = mv + (1 - M_SQRT2)*m_b_current; + double refP = refPressure(); + double daAlphadT = daAlpha_dT(); + double coeff1 = 0; + double den1 = 2 * M_SQRT2 * m_b_current * m_b_current; + double den2 = mv * mv + 2 * mv * m_b_current - m_b_current * m_b_current; + + // Calculate sum(n_j (a alpha)_i, k * (1/alpha_k d/dT(alpha_k))) -> m_pp + // Calculate sum(n_j (a alpha)_i, k * (1/alpha_i d/dT(alpha_i))) -> m_tmpV + for (size_t k = 0; k < m_kk; k++) { + m_pp[k] = 0.0; + m_tmpV[k] = 0; + for (size_t i = 0; i < m_kk; i++) { + size_t counter = k + m_kk*i; + m_pp[k] += moleFractions_[i] * aAlpha_vec_Curr_[counter]; + m_tmpV[k] += moleFractions_[i] * a_coeff_vec(1, counter) *(dalphadT_vec_Curr_[i] / alpha_vec_Curr_[i]); + } + m_pp[k] = m_pp[k] * dalphadT_vec_Curr_[k] / alpha_vec_Curr_[k]; + } + + + for (size_t k = 0; k < m_kk; k++) { + coeff1 = m_b_current * (m_pp[k] + m_tmpV[k]) - daAlphadT * b_vec_Curr_[k]; + sbar[k] += GasConstant * log(GasConstant * TKelvin / (refP * mv)) + + GasConstant + + GasConstant * log(mv / vmb) + + GasConstant * b_vec_Curr_[k] / vmb + - coeff1* log(vpb2 / vmb2) / den1 + - b_vec_Curr_[k] * mv * daAlphadT / den2 / m_b_current; + } + pressureDerivatives(); + getPartialMolarVolumes(m_partialMolarVolumes.data()); + for (size_t k = 0; k < m_kk; k++) { + sbar[k] -= m_partialMolarVolumes[k] * dpdT_; + } +} + +void PengRobinsonMFTP::getPartialMolarIntEnergies(double* ubar) const +{ + getIntEnergy_RT(ubar); + scale(ubar, ubar+m_kk, ubar, RT()); +} + +void PengRobinsonMFTP::getPartialMolarCp(double* cpbar) const +{ + getCp_R(cpbar); + scale(cpbar, cpbar+m_kk, cpbar, GasConstant); +} + +void PengRobinsonMFTP::getPartialMolarVolumes(double* vbar) const +{ + for (size_t k = 0; k < m_kk; k++) { + m_pp[k] = 0.0; + for (size_t i = 0; i < m_kk; i++) { + size_t counter = k + m_kk*i; + m_pp[k] += moleFractions_[i] * aAlpha_vec_Curr_[counter]; + } + } + + double mv = molarVolume(); + double vmb = mv - m_b_current; + double vpb = mv + m_b_current; + double fac = mv * mv + 2 * mv * m_b_current - m_b_current * m_b_current; + double fac2 = fac * fac; + double RTkelvin = RT(); + + for (size_t k = 0; k < m_kk; k++) { + double num = (RTkelvin + RTkelvin * m_b_current/ vmb + RTkelvin * b_vec_Curr_[k] / vmb + + RTkelvin * m_b_current * b_vec_Curr_[k] /(vmb * vmb) + - 2 * mv * m_pp[k] / fac + + 2 * mv * vmb * m_aAlpha_current * b_vec_Curr_[k] / fac2 + ); + double denom = (pressure() + RTkelvin * m_b_current / (vmb * vmb) + + m_aAlpha_current/fac + - 2 * mv* vpb *m_aAlpha_current / fac2 + ); + vbar[k] = num / denom; + } +} + +double PengRobinsonMFTP::speciesCritTemperature(double a, double b) const +{ + double pc, tc, vc; + calcCriticalConditions(a, b, pc, tc, vc); + return tc; +} + +double PengRobinsonMFTP::critTemperature() const +{ + double pc, tc, vc; + calcCriticalConditions(m_a_current, m_b_current, pc, tc, vc); + return tc; +} + +double PengRobinsonMFTP::critPressure() const +{ + double pc, tc, vc; + calcCriticalConditions(m_a_current, m_b_current, pc, tc, vc); + return pc; +} + +double PengRobinsonMFTP::critVolume() const +{ + double pc, tc, vc; + calcCriticalConditions(m_a_current, m_b_current, pc, tc, vc); + return vc; +} + +double PengRobinsonMFTP::critCompressibility() const +{ + double pc, tc, vc; + calcCriticalConditions(m_a_current, m_b_current, pc, tc, vc); + return pc*vc/tc/GasConstant; +} + +double PengRobinsonMFTP::critDensity() const +{ + double pc, tc, vc; + calcCriticalConditions(m_a_current, m_b_current, pc, tc, vc); + double mmw = meanMolecularWeight(); + return mmw / vc; +} + +void PengRobinsonMFTP::setToEquilState(const double* mu_RT) +{ + double tmp, tmp2; + _updateReferenceStateThermo(); + getGibbs_RT_ref(m_tmpV.data()); + + // Within the method, we protect against inf results if the exponent is too + // high. + // + // If it is too low, we set the partial pressure to zero. This capability is + // needed by the elemental potential method. + double pres = 0.0; + double m_p0 = refPressure(); + for (size_t k = 0; k < m_kk; k++) { + tmp = -m_tmpV[k] + mu_RT[k]; + if (tmp < -600.) { + m_pp[k] = 0.0; + } else if (tmp > 500.0) { + tmp2 = tmp / 500.; + tmp2 *= tmp2; + m_pp[k] = m_p0 * exp(500.) * tmp2; + } else { + m_pp[k] = m_p0 * exp(tmp); + } + pres += m_pp[k]; + } + // set state + setState_PX(pres, &m_pp[0]); +} + +bool PengRobinsonMFTP::addSpecies(shared_ptr spec) +{ + bool added = MixtureFugacityTP::addSpecies(spec); + if (added) { + a_vec_Curr_.resize(m_kk * m_kk, 0.0); + b_vec_Curr_.push_back(0.0); + a_vec_Curr_.push_back(0.0); + aAlpha_vec_Curr_.resize(m_kk * m_kk, 0.0); + aAlpha_vec_Curr_.push_back(0.0); + kappa_vec_.push_back(0.0); + + alpha_vec_Curr_.push_back(0.0); + a_coeff_vec.resize(1, m_kk * m_kk, 0.0); + aAlpha_coeff_vec.resize(1, m_kk * m_kk, 0.0); + dalphadT_vec_Curr_.push_back(0.0); + d2alphadT2_.push_back(0.0); + + m_pp.push_back(0.0); + m_tmpV.push_back(0.0); + m_partialMolarVolumes.push_back(0.0); + dpdni_.push_back(0.0); + } + return added; +} + +vector PengRobinsonMFTP::getCoeff(const std::string& iName) +{ + vector_fp spCoeff{ NAN, NAN }; + + // Get number of species in the database + // open xml file critProperties.xml + XML_Node* doc = get_XML_File("critProperties.xml"); + size_t nDatabase = doc->nChildren(); + + // Loop through all species in the database and attempt to match supplied + // species to each. If present, calculate pureFluidParameters a_k and b_k + // based on crit properties T_c and P_c: + for (size_t isp = 0; isp < nDatabase; isp++) { + XML_Node& acNodeDoc = doc->child(isp); + std::string iNameLower = toLowerCopy(iName); + std::string dbName = toLowerCopy(acNodeDoc.attrib("name")); + + // Attempt to match provided species iName to current database species + // dbName: + if (iNameLower == dbName) { + // Read from database and calculate a and b coefficients + double vParams; + double T_crit, P_crit; + + if (acNodeDoc.hasChild("Tc")) { + vParams = 0.0; + XML_Node& xmlChildCoeff = acNodeDoc.child("Tc"); + if (xmlChildCoeff.hasAttrib("value")) { + std::string critTemp = xmlChildCoeff.attrib("value"); + vParams = strSItoDbl(critTemp); + } + if (vParams <= 0.0) { //Assuming that Pc and Tc are non zero. + throw CanteraError("PengRobinsonMFTP::getCoeff", + "Critical Temperature must be positive "); + } + T_crit = vParams; + } + if (acNodeDoc.hasChild("Pc")) { + vParams = 0.0; + XML_Node& xmlChildCoeff = acNodeDoc.child("Pc"); + if (xmlChildCoeff.hasAttrib("value")) { + std::string critPressure = xmlChildCoeff.attrib("value"); + vParams = strSItoDbl(critPressure); + } + if (vParams <= 0.0) { //Assuming that Pc and Tc are non zero. + throw CanteraError("PengRobinsonMFTP::getCoeff", + "Critical Pressure must be positive "); + } + P_crit = vParams; + } + + //Assuming no temperature dependence + spCoeff[0] = omega_a * (GasConstant* GasConstant) * (T_crit* T_crit) / P_crit; //coeff a + spCoeff[1] = omega_b * GasConstant * T_crit / P_crit; // coeff b + break; + } + } + return spCoeff; +} + +void PengRobinsonMFTP::initThermoXML(XML_Node& phaseNode, const std::string& id) +{ + if (phaseNode.hasChild("thermo")) { + XML_Node& thermoNode = phaseNode.child("thermo"); + std::string model = thermoNode["model"]; + if (model != "PengRobinson" && model != "PengRobinsonMFTP") { + throw CanteraError("PengRobinsonMFTP::initThermoXML", + "Unknown thermo model : " + model); + } + + // Go get all of the coefficients and factors in the + // activityCoefficients XML block + if (thermoNode.hasChild("activityCoefficients")) { + XML_Node& acNode = thermoNode.child("activityCoefficients"); + + // Count the number of species with parameters provided in the + // input file: + size_t nParams = 0; + + // Loop through the children and read out fluid parameters. Process + // all the pureFluidParameters, first: + for (size_t i = 0; i < acNode.nChildren(); i++) { + XML_Node& xmlACChild = acNode.child(i); + if (caseInsensitiveEquals(xmlACChild.name(), "purefluidparameters")) { + readXMLPureFluid(xmlACChild); + nParams += 1; + } + } + + // If any species exist which have undefined pureFluidParameters, + // search the database in 'critProperties.xml' to find critical + // temperature and pressure to calculate a and b. + + // Loop through all species in the CTI file + size_t iSpecies = 0; + + for (size_t i = 0; i < m_kk; i++) { + string iName = speciesName(i); + + // Get the index of the species + iSpecies = speciesIndex(iName); + + // Check if a and b are already populated (only the diagonal elements of a). + size_t counter = iSpecies + m_kk * iSpecies; + + // If not, then search the database: + if (isnan(a_coeff_vec(0, counter))) { + + vector coeffArray; + + // Search the database for the species name and calculate + // coefficients a and b, from critical properties: + // coeffArray[0] = a0, coeffArray[1] = b, coeffArray[2] = w; + coeffArray = getCoeff(iName); + + // Check if species was found in the database of critical properties, + // and assign the results + if (!isnan(coeffArray[0])) { + //Assuming no temperature dependence (i,e a1 = 0) + setSpeciesCoeffs(iName, coeffArray[0], 0.0, coeffArray[1]); + } + } + } + + // Loop back through the "activityCoefficients" children and process the + // crossFluidParameters in the XML tree: + for (size_t i = 0; i < acNode.nChildren(); i++) { + XML_Node& xmlACChild = acNode.child(i); + if (caseInsensitiveEquals(xmlACChild.name(), "crossfluidparameters")) { + readXMLCrossFluid(xmlACChild); + } + } + } + } + + MixtureFugacityTP::initThermoXML(phaseNode, id); +} + +void PengRobinsonMFTP::readXMLPureFluid(XML_Node& pureFluidParam) +{ + string xname = pureFluidParam.name(); + if (xname != "pureFluidParameters") { + throw CanteraError("PengRobinsonMFTP::readXMLPureFluid", + "Incorrect name for processing this routine: " + xname); + } + + double a0 = 0.0; + double a1 = 0.0; + double b = 0.0; + double w = 0.0; + for (size_t iChild = 0; iChild < pureFluidParam.nChildren(); iChild++) { + XML_Node& xmlChild = pureFluidParam.child(iChild); + string nodeName = toLowerCopy(xmlChild.name()); + + if (nodeName == "a_coeff") { + vector_fp vParams; + string iModel = toLowerCopy(xmlChild.attrib("model")); + getFloatArray(xmlChild, vParams, true, "Pascal-m6/kmol2", "a_coeff"); + + if (vParams.size() == 1) { + a0 = vParams[0]; + } else if (vParams.size() == 2) { + a0 = vParams[0]; + a1 = vParams[1]; + } else { + throw CanteraError("PengRobinsonMFTP::readXMLPureFluid", + "unknown model or incorrect number of parameters"); + } + } else if (nodeName == "b_coeff") { + b = getFloatCurrent(xmlChild, "toSI"); + } else if (nodeName == "acentric_factor") { + w = getFloatCurrent(xmlChild); + } + } + calculateAlpha(pureFluidParam.attrib("species"), a0, b, w); + setSpeciesCoeffs(pureFluidParam.attrib("species"), a0, b, w); +} + +void PengRobinsonMFTP::readXMLCrossFluid(XML_Node& CrossFluidParam) +{ + string xname = CrossFluidParam.name(); + if (xname != "crossFluidParameters") { + throw CanteraError("PengRobinsonMFTP::readXMLCrossFluid", + "Incorrect name for processing this routine: " + xname); + } + + string iName = CrossFluidParam.attrib("species1"); + string jName = CrossFluidParam.attrib("species2"); + + size_t num = CrossFluidParam.nChildren(); + for (size_t iChild = 0; iChild < num; iChild++) { + XML_Node& xmlChild = CrossFluidParam.child(iChild); + string nodeName = toLowerCopy(xmlChild.name()); + + if (nodeName == "a_coeff") { + vector_fp vParams; + getFloatArray(xmlChild, vParams, true, "Pascal-m6/kmol2", "a_coeff"); + string iModel = toLowerCopy(xmlChild.attrib("model")); + if (iModel == "constant" && vParams.size() == 1) { + setBinaryCoeffs(iName, jName, vParams[0], 0.0); + } else if (iModel == "linear_a") { + setBinaryCoeffs(iName, jName, vParams[0], vParams[1]); + } else { + throw CanteraError("PengRobinsonMFTP::readXMLCrossFluid", + "unknown model ({}) or wrong number of parameters ({})", + iModel, vParams.size()); + } + } + } +} + +void PengRobinsonMFTP::setParametersFromXML(const XML_Node& thermoNode) +{ + MixtureFugacityTP::setParametersFromXML(thermoNode); +} + +double PengRobinsonMFTP::sresid() const +{ + double molarV = molarVolume(); + double hh = m_b_current / molarV; + double zz = z(); + double alpha_1 = daAlpha_dT(); + double T = temperature(); + double vpb = molarV + (1.0 + M_SQRT2) *m_b_current; + double vmb = molarV + (1.0 - M_SQRT2) *m_b_current; + double fac = alpha_1 / (2.0 * M_SQRT2 * m_b_current); + double sresid_mol_R = log(zz*(1.0 - hh)) + fac * log(vpb / vmb) / GasConstant; + return GasConstant * sresid_mol_R; +} + +double PengRobinsonMFTP::hresid() const +{ + double molarV = molarVolume(); + double zz = z(); + double aAlpha_1 = daAlpha_dT(); + double T = temperature(); + double vpb = molarV + (1 + M_SQRT2) *m_b_current; + double vmb = molarV + (1 - M_SQRT2) *m_b_current; + double fac = 1 / (2.0 * M_SQRT2 * m_b_current); + return GasConstant * T * (zz - 1.0) + fac * log(vpb / vmb) *(T * aAlpha_1 - m_aAlpha_current); +} + +double PengRobinsonMFTP::liquidVolEst(double TKelvin, double& presGuess) const +{ + double v = m_b_current * 1.1; + double atmp; + double btmp; + double aAlphatmp; + calculateAB(TKelvin, atmp, btmp, aAlphatmp); + double pres = std::max(psatEst(TKelvin), presGuess); + double Vroot[3]; + bool foundLiq = false; + int m = 0; + while (m < 100 && !foundLiq) { + int nsol = NicholsSolve(TKelvin, 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 PengRobinsonMFTP::densityCalc(double TKelvin, double presPa, int phaseRequested, double rhoGuess) +{ + // It's necessary to set the temperature so that m_aAlpha_current is set correctly. + setTemperature(TKelvin); + double tcrit = critTemperature(); + double mmw = meanMolecularWeight(); + if (rhoGuess == -1.0) { + if (phaseRequested != FLUID_GAS) { + if (TKelvin > tcrit) { + rhoGuess = presPa * mmw / (GasConstant * TKelvin); + } else { + if (phaseRequested == FLUID_UNDEFINED || phaseRequested == FLUID_SUPERCRIT) { + rhoGuess = presPa * mmw / (GasConstant * TKelvin); + } else if (phaseRequested >= FLUID_LIQUID_0) { + double lqvol = liquidVolEst(TKelvin, presPa); + rhoGuess = mmw / lqvol; + } + } + } else { + // Assume the Gas phase initial guess, if nothing is specified to the routine + rhoGuess = presPa * mmw / (GasConstant * TKelvin); + } + } + + double volGuess = mmw / rhoGuess; + NSolns_ = NicholsSolve(TKelvin, presPa, m_a_current, m_b_current, m_aAlpha_current, Vroot_); + + double molarVolLast = Vroot_[0]; + if (NSolns_ >= 2) { + if (phaseRequested >= FLUID_LIQUID_0) { + molarVolLast = Vroot_[0]; + } else if (phaseRequested == FLUID_GAS || phaseRequested == FLUID_SUPERCRIT) { + molarVolLast = Vroot_[2]; + } else { + if (volGuess > Vroot_[1]) { + molarVolLast = Vroot_[2]; + } else { + molarVolLast = Vroot_[0]; + } + } + } else if (NSolns_ == 1) { + if (phaseRequested == FLUID_GAS || phaseRequested == FLUID_SUPERCRIT || phaseRequested == FLUID_UNDEFINED) { + molarVolLast = Vroot_[0]; + } else { + return -2.0; + } + } else if (NSolns_ == -1) { + if (phaseRequested >= FLUID_LIQUID_0 || phaseRequested == FLUID_UNDEFINED || phaseRequested == FLUID_SUPERCRIT) { + molarVolLast = Vroot_[0]; + } else if (TKelvin > tcrit) { + molarVolLast = Vroot_[0]; + } else { + return -2.0; + } + } else { + molarVolLast = Vroot_[0]; + return -1.0; + } + return mmw / molarVolLast; +} + +double PengRobinsonMFTP::densSpinodalLiquid() const +{ + double Vroot[3]; + double T = temperature(); + int nsol = NicholsSolve(T, pressure(), m_a_current, m_b_current, m_aAlpha_current, 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 PengRobinsonMFTP::densSpinodalGas() const +{ + double Vroot[3]; + double T = temperature(); + int nsol = NicholsSolve(T, pressure(), m_a_current, m_b_current, m_aAlpha_current, 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 PengRobinsonMFTP::pressureCalc(double TKelvin, double molarVol) const +{ + double den = molarVol * molarVol + 2 * molarVol * m_b_current - m_b_current * m_b_current; + double pres = GasConstant * TKelvin / (molarVol - m_b_current) - m_aAlpha_current / den; + return pres; +} + +double PengRobinsonMFTP::dpdVCalc(double TKelvin, double molarVol, double& presCalc) const +{ + double den = molarVol * molarVol + 2 * molarVol * m_b_current - m_b_current * m_b_current; + presCalc = GasConstant * TKelvin / (molarVol - m_b_current) - m_aAlpha_current/ den; + + double vpb = molarVol + m_b_current; + double vmb = molarVol - m_b_current; + double dpdv = -GasConstant * TKelvin / (vmb * vmb) + 2 *m_aAlpha_current * vpb / (den*den); + return dpdv; +} + +void PengRobinsonMFTP::pressureDerivatives() const +{ + double TKelvin = temperature(); + double mv = molarVolume(); + double pres; + + dpdV_ = dpdVCalc(TKelvin, mv, pres); + double vmb = mv - m_b_current; + double den = mv * mv + 2 * mv * m_b_current - m_b_current * m_b_current; + dpdT_ = (GasConstant / vmb - daAlpha_dT() / den); +} + +void PengRobinsonMFTP::updateMixingExpressions() +{ + updateAB(); +} + +void PengRobinsonMFTP::updateAB() +{ + double temp = temperature(); + //Update aAlpha_i + double sqt_alpha; + double criTemp = critTemperature(); + double sqt_T_reduced = sqrt(temp / criTemp); + + // Update indiviual alpha + for (size_t j = 0; j < m_kk; j++) { + sqt_alpha = 1 + kappa_vec_[j] * (1 - sqt_T_reduced); + alpha_vec_Curr_[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++) { + size_t counter = i * m_kk + j; + a_vec_Curr_[counter] = a_coeff_vec(0, counter); + aAlpha_vec_Curr_[counter] = sqrt(alpha_vec_Curr_[i] * alpha_vec_Curr_[j]) * a_coeff_vec(0, counter); + } + } + + m_b_current = 0.0; + m_a_current = 0.0; + m_aAlpha_current = 0.0; + + for (size_t i = 0; i < m_kk; i++) { + m_b_current += moleFractions_[i] * b_vec_Curr_[i]; + for (size_t j = 0; j < m_kk; j++) { + m_a_current += a_vec_Curr_[i * m_kk + j] * moleFractions_[i] * moleFractions_[j]; + m_aAlpha_current += aAlpha_vec_Curr_[i * m_kk + j] * moleFractions_[i] * moleFractions_[j]; + } + } +} + +void PengRobinsonMFTP::calculateAB(double temp, double& aCalc, double& bCalc, double& aAlphaCalc) const +{ + bCalc = 0.0; + aCalc = 0.0; + aAlphaCalc = 0.0; + for (size_t i = 0; i < m_kk; i++) { + bCalc += moleFractions_[i] * b_vec_Curr_[i]; + for (size_t j = 0; j < m_kk; j++) { + size_t counter = i * m_kk + j; + double a_vec_Curr = a_coeff_vec(0, counter); + aCalc += a_vec_Curr * moleFractions_[i] * moleFractions_[j]; + aAlphaCalc += aAlpha_vec_Curr_[counter] * moleFractions_[i] * moleFractions_[j]; + } + } +} + +double PengRobinsonMFTP::daAlpha_dT() const +{ + double daAlphadT = 0.0, temp, k, Tc = 0.0, sqtTr = 0.0; + double coeff1, coeff2; + for (size_t i = 0; i < m_kk; i++) { + size_t counter = i + m_kk * i; + // Calculate first derivative of alpha for individual species + Tc = speciesCritTemperature(a_vec_Curr_[counter], b_vec_Curr_[i]); + sqtTr = sqrt(temperature() / Tc); //we need species critical temperature + coeff1 = 1 / (Tc*sqtTr); + coeff2 = sqtTr - 1; + k = kappa_vec_[i]; + dalphadT_vec_Curr_[i] = coeff1 *(k* k*coeff2 - k); + } + //Calculate mixture derivative + for (size_t i = 0; i < m_kk; i++) { + size_t counter1 = i + m_kk * i; + for (size_t j = 0; j < m_kk; j++) { + size_t counter2 = j * m_kk + j; + temp = 0.5 * sqrt((a_vec_Curr_[counter1] * a_vec_Curr_[counter2]) / (alpha_vec_Curr_[i] * alpha_vec_Curr_[j])); + daAlphadT += moleFractions_[i] * moleFractions_[j] * temp + * (dalphadT_vec_Curr_[j] * alpha_vec_Curr_[i] + dalphadT_vec_Curr_[i] * alpha_vec_Curr_[j]); + } + } + return daAlphadT; +} + +double PengRobinsonMFTP::d2aAlpha_dT2() const +{ + double daAlphadT = 0.0, temp, fac1, fac2, alphaij, alphai, alphaj, d2aAlphadT2 = 0.0, num; + double k; + double sqt_Tr = sqrt(temperature() / critTemperature()); //we need species critical temperature + double coeff1 = 1 / (critTemperature()*critTemperature()*sqt_Tr); + double coeff2 = sqt_Tr - 1; + for (size_t i = 0; i < m_kk; i++) { + // Calculate first and second derivatives of alpha for individual species + size_t counter = i + m_kk * i; + k = kappa_vec_[i]; + dalphadT_vec_Curr_[i] = coeff1 *(k* k*coeff2 - k); + d2alphadT2_[i] = (k*k + k) * coeff1 / (2 * sqt_Tr*sqt_Tr); + } + + //Calculate mixture derivative + for (size_t i = 0; i < m_kk; i++) { + size_t counter1 = i + m_kk * i; + alphai = alpha_vec_Curr_[i]; + for (size_t j = 0; j < m_kk; j++) { + size_t counter2 = j + m_kk * j; + alphaj = alpha_vec_Curr_[j]; + alphaij = alphai * alphaj; + temp = 0.5 * sqrt((a_vec_Curr_[counter1] * a_vec_Curr_[counter2]) / (alphaij)); + num = (dalphadT_vec_Curr_[j] * alphai + dalphadT_vec_Curr_[i] * alphaj); + fac1 = -(0.5 / alphaij)*num*num; + fac2 = alphaj * d2alphadT2_[counter1] + alphai *d2alphadT2_[counter2] + 2 * dalphadT_vec_Curr_[i] * dalphadT_vec_Curr_[j]; + d2aAlphadT2 += moleFractions_[i] * moleFractions_[j] * temp *(fac1 + fac2); + } + } + return d2aAlphadT2; +} + +void PengRobinsonMFTP::calcCriticalConditions(double a, double b, + double& pc, double& tc, double& vc) const +{ + if (b <= 0.0) { + tc = 1000000.; + pc = 1.0E13; + vc = omega_vc * GasConstant * tc / pc; + return; + } + if (a <= 0.0) { + tc = 0.0; + pc = 0.0; + vc = 2.0 * b; + return; + } + tc = a * omega_b / (b * omega_a * GasConstant); + pc = omega_b * GasConstant * tc / b; + vc = omega_vc * GasConstant * tc / pc; +} + +int PengRobinsonMFTP::NicholsSolve(double TKelvin, double pres, double a, double b, double aAlpha, + double Vroot[3]) const +{ + double tmp; + fill_n(Vroot, 3, 0.0); + if (TKelvin <= 0.0) { + throw CanteraError("PengRobinsonMFTP::NicholsSolve()", "negative temperature T = {}", TKelvin); + } + + // Derive the coefficients of the cubic polynomial (in terms of molar volume v) to solve. + double bsqr = b * b; + double RT_p = GasConstant * TKelvin / pres; + double aAlpha_p = aAlpha / pres; + double an = 1.0; + double bn = (b - RT_p); + double cn = -(2 * RT_p * b - aAlpha_p + 3 * bsqr); + double dn = (bsqr * RT_p + bsqr * b - aAlpha_p * b); + + double tc = a * omega_b / (b * omega_a * GasConstant); + double pc = omega_b * GasConstant * tc / b; + double vc = omega_vc * GasConstant * tc / pc; + + // Derive the center of the cubic, x_N + double xN = - bn /(3 * an); + + // Derive the value of delta**2. This is a key quantity that determines the number of turning points + double delta2 = (bn * bn - 3 * an * cn) / (9 * an * an); + double delta = 0.0; + + // Calculate a couple of ratios + // Cubic equation in z : z^3 - (1-B) z^2 + (A -2B -3B^2)z - (AB- B^2- B^3) = 0 + double ratio1 = 3.0 * an * cn / (bn * bn); + double ratio2 = pres * b / (GasConstant * TKelvin); // B + if (fabs(ratio1) < 1.0E-7) { + double ratio3 = aAlpha / (GasConstant * TKelvin) * pres / (GasConstant * TKelvin); // A + if (fabs(ratio2) < 1.0E-5 && fabs(ratio3) < 1.0E-5) { + // A and B terms in cubic equation for z are almost zero, then z is near to 1 + double zz = 1.0; + for (int i = 0; i < 10; i++) { + double znew = zz / (zz - ratio2) - ratio3 / (zz + ratio1); + double deltaz = znew - zz; + zz = znew; + if (fabs(deltaz) < 1.0E-14) { + break; + } + } + double v = zz * GasConstant * TKelvin / pres; + Vroot[0] = v; + return 1; + } + } + + int nSolnValues; // Represents number of solutions to the cubic equation + double h2 = 4. * an * an * delta2 * delta2 * delta2; // h^2 + if (delta2 > 0.0) { + delta = sqrt(delta2); + } + + double h = 2.0 * an * delta * delta2; + double yN = 2.0 * bn * bn * bn / (27.0 * an * an) - bn * cn / (3.0 * an) + dn; // y_N term + double disc = yN * yN - h2; // discriminant + + //check if y = h + if (fabs(fabs(h) - fabs(yN)) < 1.0E-10) { + if (disc > 1e-10) { + throw CanteraError("PengRobinsonMFTP::NicholsSolve()", "value of yN and h are too high, unrealistic roots may be obtained"); + } + disc = 0.0; + } + + if (disc < -1e-14) { + // disc<0 then we have three distinct roots. + nSolnValues = 3; + } else if (fabs(disc) < 1e-14) { + // disc=0 then we have two distinct roots (third one is repeated root) + nSolnValues = 2; + // We are here as p goes to zero. + } else if (disc > 1e-14) { + // disc> 0 then we have one real root. + nSolnValues = 1; + } + + // One real root -> have to determine whether gas or liquid is the root + if (disc > 0.0) { + double tmpD = sqrt(disc); + double tmp1 = (- yN + tmpD) / (2.0 * an); + double sgn1 = 1.0; + if (tmp1 < 0.0) { + sgn1 = -1.0; + tmp1 = -tmp1; + } + double tmp2 = (- yN - tmpD) / (2.0 * an); + double sgn2 = 1.0; + if (tmp2 < 0.0) { + sgn2 = -1.0; + tmp2 = -tmp2; + } + double p1 = pow(tmp1, 1./3.); + double p2 = pow(tmp2, 1./3.); + double alpha = xN + sgn1 * p1 + sgn2 * p2; + Vroot[0] = alpha; + Vroot[1] = 0.0; + Vroot[2] = 0.0; + } else if (disc < 0.0) { + // Three real roots alpha, beta, gamma are obtained. + double val = acos(-yN / h); + double theta = val / 3.0; + double twoThirdPi = 2. * Pi / 3.; + double alpha = xN + 2. * delta * cos(theta); + double beta = xN + 2. * delta * cos(theta + twoThirdPi); + double gamma = xN + 2. * delta * cos(theta + 2.0 * twoThirdPi); + Vroot[0] = beta; + Vroot[1] = gamma; + Vroot[2] = alpha; + + for (int i = 0; i < 3; i++) { + tmp = an * Vroot[i] * Vroot[i] * Vroot[i] + bn * Vroot[i] * Vroot[i] + cn * Vroot[i] + dn; + if (fabs(tmp) > 1.0E-4) { + for (int j = 0; j < 3; j++) { + if (j != i && fabs(Vroot[i] - Vroot[j]) < 1.0E-4 * (fabs(Vroot[i]) + fabs(Vroot[j]))) { + writelog("PengRobinsonMFTP::NicholsSolve(T = {}, p = {}):" + " WARNING roots have merged: {}, {}\n", + TKelvin, pres, Vroot[i], Vroot[j]); + } + } + } + } + } else if (disc == 0.0) { + //Three equal roots are obtained, i.e. alpha = beta = gamma + if (yN < 1e-18 && h < 1e-18) { + // yN = 0.0 and h = 0 i.e. disc = 0 + Vroot[0] = xN; + Vroot[1] = xN; + Vroot[2] = xN; + } else { + // h and yN need to figure out whether delta^3 is positive or negative + if (yN > 0.0) { + tmp = pow(yN/(2*an), 1./3.); + // In this case, tmp and delta must be equal. + if (fabs(tmp - delta) > 1.0E-9) { + throw CanteraError("PengRobinsonMFTP::NicholsSolve()", "Inconsistancy in cubic solver : solver is bad conditioned."); + } + Vroot[1] = xN + delta; + Vroot[0] = xN - 2.0*delta; // liquid phase root + } else { + tmp = pow(yN/(2*an), 1./3.); + // In this case, tmp and delta must be equal. + if (fabs(tmp - delta) > 1.0E-9) { + throw CanteraError("PengRobinsonMFTP::NicholsSolve()", "Inconsistancy in cubic solver : solver is bad conditioned."); + } + delta = -delta; + Vroot[0] = xN + delta; + Vroot[1] = xN - 2.0*delta; // gas phase root + } + } + } + + // Find an accurate root, since there might be a heavy amount of roundoff error due to bad conditioning in this solver. + double res, dresdV = 0.0; + for (int i = 0; i < nSolnValues; i++) { + for (int n = 0; n < 20; n++) { + res = an * Vroot[i] * Vroot[i] * Vroot[i] + bn * Vroot[i] * Vroot[i] + cn * Vroot[i] + dn; + if (fabs(res) < 1.0E-14) { // accurate root is obtained + break; + } + dresdV = 3.0 * an * Vroot[i] * Vroot[i] + 2.0 * bn * Vroot[i] + cn; // derivative of the residual + double del = - res / dresdV; + Vroot[i] += del; + if (fabs(del) / (fabs(Vroot[i]) + fabs(del)) < 1.0E-14) { + break; + } + double res2 = an * Vroot[i] * Vroot[i] * Vroot[i] + bn * Vroot[i] * Vroot[i] + cn * Vroot[i] + dn; + if (fabs(res2) < fabs(res)) { + continue; + } else { + Vroot[i] -= del; // Go back to previous value of Vroot. + Vroot[i] += 0.1 * del; // under-relax by 0.1 + } + } + if ((fabs(res) > 1.0E-14) && (fabs(res) > 1.0E-14 * fabs(dresdV) * fabs(Vroot[i]))) { + writelog("PengRobinsonMFTP::NicholsSolve(T = {}, p = {}): " + "WARNING root didn't converge V = {}", TKelvin, pres, Vroot[i]); + writelogendl(); + } + } + + if (nSolnValues == 1) { + if (TKelvin > tc) { + if (Vroot[0] < vc) { + // Liquid phase root + nSolnValues = -1; + } + } else { + if (Vroot[0] < xN) { + nSolnValues = -1; + } + } + } else { + if (nSolnValues == 2 && delta > 1e-14) { + nSolnValues = -2; + } + } + return nSolnValues; +} + +} diff --git a/src/thermo/ThermoFactory.cpp b/src/thermo/ThermoFactory.cpp index 0f634757997..c899bf5f04d 100644 --- a/src/thermo/ThermoFactory.cpp +++ b/src/thermo/ThermoFactory.cpp @@ -23,6 +23,7 @@ #include "cantera/thermo/IonsFromNeutralVPSSTP.h" #include "cantera/thermo/PureFluidPhase.h" #include "cantera/thermo/RedlichKwongMFTP.h" +#include "cantera/thermo/PengRobinsonMFTP.h" #include "cantera/thermo/ConstDensityThermo.h" #include "cantera/thermo/SurfPhase.h" #include "cantera/thermo/EdgePhase.h" @@ -88,6 +89,9 @@ ThermoFactory::ThermoFactory() reg("RedlichKwong", []() { return new RedlichKwongMFTP(); }); m_synonyms["RedlichKwongMFTP"] = "RedlichKwong"; m_synonyms["Redlich-Kwong"] = "RedlichKwong"; + reg("PengRobinson", []() { return new PengRobinsonMFTP(); }); + m_synonyms["PengRobinsonMFTP"] = "PengRobinson"; + m_synonyms["Peng-Robinson"] = "PengRobinson"; reg("MaskellSolidSolnPhase", []() { return new MaskellSolidSolnPhase(); }); m_synonyms["Maskell-solid-solution"] = "MaskellSolidSolnPhase"; reg("PureLiquidWater", []() { return new WaterSSTP(); }); diff --git a/test/data/co2_PR_example.cti b/test/data/co2_PR_example.cti new file mode 100644 index 00000000000..8692fef113a --- /dev/null +++ b/test/data/co2_PR_example.cti @@ -0,0 +1,168 @@ +# Transport data from file ../transport/gri30_tran.dat. + +units(length ="cm", time ="s", quantity ="mol", act_energy ="cal/mol") + + +PengRobinsonMFTP(name ="carbondioxide", + elements ="C O H N", + species ="""CO2 H2O H2 CO CH4 O2 N2""", + activity_coefficients = (pureFluidParameters(species="CO2", a_coeff = [3.958134E+11, 0], b_coeff = 26.6275, acentric_factor = 0.228), + pureFluidParameters(species="H2O", a_coeff = [5.998873E+11, 0], b_coeff = 18.9714, acentric_factor = 0.344), + pureFluidParameters(species="H2", a_coeff = [2.668423E+10, 0], b_coeff = 16.5478, acentric_factor = -0.22), + pureFluidParameters(species="CO", a_coeff = [1.607164E+11, 0], b_coeff = 24.6549, acentric_factor = 0.049), + pureFluidParameters(species="CH4", a_coeff = [2.496344E+11, 0], b_coeff = 26.8028, acentric_factor = 0.01), + pureFluidParameters(species="O2", a_coeff = [1.497732E+11, 0], b_coeff = 19.8281, acentric_factor = 0.022), + pureFluidParameters(species="N2", a_coeff = [1.485031E+11, 0], b_coeff = 28.0810, acentric_factor = 0.04)), + transport ="Multi", + reactions ="all", + initial_state = state(temperature = 300.0, + pressure = OneAtm, + mole_fractions = 'CO2:0.99, H2:0.01')) + +#------------------------------------------------------------------------------- +# Species data +#------------------------------------------------------------------------------- + +species(name ="H2", + atoms ="H:2", + thermo = ( + NASA([200.00, 1000.00], [2.344331120E+00, 7.980520750E-03, + -1.947815100E-05, 2.015720940E-08, -7.376117610E-12, + -9.179351730E+02, 6.830102380E-01]), + NASA([1000.00, 3500.00], [3.337279200E+00, -4.940247310E-05, + 4.994567780E-07, -1.795663940E-10, 2.002553760E-14, + -9.501589220E+02, -3.205023310E+00]) + ), + transport = gas_transport( + geom ="linear", + diam = 2.92, + well_depth = 38.00, + polar = 0.79, + rot_relax = 280.00), + note ="TPIS78" + ) + +species(name ="CO", + atoms ="C:1 O:1", + thermo = ( + NASA([200.00, 1000.00], [3.579533470E+00, -6.103536800E-04, + 1.016814330E-06, 9.070058840E-10, -9.044244990E-13, + -1.434408600E+04, 3.508409280E+00]), + NASA([1000.00, 3500.00], [2.715185610E+00, 2.062527430E-03, + -9.988257710E-07, 2.300530080E-10, -2.036477160E-14, + -1.415187240E+04, 7.818687720E+00]) + ), + transport = gas_transport( + geom ="linear", + diam = 3.65, + well_depth = 98.10, + polar = 1.95, + rot_relax = 1.80), + note ="TPIS79" + ) + + +species(name ="N2", + atoms ="N:2", + thermo = ( + NASA([300.00, 1000.00], [3.298677000E+00, 1.408240400E-03, + -3.963222000E-06, 5.641515000E-09, -2.444854000E-12, + -1.020899900E+03, 3.950372000E+00]), + NASA([1000.00, 5000.00], [2.926640000E+00, 1.487976800E-03, + -5.684760000E-07, 1.009703800E-10, -6.753351000E-15, + -9.227977000E+02, 5.980528000E+00]) + ), + transport = gas_transport( + geom ="linear", + diam = 3.62, + well_depth = 97.53, + polar = 1.76, + rot_relax = 4.00), + note ="121286" + ) + +species(name ="O2", + atoms ="O:2", + thermo = ( + NASA([200.00, 1000.00], [3.782456360E+00, -2.996734160E-03, + 9.847302010E-06, -9.681295090E-09, 3.243728370E-12, + -1.063943560E+03, 3.657675730E+00]), + NASA([1000.00, 3500.00], [3.282537840E+00, 1.483087540E-03, + -7.579666690E-07, 2.094705550E-10, -2.167177940E-14, + -1.088457720E+03, 5.453231290E+00]) + ), + transport = gas_transport( + geom ="linear", + diam = 3.46, + well_depth = 107.40, + polar = 1.60, + rot_relax = 3.80), + note ="TPIS89" + ) + + +species(name ="CO2", + atoms ="C:1 O:2", + thermo = ( + NASA([200.00, 1000.00], [2.356773520E+00, 8.984596770E-03, + -7.123562690E-06, 2.459190220E-09, -1.436995480E-13, + -4.837196970E+04, 9.901052220E+00]), + NASA([1000.00, 3500.00], [3.857460290E+00, 4.414370260E-03, + -2.214814040E-06, 5.234901880E-10, -4.720841640E-14, + -4.875916600E+04, 2.271638060E+00]) + ), + transport = gas_transport( + geom ="linear", + diam = 3.76, + well_depth = 244.00, + polar = 2.65, + rot_relax = 2.10), + note ="L 7/88" + ) + +species(name ="CH4", + atoms ="C:1 H:4", + thermo = ( + NASA([200.00, 1000.00], [5.149876130E+00, -1.367097880E-02, + 4.918005990E-05, -4.847430260E-08, 1.666939560E-11, + -1.024664760E+04, -4.641303760E+00]), + NASA([1000.00, 3500.00], [7.485149500E-02, 1.339094670E-02, + -5.732858090E-06, 1.222925350E-09, -1.018152300E-13, + -9.468344590E+03, 1.843731800E+01]) + ), + transport = gas_transport( + geom ="nonlinear", + diam = 3.75, + well_depth = 141.40, + polar = 2.60, + rot_relax = 13.00), + note ="L 8/88" + ) + + +species(name ="H2O", + atoms ="H:2 O:1", + thermo = ( + NASA([200.00, 1000.00], [4.198640560E+00, -2.036434100E-03, + 6.520402110E-06, -5.487970620E-09, 1.771978170E-12, + -3.029372670E+04, -8.490322080E-01]), + NASA([1000.00, 3500.00], [3.033992490E+00, 2.176918040E-03, + -1.640725180E-07, -9.704198700E-11, 1.682009920E-14, + -3.000429710E+04, 4.966770100E+00]) + ), + transport = gas_transport( + geom ="nonlinear", + diam = 2.60, + well_depth = 572.40, + dipole = 1.85, + rot_relax = 4.00), + note ="L 8/89" + ) + +#——————————————————————————————————————————— +# Reaction data +#——————————————————————————————————————————— + +# Reaction 1 +reaction("CO2 + H2 <=> CO + H2O", [1.2E+3, 0, 0]) + diff --git a/test/thermo/PengRobinsonMFTP_Test.cpp b/test/thermo/PengRobinsonMFTP_Test.cpp new file mode 100644 index 00000000000..9d8043cb5dc --- /dev/null +++ b/test/thermo/PengRobinsonMFTP_Test.cpp @@ -0,0 +1,205 @@ +#include "gtest/gtest.h" +#include "cantera/thermo/PengRobinsonMFTP.h" +#include "cantera/thermo/ThermoFactory.h" + + +namespace Cantera +{ + +class PengRobinsonMFTP_Test : public testing::Test +{ +public: + PengRobinsonMFTP_Test() { + test_phase.reset(newPhase("../data/co2_PR_example.cti")); + } + + //vary the composition of a co2-h2 mixture: + void set_r(const double r) { + vector_fp moleFracs(7); + moleFracs[0] = r; + moleFracs[2] = 1-r; + test_phase->setMoleFractions(&moleFracs[0]); + } + + std::unique_ptr test_phase; +}; + +TEST_F(PengRobinsonMFTP_Test, construct_from_cti) +{ + PengRobinsonMFTP* peng_robinson_phase = dynamic_cast(test_phase.get()); + EXPECT_TRUE(peng_robinson_phase != NULL); +} + +TEST_F(PengRobinsonMFTP_Test, chem_potentials) +{ + test_phase->setState_TP(298.15, 101325.); + /* Chemical potential should increase with increasing co2 mole fraction: + * mu = mu_0 + RT ln(gamma_k*X_k). + * where gamma_k is the activity coefficient. Run regression test against values + * calculated using the model. + */ + const double expected_result[9] = { + -4.5736182681761962e+008, + -4.5733771904416579e+008, + -4.5732943831449223e+008, + -4.5732206687414169e+008, + -4.5731546826955432e+008, + -4.5730953161186475e+008, + -4.5730416590547645e+008, + -4.5729929581635743e+008, + -4.5729485847173005e+008 + }; + + double xmin = 0.6; + double xmax = 0.9; + int numSteps = 9; + double dx = (xmax-xmin)/(numSteps-1); + vector_fp chemPotentials(7); + for(int i=0; i < numSteps; ++i) + { + set_r(xmin + i*dx); + test_phase->getChemPotentials(&chemPotentials[0]); + EXPECT_NEAR(expected_result[i], chemPotentials[0], 1.e-6); + } +} + +TEST_F(PengRobinsonMFTP_Test, activityCoeffs) +{ + test_phase->setState_TP(298., 1.); + + // Test that mu0 + RT log(activityCoeff * MoleFrac) == mu + const double RT = GasConstant * 298.; + vector_fp mu0(7); + vector_fp activityCoeffs(7); + vector_fp chemPotentials(7); + double xmin = 0.6; + double xmax = 0.9; + int numSteps = 9; + double dx = (xmax-xmin)/(numSteps-1); + + for(int i=0; i < numSteps; ++i) + { + const double r = xmin + i*dx; + set_r(r); + test_phase->getChemPotentials(&chemPotentials[0]); + test_phase->getActivityCoefficients(&activityCoeffs[0]); + test_phase->getStandardChemPotentials(&mu0[0]); + EXPECT_NEAR(chemPotentials[0], mu0[0] + RT*std::log(activityCoeffs[0] * r), 1.e-6); + EXPECT_NEAR(chemPotentials[2], mu0[2] + RT*std::log(activityCoeffs[2] * (1-r)), 1.e-6); + } +} + +TEST_F(PengRobinsonMFTP_Test, standardConcentrations) +{ + EXPECT_DOUBLE_EQ(test_phase->pressure()/(test_phase->temperature()*GasConstant), test_phase->standardConcentration(0)); + EXPECT_DOUBLE_EQ(test_phase->pressure()/(test_phase->temperature()*GasConstant), test_phase->standardConcentration(1)); +} + +TEST_F(PengRobinsonMFTP_Test, activityConcentrations) +{ + // Check to make sure activityConcentration_i == standardConcentration_i * gamma_i * X_i + vector_fp standardConcs(7); + vector_fp activityCoeffs(7); + vector_fp activityConcentrations(7); + double xmin = 0.6; + double xmax = 0.9; + int numSteps = 9; + double dx = (xmax-xmin)/(numSteps-1); + + for(int i=0; i < numSteps; ++i) + { + const double r = xmin + i*dx; + set_r(r); + test_phase->getActivityCoefficients(&activityCoeffs[0]); + standardConcs[0] = test_phase->standardConcentration(0); + standardConcs[2] = test_phase->standardConcentration(2); + test_phase->getActivityConcentrations(&activityConcentrations[0]); + + EXPECT_NEAR(standardConcs[0] * r * activityCoeffs[0], activityConcentrations[0], 1.e-6); + EXPECT_NEAR(standardConcs[2] * (1-r) * activityCoeffs[2], activityConcentrations[2], 1.e-6); + } +} + +TEST_F(PengRobinsonMFTP_Test, setTP) +{ + // Check to make sure that the phase diagram is accurately reproduced for a few select isobars + + // All sub-cooled liquid: + const double p1[6] = { + 1.7474528924963985e+002, + 1.6800540828415956e+002, + 1.62278413743154e+002, + 1.5728963799103039e+002, + 1.5286573762819748e+002, + 1.4888956030449546e+002 + }; + // Phase change between temperatures 4 & 5: + const double p2[6] = { + 7.5565889855724288e+002, + 7.2577747673480337e+002, + 6.913183942651284e+002, + 6.494661249672663e+002, + 5.9240469307757724e+002, + 3.645826047440932e+002 + }; + // Supercritical; no discontinuity in rho values: + const double p3[6] = { + 8.047430802847415e+002, + 7.8291565113595595e+002, + 7.5958477920749681e+002, + 7.3445460137134626e+002, + 7.0712433093853724e+002, + 6.77034438769492e+002 + }; + + for(int i=0; i<6; ++i) + { + const double temp = 294 + i*2; + set_r(0.999); + test_phase->setState_TP(temp, 5542027.5); + EXPECT_NEAR(test_phase->density(),p1[i],1.e-8); + + test_phase->setState_TP(temp, 7389370.); + EXPECT_NEAR(test_phase->density(),p2[i],1.e-8); + + test_phase->setState_TP(temp, 9236712.5); + EXPECT_NEAR(test_phase->density(),p3[i],1.e-8); + } +} + +TEST_F(PengRobinsonMFTP_Test, getPressure) +{ + // Check to make sure that the P-R equation is accurately reproduced for a few selected values + + /* This test uses CO2 as the only species. + * Values of a_coeff, b_coeff are calculated based on the the critical temperature and pressure values of CO2 as follows: + * a_coeff = 0.457235(RT_crit)^2/p_crit + * b_coeff = 0.077796(RT_crit)/p_crit + * The temperature dependent parameter in P-R EoS is calculated as + * \alpha = [1 + \kappa(1 - sqrt{T/T_crit}]^2 + * kappa is a function calulated based on the accentric factor. + */ + + double a_coeff = 3.958134E+5; + double b_coeff = 26.6275/1000; + double acc_factor = 0.228; + double pres_theoretical, kappa, alpha, mv; + const double rho = 1.0737; + const double Tcrit = test_phase->critTemperature(); + + //Calculate kappa value + kappa = 0.37464 + 1.54226*acc_factor - 0.26992*acc_factor*acc_factor; + + for (int i = 0; i<10; i++) + { + const double temp = 296 + i * 2; + set_r(0.999); + test_phase->setState_TR(temp, 1.0737); + mv = 1 / rho * test_phase->meanMolecularWeight(); + //Calculate pressure using Peng-Robinson EoS + alpha = pow(1 + kappa*(1 - sqrt(temp / Tcrit)), 2); + pres_theoretical = GasConstant*temp / (mv - b_coeff) - a_coeff*alpha / (mv*mv + 2*b_coeff*mv - b_coeff*b_coeff); + EXPECT_NEAR(test_phase->pressure(), pres_theoretical, 2); + } +} +};