12#include <boost/algorithm/string.hpp>
13#include <boost/math/tools/roots.hpp>
15namespace bmt = boost::math::tools;
34 "Unknown species '{}'.",
species);
40 m_kappa[k] = 0.37464 + 1.54226*w - 0.26992*w*w;
42 m_kappa[k] = 0.374642 + 1.487503*w - 0.164423*w*w + 0.016666*w*w*w;
49 double sqt_alpha = 1 + m_kappa[k] * (1 - sqt_T_r);
50 m_alpha[k] = sqt_alpha*sqt_alpha;
53 double aAlpha_k = a*m_alpha[k];
54 m_aAlpha_binary(k,k) = aAlpha_k;
57 for (
size_t j = 0; j <
m_kk; j++) {
61 double a0kj = sqrt(m_a_coeffs(j,j) * a);
62 double aAlpha_j = a*m_alpha[j];
63 double a_Alpha = sqrt(aAlpha_j*aAlpha_k);
64 if (m_a_coeffs(j, k) == 0) {
65 m_a_coeffs(j, k) = a0kj;
66 m_aAlpha_binary(j, k) = a_Alpha;
67 m_a_coeffs(k, j) = a0kj;
68 m_aAlpha_binary(k, j) = a_Alpha;
75 const string& species_j,
double a0)
80 "Unknown species '{}'.", species_i);
85 "Unknown species '{}'.", species_j);
88 m_a_coeffs(ki, kj) = m_a_coeffs(kj, ki) = a0;
92 double alpha_ij = m_alpha[ki] * m_alpha[kj];
93 m_aAlpha_binary(ki, kj) = m_aAlpha_binary(kj, ki) = a0*alpha_ij;
103 double vpb = mv + (1 +
Sqrt2) *
m_b;
104 double vmb = mv + (1 -
Sqrt2) *
m_b;
126 double denom = mv * mv + 2 * mv *
m_b -
m_b *
m_b;
139 double vpb2 = mv + (1 +
Sqrt2) *
m_b;
140 double vmb2 = mv + (1 -
Sqrt2) *
m_b;
141 double vmb = mv -
m_b;
144 for (
size_t k = 0; k <
m_kk; k++) {
146 for (
size_t i = 0; i <
m_kk; i++) {
152 double denom2 =
m_b * (mv * mv + 2 * mv *
m_b -
m_b *
m_b);
154 for (
size_t k = 0; k <
m_kk; k++) {
156 ac[k] = (-RT_ * log(pres * mv/ RT_) + RT_ * log(mv / vmb)
157 + RT_ * m_b_coeffs[k] / vmb
158 - (num /denom) * log(vpb2/vmb2)
162 for (
size_t k = 0; k <
m_kk; k++) {
163 ac[k] = exp(ac[k]/ RT_);
173 for (
size_t k = 0; k <
m_kk; k++) {
175 mu[k] += RT_ * (log(xx));
179 double vmb = mv -
m_b;
180 double vpb2 = mv + (1 +
Sqrt2) *
m_b;
181 double vmb2 = mv + (1 -
Sqrt2) *
m_b;
183 for (
size_t k = 0; k <
m_kk; k++) {
185 for (
size_t i = 0; i <
m_kk; i++) {
192 double denom2 =
m_b * (mv * mv + 2 * mv *
m_b -
m_b *
m_b);
194 for (
size_t k = 0; k <
m_kk; k++) {
197 mu[k] += RT_ * log(pres/refP) - RT_ * log(pres * mv / RT_)
198 + RT_ * log(mv / vmb) + RT_ * m_b_coeffs[k] / vmb
199 - (num /denom) * log(vpb2/vmb2)
210 tmp.resize(
m_kk,0.0);
215 double vmb = mv -
m_b;
216 double vpb2 = mv + (1 +
Sqrt2) *
m_b;
217 double vmb2 = mv + (1 -
Sqrt2) *
m_b;
220 for (
size_t k = 0; k <
m_kk; k++) {
223 for (
size_t i = 0; i <
m_kk; i++) {
224 double grad_aAlpha = m_dalphadT[i]/m_alpha[i] + m_dalphadT[k]/m_alpha[k];
230 double denom = mv * mv + 2 * mv *
m_b -
m_b *
m_b;
231 double denom2 = denom * denom;
233 for (
size_t k = 0; k <
m_kk; k++) {
234 m_dpdni[k] = RT_ / vmb + RT_ * m_b_coeffs[k] / (vmb * vmb)
243 for (
size_t k = 0; k <
m_kk; k++) {
244 fac4 = T*tmp[k] -2 *
m_pp[k];
245 double hE_v = mv *
m_dpdni[k] - RT_
246 - m_b_coeffs[k] / fac3 * log(vpb2 / vmb2) * fac
247 + (mv * m_b_coeffs[k]) / (
m_b * denom) * fac
248 + 1 / (2 *
Sqrt2 *
m_b) * log(vpb2 / vmb2) * fac4;
249 hbar[k] = hbar[k] + hE_v;
260 for (
size_t k = 0; k <
m_kk; k++) {
261 sbar[k] = (sbar[k] -
m_tmpV[k])/T;
271 for (
size_t k = 0; k <
m_kk; k++) {
272 ubar[k] = ubar[k] - p*
m_tmpV[k];
283 for (
size_t k = 0; k <
m_kk; k++) {
285 for (
size_t i = 0; i <
m_kk; i++) {
291 double vmb = mv -
m_b;
292 double vpb = mv +
m_b;
293 double fac = mv * mv + 2 * mv *
m_b -
m_b *
m_b;
294 double fac2 = fac * fac;
297 for (
size_t k = 0; k <
m_kk; k++) {
298 double num = RT_ + RT_ *
m_b/ vmb + RT_ * m_b_coeffs[k] / vmb
299 + RT_ *
m_b * m_b_coeffs[k] /(vmb * vmb) - 2 * mv *
m_pp[k] / fac
303 vbar[k] = num / denom;
311 }
else if (a <= 0.0) {
323 m_b_coeffs.push_back(0.0);
325 m_kappa.push_back(0.0);
327 m_alpha.push_back(0.0);
328 m_dalphadT.push_back(0.0);
329 m_d2alphadT2.push_back(0.0);
331 m_partialMolarVolumes.push_back(0.0);
342 std::unordered_map<string, AnyMap*> dbSpecies;
347 if (m_a_coeffs(k, k) != 0.0) {
350 bool foundCoeffs =
false;
351 if (data.hasKey(
"equation-of-state") &&
352 data[
"equation-of-state"].hasMapWhere(
"model",
"Peng-Robinson"))
356 auto eos = data[
"equation-of-state"].getMapWhere(
357 "model",
"Peng-Robinson");
358 if (eos.hasKey(
"a") && eos.hasKey(
"b") && eos.hasKey(
"acentric-factor")) {
359 double a0 = eos.convert(
"a",
"Pa*m^6/kmol^2");
360 double b = eos.convert(
"b",
"m^3/kmol");
362 double w = eos[
"acentric-factor"].asDouble();
367 if (eos.hasKey(
"binary-a")) {
370 for (
auto& [name2, coeff] : binary_a) {
371 double a0 = units.
convert(coeff,
"Pa*m^6/kmol^2");
382 double Tc = NAN, Pc = NAN, omega_ac = NAN;
383 if (data.hasKey(
"critical-parameters")) {
386 auto& critProps = data[
"critical-parameters"].as<
AnyMap>();
387 Tc = critProps.
convert(
"critical-temperature",
"K");
388 Pc = critProps.convert(
"critical-pressure",
"Pa");
389 omega_ac = critProps[
"acentric-factor"].asDouble();
393 if (critPropsDb.
empty()) {
395 dbSpecies = critPropsDb[
"species"].asMap(
"name");
399 auto ucName = boost::algorithm::to_upper_copy(
name);
400 if (dbSpecies.count(ucName)) {
401 auto& spec = *dbSpecies.at(ucName);
402 auto& critProps = spec[
"critical-parameters"].as<
AnyMap>();
403 Tc = critProps.
convert(
"critical-temperature",
"K");
404 Pc = critProps.convert(
"critical-pressure",
"Pa");
405 omega_ac = critProps[
"acentric-factor"].asDouble();
417 "No Peng-Robinson model parameters or critical properties found for "
418 "species '{}'",
name);
431 auto& eosNode = speciesNode[
"equation-of-state"].getMapWhere(
432 "model",
"Peng-Robinson",
true);
433 eosNode[
"a"].setQuantity(m_a_coeffs(k, k),
"Pa*m^6/kmol^2");
434 eosNode[
"b"].setQuantity(m_b_coeffs[k],
"m^3/kmol");
437 auto& critProps = speciesNode[
"critical-parameters"];
440 critProps[
"critical-temperature"].setQuantity(Tc,
"K");
441 critProps[
"critical-pressure"].setQuantity(Pc,
"Pa");
449 auto& eosNode = speciesNode[
"equation-of-state"].getMapWhere(
450 "model",
"Peng-Robinson",
true);
453 bin_a[
species].setQuantity(coeff,
"Pa*m^6/kmol^2");
455 eosNode[
"binary-a"] = std::move(bin_a);
462 double hh =
m_b / molarV;
465 double vpb = molarV + (1.0 +
Sqrt2) *
m_b;
466 double vmb = molarV + (1.0 -
Sqrt2) *
m_b;
467 double fac = alpha_1 / (2.0 *
Sqrt2 *
m_b);
468 double sresid_mol_R = log(zz*(1.0 - hh)) + fac * log(vpb / vmb) /
GasConstant;
478 double vpb = molarV + (1 +
Sqrt2) *
m_b;
479 double vmb = molarV + (1 -
Sqrt2) *
m_b;
480 double fac = 1 / (2.0 *
Sqrt2 *
m_b);
487 double v =
m_b * 1.1;
488 double atmp, btmp, aAlphatmp;
490 double pres = std::max(
psatEst(T), presGuess);
492 bool foundLiq =
false;
494 while (m < 100 && !foundLiq) {
495 int nsol =
solveCubic(T, pres, atmp, btmp, aAlphatmp, Vroot);
496 if (nsol == 1 || nsol == 2) {
523 if (rhoGuess == -1.0) {
524 if (phaseRequested >= FLUID_LIQUID_0) {
526 rhoGuess = mmw / lqvol;
533 double volGuess = mmw / rhoGuess;
536 double molarVolLast = m_Vroot[0];
538 if (phaseRequested >= FLUID_LIQUID_0) {
539 molarVolLast = m_Vroot[0];
540 }
else if (phaseRequested == FLUID_GAS || phaseRequested == FLUID_SUPERCRIT) {
541 molarVolLast = m_Vroot[2];
543 if (volGuess > m_Vroot[1]) {
544 molarVolLast = m_Vroot[2];
546 molarVolLast = m_Vroot[0];
549 }
else if (m_NSolns == 1) {
550 if (phaseRequested == FLUID_GAS || phaseRequested == FLUID_SUPERCRIT
551 || phaseRequested == FLUID_UNDEFINED)
553 molarVolLast = m_Vroot[0];
557 }
else if (m_NSolns == -1) {
558 if (phaseRequested >= FLUID_LIQUID_0 || phaseRequested == FLUID_UNDEFINED
559 || phaseRequested == FLUID_SUPERCRIT)
561 molarVolLast = m_Vroot[0];
562 }
else if (T > tcrit) {
563 molarVolLast = m_Vroot[0];
568 molarVolLast = m_Vroot[0];
571 return mmw / molarVolLast;
583 auto resid = [
this, T](
double v) {
588 boost::uintmax_t maxiter = 100;
589 auto [lower, upper] = bmt::toms748_solve(
590 resid, Vroot[0], Vroot[1], bmt::eps_tolerance<double>(48), maxiter);
593 return mmw / (0.5 * (lower + upper));
605 auto resid = [
this, T](
double v) {
610 boost::uintmax_t maxiter = 100;
611 auto [lower, upper] = bmt::toms748_solve(
612 resid, Vroot[1], Vroot[2], bmt::eps_tolerance<double>(48), maxiter);
615 return mmw / (0.5 * (lower + upper));
620 double denom = molarVol * molarVol + 2 * molarVol *
m_b -
m_b *
m_b;
621 double vpb = molarVol +
m_b;
622 double vmb = molarVol -
m_b;
651 double vmb = mv -
m_b;
652 double denom = mv * mv + 2 * mv *
m_b -
m_b *
m_b;
661 for (
size_t j = 0; j <
m_kk; j++) {
663 double sqt_alpha = 1 + m_kappa[j] * (1 - sqrt(temp / critTemp_j));
664 m_alpha[j] = sqt_alpha*sqt_alpha;
668 for (
size_t i = 0; i <
m_kk; i++) {
669 for (
size_t j = 0; j <
m_kk; j++) {
670 m_aAlpha_binary(i, j) = sqrt(m_alpha[i] * m_alpha[j]) * m_a_coeffs(i,j);
681 for (
size_t i = 0; i <
m_kk; i++) {
683 for (
size_t j = 0; j <
m_kk; j++) {
692 double daAlphadT = 0.0, k, Tc, sqtTr, coeff1, coeff2;
693 for (
size_t i = 0; i <
m_kk; i++) {
697 coeff1 = 1 / (Tc*sqtTr);
700 m_dalphadT[i] = coeff1 * (k*k*coeff2 - k);
703 for (
size_t i = 0; i <
m_kk; i++) {
704 for (
size_t j = 0; j <
m_kk; j++) {
706 * m_aAlpha_binary(i, j)
707 * (m_dalphadT[i] / m_alpha[i] + m_dalphadT[j] / m_alpha[j]);
715 for (
size_t i = 0; i <
m_kk; i++) {
718 double coeff1 = 1 / (Tcrit_i*sqt_Tr);
719 double coeff2 = sqt_Tr - 1;
721 double k = m_kappa[i];
722 m_dalphadT[i] = coeff1 * (k*k*coeff2 - k);
723 m_d2alphadT2[i] = (k*k + k) * coeff1 / (2*sqt_Tr*sqt_Tr*Tcrit_i);
727 double d2aAlphadT2 = 0.0;
728 for (
size_t i = 0; i <
m_kk; i++) {
729 double alphai = m_alpha[i];
730 for (
size_t j = 0; j <
m_kk; j++) {
731 double alphaj = m_alpha[j];
732 double alphaij = alphai * alphaj;
733 double term1 = m_d2alphadT2[i] / alphai + m_d2alphadT2[j] / alphaj;
734 double term2 = 2 * m_dalphadT[i] * m_dalphadT[j] / alphaij;
735 double term3 = m_dalphadT[i] / alphai + m_dalphadT[j] / alphaj;
737 * m_aAlpha_binary(i, j)
738 * (term1 + term2 - 0.5 * term3 * term3);
744void PengRobinson::calcCriticalConditions(
double& pc,
double& tc,
double& vc)
const
764 double Vroot[3])
const
769 double aAlpha_p = aAlpha / pres;
771 double bn = (b - RT_p);
772 double cn = -(2 * RT_p * b - aAlpha_p + 3 * bsqr);
773 double dn = (bsqr * RT_p + bsqr * b - aAlpha_p * b);
780 an, bn, cn, dn, tc, vc);
Declaration for class Cantera::Species.
Headers for the factory class that can create known ThermoPhase objects (see Thermodynamic Properties...
A map of string keys to values whose type can vary at runtime.
const UnitSystem & units() const
Return the default units that should be used to convert stored values.
bool empty() const
Return boolean indicating whether AnyMap is empty.
double convert(const string &key, const string &units) const
Convert the item stored by the given key to the units specified in units.
static AnyMap fromYamlFile(const string &name, const string &parent_name="")
Create an AnyMap from a YAML file.
virtual void resize(size_t n, size_t m, double v=0.0)
Resize the array, and fill the new entries with 'v'.
Base class for exceptions thrown by Cantera classes.
double critPressure() const override
Critical pressure (Pa).
double critDensity() const override
Critical density (kg/m3).
void getGibbs_ref(double *g) const override
Returns the vector of the Gibbs function of the reference state at the current temperature of the sol...
double critTemperature() const override
Critical temperature (K).
virtual void _updateReferenceStateThermo() const
Updates the reference state thermodynamic functions at the current T of the solution.
vector< double > m_tmpV
Temporary storage - length = m_kk.
int solveCubic(double T, double pres, double a, double b, double aAlpha, double Vroot[3], double an, double bn, double cn, double dn, double tc, double vc) const
Solve the cubic equation of state.
vector< double > moleFractions_
Storage for the current values of the mole fractions of the species.
void setTemperature(const double temp) override
Set the temperature of the phase.
virtual double psatEst(double TKelvin) const
Estimate for the saturation pressure.
void getStandardVolumes(double *vol) const override
Get the molar volumes of each species in their standard states at the current T and P of the solution...
vector< double > m_cp0_R
Temporary storage for dimensionless reference state heat capacities.
bool addSpecies(shared_ptr< Species > spec) override
Add a Species to this Phase.
double z() const
Calculate the value of z.
void getEnthalpy_RT_ref(double *hrt) const override
Returns the vector of nondimensional enthalpies of the reference state at the current temperature of ...
An error indicating that an unimplemented function has been called.
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.
PengRobinson(const string &infile="", const string &id="")
Construct and initialize a PengRobinson object directly from an input file.
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
void checkSpeciesIndex(size_t k) const
Check that the specified species index is in range.
size_t m_kk
Number of species in the phase.
double temperature() const
Temperature (K).
double meanMolecularWeight() const
The mean molecular weight. Units: (kg/kmol)
size_t speciesIndex(const string &name) const
Returns the index of a species named 'name' within the Phase object.
double moleFraction(size_t k) const
Return the mole fraction of a single species.
double mean_X(const double *const Q) const
Evaluate the mole-fraction-weighted mean of an array Q.
shared_ptr< Species > species(const string &name) const
Return the Species object for the named species.
virtual double molarVolume() const
Molar volume (m^3/kmol).
string name() const
Return the name of the phase.
double RT() const
Return the Gas Constant multiplied by the current temperature.
void initThermoFile(const string &inputFile, const string &id)
Initialize a ThermoPhase object using an input file.
virtual void getSpeciesParameters(const string &name, AnyMap &speciesNode) const
Get phase-specific parameters of a Species object such that an identical one could be reconstructed a...
virtual double refPressure() const
Returns the reference pressure in Pa.
double convert(double value, const string &src, const string &dest) const
Convert value from the units of src to the units of dest.
void scale(InputIter begin, InputIter end, OutputIter out, S scale_factor)
Multiply elements of an array by a scale factor.
const double GasConstant
Universal Gas Constant [J/kmol/K].
const double Sqrt2
Sqrt(2)
Namespace for the Cantera kernel.
const size_t npos
index returned by functions to indicate "no position"
const double SmallNumber
smallest number to compare to zero.
Contains declarations for string manipulation functions within Cantera.
Various templated functions that carry out common vector and polynomial operations (see Templated Arr...