Cantera  3.1.0a1
Loading...
Searching...
No Matches
PengRobinson.h
Go to the documentation of this file.
1//! @file PengRobinson.h
2
3// This file is part of Cantera. See License.txt in the top-level directory or
4// at https://cantera.org/license.txt for license and copyright information.
5
6#ifndef CT_PENGROBINSON_H
7#define CT_PENGROBINSON_H
8
9#include "MixtureFugacityTP.h"
10#include "cantera/base/Array.h"
11
12namespace Cantera
13{
14/**
15 * Implementation of a multi-species Peng-Robinson equation of state
16 *
17 * @ingroup thermoprops
18 */
20{
21public:
22
23 //! Construct and initialize a PengRobinson object directly from an
24 //! input file
25 /*!
26 * @param infile Name of the input file containing the phase YAML data.
27 * If blank, an empty phase will be created.
28 * @param id ID of the phase in the input file. If empty, the
29 * first phase definition in the input file will be used.
30 */
31 explicit PengRobinson(const string& infile="", const string& id="");
32
33 string type() const override {
34 return "Peng-Robinson";
35 }
36
37 //! @name Molar Thermodynamic properties
38 //! @{
39
40 double cp_mole() const override;
41 double cv_mole() const override;
42
43 //! @}
44 //! @name Mechanical Properties
45 //! @{
46
47 //! Return the thermodynamic pressure (Pa).
48 /*!
49 * Since the mass density, temperature, and mass fractions are stored,
50 * this method uses these values to implement the
51 * mechanical equation of state @f$ P(T, \rho, Y_1, \dots, Y_K) @f$.
52 *
53 * @f[
54 * P = \frac{RT}{v-b_{mix}}
55 * - \frac{\left(\alpha a\right)_{mix}}{v^2 + 2b_{mix}v - b_{mix}^2}
56 * @f]
57 *
58 * where:
59 *
60 * @f[
61 * \alpha = \left[ 1 + \kappa \left(1-T_r^{0.5}\right)\right]^2
62 * @f]
63 *
64 * and
65 *
66 * @f[
67 * \kappa = \left(0.37464 + 1.54226\omega - 0.26992\omega^2\right),
68 * \qquad \qquad \text{For } \omega <= 0.491 \\
69 *
70 * \kappa = \left(0.379642 + 1.487503\omega - 0.164423\omega^2 + 0.016667\omega^3 \right),
71 * \qquad \text{For } \omega > 0.491
72 * @f]
73 *
74 * Coefficients @f$ a_{mix}, b_{mix} @f$ and @f$ (a \alpha)_{mix} @f$ are calculated as
75 *
76 * @f[
77 * 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}
78 * @f]
79 *
80 * @f[
81 * b_{mix} = \sum_i X_i b_i
82 * @f]
83 *
84 * @f[
85 * {a \alpha}_{mix} = \sum_i \sum_j X_i X_j {a \alpha}_{i, j}
86 * = \sum_i \sum_j X_i X_j \sqrt{a_i a_j} \sqrt{\alpha_i \alpha_j}
87 * @f]
88 */
89 double pressure() const override;
90
91 //! @}
92
93 //! Returns the standard concentration @f$ C^0_k @f$, which is used to
94 //! normalize the generalized concentration.
95 /*!
96 * This is defined as the concentration by which the generalized
97 * concentration is normalized to produce the activity.
98 * The ideal gas mixture is considered as the standard or reference state here.
99 * Since the activity for an ideal gas mixture is simply the mole fraction,
100 * for an ideal gas, @f$ C^0_k = P/\hat R T @f$.
101 *
102 * @param k Optional parameter indicating the species. The default is to
103 * assume this refers to species 0.
104 * @return
105 * Returns the standard Concentration in units of m^3 / kmol.
106 */
107 double standardConcentration(size_t k=0) const override;
108
109 //! Get the array of non-dimensional activity coefficients at the current
110 //! solution temperature, pressure, and solution concentration.
111 /*!
112 * For all objects with the Mixture Fugacity approximation, we define the
113 * standard state as an ideal gas at the current temperature and pressure of
114 * the solution. The activities are based on this standard state.
115 *
116 * @param ac Output vector of activity coefficients. Length: m_kk.
117 */
118 void getActivityCoefficients(double* ac) const override;
119
120 //! @name Partial Molar Properties of the Solution
121 //! @{
122
123 void getChemPotentials(double* mu) const override;
124 void getPartialMolarEnthalpies(double* hbar) const override;
125 void getPartialMolarEntropies(double* sbar) const override;
126 void getPartialMolarIntEnergies(double* ubar) const override;
127 //! Calculate species-specific molar specific heats
128 /*!
129 * This function is currently not implemented for Peng-Robinson phase.
130 */
131 void getPartialMolarCp(double* cpbar) const override;
132 void getPartialMolarVolumes(double* vbar) const override;
133 //! @}
134
135 //! Calculate species-specific critical temperature
136 /*!
137 * The temperature dependent parameter in P-R EoS is calculated as
138 * @f[ T_{crit} = (0.0778 a)/(0.4572 b R) @f]
139 * Units: Kelvin
140 *
141 * @param a species-specific coefficients used in P-R EoS
142 * @param b species-specific coefficients used in P-R EoS
143 */
144 double speciesCritTemperature(double a, double b) const;
145
146 //! @name Initialization Methods - For Internal use
147 //!
148 //! The following methods are used in the process of constructing
149 //! the phase and setting its parameters from a specification in an
150 //! input file. They are not normally used in application programs.
151 //! To see how they are used, see importPhase().
152 //! @{
153
154 bool addSpecies(shared_ptr<Species> spec) override;
155 void initThermo() override;
156 void getSpeciesParameters(const string& name, AnyMap& speciesNode) const override;
157
158 //! Set the pure fluid interaction parameters for a species
159 /*!
160 * @param species Name of the species
161 * @param a @f$ a @f$ parameter in the Peng-Robinson model [Pa-m^6/kmol^2]
162 * @param b @f$ a @f$ parameter in the Peng-Robinson model [m^3/kmol]
163 * @param w acentric factor
164 */
165 void setSpeciesCoeffs(const string& species, double a, double b, double w);
166
167 //! Set values for the interaction parameter between two species
168 /*!
169 * @param species_i Name of one species
170 * @param species_j Name of the other species
171 * @param a @f$ a @f$ parameter in the Peng-Robinson model [Pa-m^6/kmol^2]
172 */
173 void setBinaryCoeffs(const string& species_i, const string& species_j, double a);
174 //! @}
175
176protected:
177 // Special functions inherited from MixtureFugacityTP
178 double sresid() const override;
179 double hresid() const override;
180
181 double liquidVolEst(double T, double& pres) const override;
182 double densityCalc(double T, double pressure, int phase, double rhoguess) override;
183
184 double densSpinodalLiquid() const override;
185 double densSpinodalGas() const override;
186 double dpdVCalc(double T, double molarVol, double& presCalc) const override;
187
188 // Special functions not inherited from MixtureFugacityTP
189
190 //! Calculate temperature derivative @f$ d(a \alpha)/dT @f$
191 /*!
192 * These are stored internally.
193 */
194 double daAlpha_dT() const;
195
196 //! Calculate second derivative @f$ d^2(a \alpha)/dT^2 @f$
197 /*!
198 * These are stored internally.
199 */
200 double d2aAlpha_dT2() const;
201
202public:
203
204 double isothermalCompressibility() const override;
205 double thermalExpansionCoeff() const override;
206 double soundSpeed() const override;
207
208 //! Calculate @f$ dp/dV @f$ and @f$ dp/dT @f$ at the current conditions
209 /*!
210 * These are stored internally.
211 */
212 void calculatePressureDerivatives() const;
213
214 //! Update the @f$ a @f$, @f$ b @f$, and @f$ \alpha @f$ parameters
215 /*!
216 * The @f$ a @f$ and the @f$ b @f$ parameters depend on the mole fraction and the
217 * parameter @f$ \alpha @f$ depends on the temperature. This function updates
218 * the internal numbers based on the state of the object.
219 */
220 void updateMixingExpressions() override;
221
222 //! Calculate the @f$ a @f$, @f$ b @f$, and @f$ \alpha @f$ parameters given the temperature
223 /*!
224 * This function doesn't change the internal state of the object, so it is a
225 * const function. It does use the stored mole fractions in the object.
226 *
227 * @param aCalc (output) Returns the a value
228 * @param bCalc (output) Returns the b value.
229 * @param aAlpha (output) Returns the (a*alpha) value.
230 */
231 void calculateAB(double& aCalc, double& bCalc, double& aAlpha) const;
232
233 void calcCriticalConditions(double& pc, double& tc, double& vc) const override;
234
235 //! Prepare variables and call the function to solve the cubic equation of state
236 int solveCubic(double T, double pres, double a, double b, double aAlpha,
237 double Vroot[3]) const;
238protected:
239 //! Value of @f$ b @f$ in the equation of state
240 /*!
241 * `m_b` is a function of the mole fractions and species-specific b values.
242 */
243 double m_b = 0.0;
244
245 //! Value of @f$ a @f$ in the equation of state
246 /*!
247 * `m_a` depends only on the mole fractions.
248 */
249 double m_a = 0.0;
250
251 //! Value of @f$ a \alpha @f$ in the equation of state
252 /*!
253 * `m_aAlpha_mix` is a function of the temperature and the mole fractions.
254 */
255 double m_aAlpha_mix = 0.0;
256
257 // Vectors required to store species-specific a_coeff, b_coeff, alpha, kappa
258 // and other derivatives. Length = m_kk.
259 vector<double> m_b_coeffs;
260 vector<double> m_kappa;
261 vector<double> m_acentric; //!< acentric factor for each species, length #m_kk
262 mutable vector<double> m_dalphadT;
263 mutable vector<double> m_d2alphadT2;
264 vector<double> m_alpha;
265
266 // Matrices for Binary coefficients a_{i,j} and {a*alpha}_{i.j} are saved in an
267 // array form. Size = (m_kk, m_kk).
268 Array2D m_a_coeffs;
269 Array2D m_aAlpha_binary;
270
271 //! Explicitly-specified binary interaction parameters, to enable serialization
272 map<string, map<string, double>> m_binaryParameters;
273
274 int m_NSolns = 0;
275
276 double m_Vroot[3] = {0.0, 0.0, 0.0};
277
278 //! Temporary storage - length = m_kk.
279 mutable vector<double> m_pp;
280
281 // Partial molar volumes of the species
282 mutable vector<double> m_partialMolarVolumes;
283
284 //! The derivative of the pressure with respect to the volume
285 /*!
286 * Calculated at the current conditions. temperature and mole number kept
287 * constant
288 */
289 mutable double m_dpdV = 0.0;
290
291 //! The derivative of the pressure with respect to the temperature
292 /*!
293 * Calculated at the current conditions. Total volume and mole number kept
294 * constant
295 */
296 mutable double m_dpdT = 0.0;
297
298 //! Vector of derivatives of pressure with respect to mole number
299 /*!
300 * Calculated at the current conditions. Total volume, temperature and
301 * other mole number kept constant
302 */
303 mutable vector<double> m_dpdni;
304
305 enum class CoeffSource { EoS, CritProps, Database };
306 //! For each species, specifies the source of the a, b, and omega coefficients
307 vector<CoeffSource> m_coeffSource;
308
309private:
310 //! Omega constant: a0 (= omega_a) used in Peng-Robinson equation of state
311 /*!
312 * This value is calculated by solving P-R cubic equation at the critical point.
313 */
314 static const double omega_a;
315
316 //! Omega constant: b0 (= omega_b) used in Peng-Robinson equation of state
317 /*!
318 * This value is calculated by solving P-R cubic equation at the critical point.
319 */
320 static const double omega_b;
321
322 //! Omega constant for the critical molar volume
323 static const double omega_vc;
324};
325}
326
327#endif
Header file for class Cantera::Array2D.
Header file for a derived class of ThermoPhase that handles non-ideal mixtures based on the fugacity ...
A map of string keys to values whose type can vary at runtime.
Definition AnyMap.h:427
A class for 2D arrays stored in column-major (Fortran-compatible) form.
Definition Array.h:32
This is a filter class for ThermoPhase that implements some preparatory steps for efficiently handlin...
Implementation of a multi-species Peng-Robinson equation of state.
double thermalExpansionCoeff() const override
Return the volumetric thermal expansion coefficient. Units: 1/K.
double densSpinodalLiquid() const override
Return the value of the density at the liquid spinodal point (on the liquid side) for the current tem...
void getPartialMolarEnthalpies(double *hbar) const override
Returns an array of partial molar enthalpies for the species in the mixture.
void getChemPotentials(double *mu) const override
Get the species chemical potentials. Units: J/kmol.
void setSpeciesCoeffs(const string &species, double a, double b, double w)
Set the pure fluid interaction parameters for a species.
double sresid() const override
Calculate the deviation terms for the total entropy of the mixture from the ideal gas mixture.
double soundSpeed() const override
Return the speed of sound. Units: m/s.
double pressure() const override
Return the thermodynamic pressure (Pa).
map< string, map< string, double > > m_binaryParameters
Explicitly-specified binary interaction parameters, to enable serialization.
void getSpeciesParameters(const string &name, AnyMap &speciesNode) const override
Get phase-specific parameters of a Species object such that an identical one could be reconstructed a...
static const double omega_b
Omega constant: b0 (= omega_b) used in Peng-Robinson equation of state.
vector< double > m_pp
Temporary storage - length = m_kk.
void setBinaryCoeffs(const string &species_i, const string &species_j, double a)
Set values for the interaction parameter between two species.
vector< double > m_dpdni
Vector of derivatives of pressure with respect to mole number.
string type() const override
String indicating the thermodynamic model implemented.
void initThermo() override
Initialize the ThermoPhase object after all species have been set up.
double densSpinodalGas() const override
Return the value of the density at the gas spinodal point (on the gas side) for the current temperatu...
double m_dpdT
The derivative of the pressure with respect to the temperature.
void getPartialMolarVolumes(double *vbar) const override
Return an array of partial molar volumes for the species in the mixture.
double cv_mole() const override
Molar heat capacity at constant volume. Units: J/kmol/K.
double liquidVolEst(double T, double &pres) const override
Estimate for the molar volume of the liquid.
int solveCubic(double T, double pres, double a, double b, double aAlpha, double Vroot[3]) const
Prepare variables and call the function to solve the cubic equation of state.
double daAlpha_dT() const
Calculate temperature derivative .
void calculatePressureDerivatives() const
Calculate and at the current conditions.
double isothermalCompressibility() const override
Returns the isothermal compressibility. Units: 1/Pa.
double hresid() const override
Calculate the deviation terms for the total enthalpy of the mixture from the ideal gas mixture.
double m_aAlpha_mix
Value of in the equation of state.
double m_a
Value of in the equation of state.
static const double omega_a
Omega constant: a0 (= omega_a) used in Peng-Robinson equation of state.
static const double omega_vc
Omega constant for the critical molar volume.
vector< CoeffSource > m_coeffSource
For each species, specifies the source of the a, b, and omega coefficients.
void calculateAB(double &aCalc, double &bCalc, double &aAlpha) const
Calculate the , , and parameters given the temperature.
void getPartialMolarIntEnergies(double *ubar) const override
Return an array of partial molar internal energies for the species in the mixture.
void updateMixingExpressions() override
Update the , , and parameters.
double cp_mole() const override
Molar heat capacity at constant pressure. Units: J/kmol/K.
double m_b
Value of in the equation of state.
void getPartialMolarCp(double *cpbar) const override
Calculate species-specific molar specific heats.
double standardConcentration(size_t k=0) const override
Returns the standard concentration , which is used to normalize the generalized concentration.
bool addSpecies(shared_ptr< Species > spec) override
Add a Species to this Phase.
double dpdVCalc(double T, double molarVol, double &presCalc) const override
Calculate the pressure and the pressure derivative given the temperature and the molar volume.
void getActivityCoefficients(double *ac) const override
Get the array of non-dimensional activity coefficients at the current solution temperature,...
double d2aAlpha_dT2() const
Calculate second derivative .
double densityCalc(double T, double pressure, int phase, double rhoguess) override
Calculates the density given the temperature and the pressure and a guess at the density.
void getPartialMolarEntropies(double *sbar) const override
Returns an array of partial molar entropies of the species in the solution.
double m_dpdV
The derivative of the pressure with respect to the volume.
double speciesCritTemperature(double a, double b) const
Calculate species-specific critical temperature.
vector< double > m_acentric
acentric factor for each species, length m_kk
shared_ptr< Species > species(const string &name) const
Return the Species object for the named species.
Definition Phase.cpp:856
string name() const
Return the name of the phase.
Definition Phase.cpp:20
Namespace for the Cantera kernel.
Definition AnyMap.cpp:564