12#include <boost/algorithm/string.hpp>
13#include <boost/math/tools/roots.hpp>
15namespace bmt = boost::math::tools;
30 double a0,
double a1,
double b)
35 "Unknown species '{}'.",
species);
42 size_t counter = k +
m_kk * k;
43 a_coeff_vec(0, counter) = a0;
44 a_coeff_vec(1, counter) = a1;
47 for (
size_t j = 0; j <
m_kk; j++) {
53 if (isnan(a_coeff_vec(0, j +
m_kk * j))) {
56 }
else if (isnan(a_coeff_vec(0, j +
m_kk * k))) {
59 double a0kj = sqrt(a_coeff_vec(0, j +
m_kk * j) * a0);
60 double a1kj = sqrt(a_coeff_vec(1, j +
m_kk * j) * a1);
61 a_coeff_vec(0, j +
m_kk * k) = a0kj;
62 a_coeff_vec(1, j +
m_kk * k) = a1kj;
63 a_coeff_vec(0, k +
m_kk * j) = a0kj;
64 a_coeff_vec(1, k +
m_kk * j) = a1kj;
67 a_coeff_vec.
getRow(0, a_vec_Curr_.data());
77 "Unknown species '{}'.", species_i);
82 "Unknown species '{}'.", species_j);
91 size_t counter1 = ki +
m_kk * kj;
92 size_t counter2 = kj +
m_kk * ki;
93 a_coeff_vec(0, counter1) = a_coeff_vec(0, counter2) = a0;
94 a_coeff_vec(1, counter1) = a_coeff_vec(1, counter2) = a1;
95 a_vec_Curr_[counter1] = a_vec_Curr_[counter2] = a0;
104 double sqt = sqrt(TKelvin);
109 double dadt = da_dt();
110 double fac = TKelvin * dadt - 3.0 *
m_a_current / 2.0;
112 +1.0/(
m_b_current * sqt) * log(vpb/mv) * (-0.5 * dadt));
120 double sqt = sqrt(TKelvin);
124 double dadt = da_dt();
125 double fac = TKelvin * dadt - 3.0 *
m_a_current / 2.0;
126 return (cvref - 1.0/(2.0 *
m_b_current * TKelvin * sqt) * log(vpb/mv)*fac
127 +1.0/(
m_b_current * sqt) * log(vpb/mv)*(-0.5*dadt));
154 for (
size_t k = 0; k <
m_kk; k++) {
156 for (
size_t i = 0; i <
m_kk; i++) {
157 size_t counter = k +
m_kk*i;
163 for (
size_t k = 0; k <
m_kk; k++) {
164 ac[k] = (-
RT() * log(pres * mv /
RT())
165 +
RT() * log(mv / vmb)
166 +
RT() * b_vec_Curr_[k] / vmb
172 for (
size_t k = 0; k <
m_kk; k++) {
173 ac[k] = exp(ac[k]/
RT());
182 for (
size_t k = 0; k <
m_kk; k++) {
184 mu[k] +=
RT()*(log(xx));
192 for (
size_t k = 0; k <
m_kk; k++) {
194 for (
size_t i = 0; i <
m_kk; i++) {
195 size_t counter = k +
m_kk*i;
202 for (
size_t k = 0; k <
m_kk; k++) {
203 mu[k] += (
RT() * log(pres/refP) -
RT() * log(pres * mv /
RT())
204 +
RT() * log(mv / vmb)
205 +
RT() * b_vec_Curr_[k] / vmb
222 double sqt = sqrt(TKelvin);
225 for (
size_t k = 0; k <
m_kk; k++) {
227 for (
size_t i = 0; i <
m_kk; i++) {
228 size_t counter = k +
m_kk*i;
232 for (
size_t k = 0; k <
m_kk; k++) {
233 dpdni_[k] =
RT()/vmb +
RT() * b_vec_Curr_[k] / (vmb * vmb) - 2.0 *
m_pp[k] / (sqt * mv * vpb)
234 +
m_a_current * b_vec_Curr_[k]/(sqt * mv * vpb * vpb);
236 double dadt = da_dt();
237 double fac = TKelvin * dadt - 3.0 *
m_a_current / 2.0;
239 for (
size_t k = 0; k <
m_kk; k++) {
241 for (
size_t i = 0; i <
m_kk; i++) {
242 size_t counter = k +
m_kk*i;
249 for (
size_t k = 0; k <
m_kk; k++) {
252 + b_vec_Curr_[k] / vpb / (
m_b_current * sqt) * fac);
253 hbar[k] = hbar[k] + hE_v;
254 hbar[k] -= fac2 *
dpdni_[k];
263 double sqt = sqrt(TKelvin);
267 for (
size_t k = 0; k <
m_kk; k++) {
271 for (
size_t k = 0; k <
m_kk; k++) {
273 for (
size_t i = 0; i <
m_kk; i++) {
274 size_t counter = k +
m_kk*i;
278 for (
size_t k = 0; k <
m_kk; k++) {
280 for (
size_t i = 0; i <
m_kk; i++) {
281 size_t counter = k +
m_kk*i;
286 double dadt = da_dt();
290 for (
size_t k = 0; k <
m_kk; k++) {
298 - 1.0 / (
m_b_current * sqt) * b_vec_Curr_[k] / vpb * fac
304 for (
size_t k = 0; k <
m_kk; k++) {
305 sbar[k] -= -m_partialMolarVolumes[k] *
dpdT_;
315 for (
size_t k = 0; k <
nSpecies(); k++) {
316 ubar[k] -= p * m_partialMolarVolumes[k];
322 for (
size_t k = 0; k <
m_kk; k++) {
324 for (
size_t i = 0; i <
m_kk; i++) {
325 size_t counter = k +
m_kk*i;
329 for (
size_t k = 0; k <
m_kk; k++) {
331 for (
size_t i = 0; i <
m_kk; i++) {
332 size_t counter = k +
m_kk*i;
341 for (
size_t k = 0; k <
m_kk; k++) {
344 - 2.0 *
m_pp[k] / (sqt * vpb)
349 vbar[k] = num / denom;
357 a_vec_Curr_.resize(
m_kk *
m_kk, 0.0);
361 b_vec_Curr_.push_back(NAN);
366 m_partialMolarVolumes.push_back(0.0);
376 std::unordered_map<string, AnyMap*> dbSpecies;
381 if (!isnan(a_coeff_vec(0, k +
m_kk * k))) {
384 bool foundCoeffs =
false;
386 if (data.hasKey(
"equation-of-state") &&
387 data[
"equation-of-state"].hasMapWhere(
"model",
"Redlich-Kwong"))
391 auto eos = data[
"equation-of-state"].getMapWhere(
392 "model",
"Redlich-Kwong");
394 if (eos.hasKey(
"a") && eos.hasKey(
"b")) {
395 double a0 = 0, a1 = 0;
396 if (eos[
"a"].isScalar()) {
397 a0 = eos.convert(
"a",
"Pa*m^6/kmol^2*K^0.5");
399 auto avec = eos[
"a"].asVector<
AnyValue>(2);
400 a0 = eos.units().convert(avec[0],
"Pa*m^6/kmol^2*K^0.5");
401 a1 = eos.units().convert(avec[1],
"Pa*m^6/kmol^2/K^0.5");
403 double b = eos.convert(
"b",
"m^3/kmol");
409 if (eos.hasKey(
"binary-a")) {
412 for (
auto& [name2, coeff] : binary_a) {
413 double a0 = 0, a1 = 0;
414 if (coeff.isScalar()) {
415 a0 = units.
convert(coeff,
"Pa*m^6/kmol^2*K^0.5");
417 auto avec = coeff.asVector<
AnyValue>(2);
418 a0 = units.convert(avec[0],
"Pa*m^6/kmol^2*K^0.5");
419 a1 = units.convert(avec[1],
"Pa*m^6/kmol^2/K^0.5");
431 double Tc = NAN, Pc = NAN;
432 if (data.hasKey(
"critical-parameters")) {
435 auto& critProps = data[
"critical-parameters"].as<
AnyMap>();
436 Tc = critProps.
convert(
"critical-temperature",
"K");
437 Pc = critProps.convert(
"critical-pressure",
"Pa");
441 if (critPropsDb.
empty()) {
443 dbSpecies = critPropsDb[
"species"].asMap(
"name");
447 auto ucName = boost::algorithm::to_upper_copy(
name);
448 if (dbSpecies.count(ucName)) {
449 auto& spec = *dbSpecies.at(ucName);
450 auto& critProps = spec[
"critical-parameters"].as<
AnyMap>();
451 Tc = critProps.
convert(
"critical-temperature",
"K");
452 Pc = critProps.convert(
"critical-pressure",
"Pa");
465 "No critical property or Redlich-Kwong parameters found "
466 "for species {}.",
name);
472 AnyMap& speciesNode)
const
478 auto& eosNode = speciesNode[
"equation-of-state"].getMapWhere(
479 "model",
"Redlich-Kwong",
true);
481 size_t counter = k +
m_kk * k;
482 if (a_coeff_vec(1, counter) != 0.0) {
483 vector<AnyValue> coeffs(2);
484 coeffs[0].setQuantity(a_coeff_vec(0, counter),
"Pa*m^6/kmol^2*K^0.5");
485 coeffs[1].setQuantity(a_coeff_vec(1, counter),
"Pa*m^6/kmol^2/K^0.5");
486 eosNode[
"a"] = std::move(coeffs);
488 eosNode[
"a"].setQuantity(a_coeff_vec(0, counter),
489 "Pa*m^6/kmol^2*K^0.5");
491 eosNode[
"b"].setQuantity(b_vec_Curr_[k],
"m^3/kmol");
493 auto& critProps = speciesNode[
"critical-parameters"];
494 double a = a_coeff_vec(0, k +
m_kk * k);
495 double b = b_vec_Curr_[k];
498 critProps[
"critical-temperature"].setQuantity(Tc,
"K");
499 critProps[
"critical-pressure"].setQuantity(Pc,
"Pa");
503 auto& eosNode = speciesNode[
"equation-of-state"].getMapWhere(
504 "model",
"Redlich-Kwong",
true);
507 if (coeffs.second == 0) {
508 bin_a[name2].setQuantity(coeffs.first,
"Pa*m^6/kmol^2*K^0.5");
510 vector<AnyValue> C(2);
511 C[0].setQuantity(coeffs.first,
"Pa*m^6/kmol^2*K^0.5");
512 C[1].setQuantity(coeffs.second,
"Pa*m^6/kmol^2/K^0.5");
513 bin_a[name2] = std::move(C);
516 eosNode[
"binary-a"] = std::move(bin_a);
525 double molarV = mmw / rho;
528 double dadt = da_dt();
530 double sqT = sqrt(T);
541 double molarV = mmw / rho;
544 double dadt = da_dt();
546 double sqT = sqrt(T);
557 double pres = std::max(
psatEst(TKelvin), presGuess);
559 bool foundLiq =
false;
561 while (m < 100 && !foundLiq) {
562 int nsol =
solveCubic(TKelvin, pres, atmp, btmp, Vroot);
563 if (nsol == 1 || nsol == 2) {
590 if (rhoguess == -1.0) {
591 if (phaseRequested >= FLUID_LIQUID_0) {
593 rhoguess = mmw / lqvol;
601 double volguess = mmw / rhoguess;
604 double molarVolLast = Vroot_[0];
606 if (phaseRequested >= FLUID_LIQUID_0) {
607 molarVolLast = Vroot_[0];
608 }
else if (phaseRequested == FLUID_GAS || phaseRequested == FLUID_SUPERCRIT) {
609 molarVolLast = Vroot_[2];
611 if (volguess > Vroot_[1]) {
612 molarVolLast = Vroot_[2];
614 molarVolLast = Vroot_[0];
617 }
else if (NSolns_ == 1) {
618 if (phaseRequested == FLUID_GAS || phaseRequested == FLUID_SUPERCRIT || phaseRequested == FLUID_UNDEFINED) {
619 molarVolLast = Vroot_[0];
623 }
else if (NSolns_ == -1) {
624 if (phaseRequested >= FLUID_LIQUID_0 || phaseRequested == FLUID_UNDEFINED || phaseRequested == FLUID_SUPERCRIT) {
625 molarVolLast = Vroot_[0];
626 }
else if (TKelvin > tcrit) {
627 molarVolLast = Vroot_[0];
632 molarVolLast = Vroot_[0];
635 return mmw / molarVolLast;
647 auto resid = [
this, T](
double v) {
652 boost::uintmax_t maxiter = 100;
653 auto [lower, upper] = bmt::toms748_solve(
654 resid, Vroot[0], Vroot[1], bmt::eps_tolerance<double>(48), maxiter);
657 return mmw / (0.5 * (lower + upper));
669 auto resid = [
this, T](
double v) {
674 boost::uintmax_t maxiter = 100;
675 auto [lower, upper] = bmt::toms748_solve(
676 resid, Vroot[1], Vroot[2], bmt::eps_tolerance<double>(48), maxiter);
679 return mmw / (0.5 * (lower + upper));
684 double sqt = sqrt(TKelvin);
690 double dpdv = (-
GasConstant * TKelvin / (vmb * vmb)
720 double sqt = sqrt(TKelvin);
723 double dadt = da_dt();
732 for (
size_t i = 0; i <
m_kk; i++) {
733 for (
size_t j = 0; j <
m_kk; j++) {
734 size_t counter = i *
m_kk + j;
735 a_vec_Curr_[counter] = a_coeff_vec(0,counter) + a_coeff_vec(1,counter) * temp;
742 for (
size_t i = 0; i <
m_kk; i++) {
744 for (
size_t j = 0; j <
m_kk; j++) {
750 fmt::memory_buffer b;
751 for (
size_t k = 0; k <
m_kk; k++) {
752 if (isnan(b_vec_Curr_[k])) {
760 throw CanteraError(
"RedlichKwongMFTP::updateMixingExpressions",
761 "Missing Redlich-Kwong coefficients for species: {}", to_string(b));
770 for (
size_t i = 0; i <
m_kk; i++) {
772 for (
size_t j = 0; j <
m_kk; j++) {
773 size_t counter = i *
m_kk + j;
774 double a_vec_Curr = a_coeff_vec(0,counter) + a_coeff_vec(1,counter) * temp;
779 for (
size_t i = 0; i <
m_kk; i++) {
781 for (
size_t j = 0; j <
m_kk; j++) {
782 size_t counter = i *
m_kk + j;
783 double a_vec_Curr = a_coeff_vec(0,counter);
790double RedlichKwongMFTP::da_dt()
const
794 for (
size_t i = 0; i <
m_kk; i++) {
795 for (
size_t j = 0; j <
m_kk; j++) {
796 size_t counter = i *
m_kk + j;
804void RedlichKwongMFTP::calcCriticalConditions(
double& pc,
double& tc,
double& vc)
const
808 for (
size_t i = 0; i <
m_kk; i++) {
809 for (
size_t j = 0; j <
m_kk; j++) {
810 size_t counter = i +
m_kk * j;
834 double sqrttc, f, dfdt, deltatc;
840 for (
int j = 0; j < 10; j++) {
844 deltatc = - f / dfdt;
848 throw CanteraError(
"RedlichKwongMFTP::calcCriticalConditions",
863 double sqt = sqrt(T);
864 double cn = - (
GasConstant * T * b / pres - a/(pres * sqt) + b * b);
865 double dn = - (a * b / (pres * sqt));
869 double tc = pow(tmp, pp);
873 int nSolnValues =
MixtureFugacityTP::solveCubic(T, pres, a, b, a, Vroot, 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.
A wrapper for a variable whose type is determined at runtime.
void getRow(size_t n, double *const rw)
Get the nth row and return it in a vector.
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 getEntropy_R_ref(double *er) const override
Returns the vector of nondimensional entropies of the reference state at the current temperature of t...
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 ...
void checkSpeciesIndex(size_t k) const
Check that the specified species index is in range.
size_t nSpecies() const
Returns the number of species in the phase.
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)
string speciesName(size_t k) const
Name of the species with index k.
map< string, shared_ptr< Species > > m_species
Map of Species objects.
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.
virtual double density() const
Density (kg/m^3).
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 dpdT_
The derivative of the pressure wrt the temperature.
void calculateAB(double temp, double &aCalc, double &bCalc) const
Calculate the a and the b parameters given the temperature.
double thermalExpansionCoeff() const override
Return the volumetric thermal expansion coefficient. Units: 1/K.
void setBinaryCoeffs(const string &species_i, const string &species_j, double a0, double a1)
Set values for the interaction parameter between two species.
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.
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).
int m_formTempParam
Form of the temperature parameterization.
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...
RedlichKwongMFTP(const string &infile="", const string &id="")
Construct a RedlichKwongMFTP object from an input file.
static const double omega_b
Omega constant for b.
vector< double > m_pp
Temporary storage - length = m_kk.
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...
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.
void pressureDerivatives() const
Calculate dpdV and dpdT 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.
static const double omega_a
Omega constant for a -> value of a in terms of critical properties.
double liquidVolEst(double TKelvin, double &pres) const override
Estimate for the molar volume of the liquid.
double dpdVCalc(double TKelvin, double molarVol, double &presCalc) const override
Calculate the pressure and the pressure derivative given the temperature and the molar volume.
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 and b coefficients.
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 a and b parameters.
double cp_mole() const override
Molar heat capacity at constant pressure. Units: J/kmol/K.
double m_a_current
Value of a in the equation of state.
double standardConcentration(size_t k=0) const override
Returns the standard concentration , which is used to normalize the generalized concentration.
double dpdV_
The derivative of the pressure wrt the volume.
bool addSpecies(shared_ptr< Species > spec) override
Add a Species to this Phase.
double m_b_current
Value of b in the equation of state.
void getActivityCoefficients(double *ac) const override
Get the array of non-dimensional activity coefficients at the current solution temperature,...
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.
vector< double > dpdni_
Vector of derivatives of pressure wrt mole number.
void getPartialMolarEntropies(double *sbar) const override
Returns an array of partial molar entropies of the species in the solution.
map< string, map< string, pair< double, double > > > m_binaryParameters
Explicitly-specified binary interaction parameters.
void setSpeciesCoeffs(const string &species, double a0, double a1, double b)
Set the pure fluid interaction parameters for a species.
int solveCubic(double T, double pres, double a, double b, double Vroot[3]) const
Prepare variables and call the function to solve the cubic equation of state.
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 fmt_append(fmt::memory_buffer &b, const std::string &tmpl, Args... args)
Versions 6.2.0 and 6.2.1 of fmtlib do not include this define before they include windows....
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].
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...