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 "To be removed after Cantera 3.0. Use getChemPotentials instead.");
184 for (
size_t k = 0; k <
m_kk; k++) {
185 muRT[k] *= 1.0 /
RT();
192 for (
size_t k = 0; k <
m_kk; k++) {
194 mu[k] +=
RT()*(log(xx));
202 for (
size_t k = 0; k <
m_kk; k++) {
204 for (
size_t i = 0; i <
m_kk; i++) {
205 size_t counter = k +
m_kk*i;
212 for (
size_t k = 0; k <
m_kk; k++) {
213 mu[k] += (
RT() * log(pres/refP) -
RT() * log(pres * mv /
RT())
214 +
RT() * log(mv / vmb)
215 +
RT() * b_vec_Curr_[k] / vmb
232 double sqt = sqrt(TKelvin);
235 for (
size_t k = 0; k <
m_kk; k++) {
237 for (
size_t i = 0; i <
m_kk; i++) {
238 size_t counter = k +
m_kk*i;
242 for (
size_t k = 0; k <
m_kk; k++) {
243 dpdni_[k] =
RT()/vmb +
RT() * b_vec_Curr_[k] / (vmb * vmb) - 2.0 *
m_pp[k] / (sqt * mv * vpb)
244 +
m_a_current * b_vec_Curr_[k]/(sqt * mv * vpb * vpb);
246 double dadt = da_dt();
247 double fac = TKelvin * dadt - 3.0 *
m_a_current / 2.0;
249 for (
size_t k = 0; k <
m_kk; k++) {
251 for (
size_t i = 0; i <
m_kk; i++) {
252 size_t counter = k +
m_kk*i;
259 for (
size_t k = 0; k <
m_kk; k++) {
262 + b_vec_Curr_[k] / vpb / (
m_b_current * sqt) * fac);
263 hbar[k] = hbar[k] + hE_v;
264 hbar[k] -= fac2 *
dpdni_[k];
273 double sqt = sqrt(TKelvin);
277 for (
size_t k = 0; k <
m_kk; k++) {
281 for (
size_t k = 0; k <
m_kk; k++) {
283 for (
size_t i = 0; i <
m_kk; i++) {
284 size_t counter = k +
m_kk*i;
288 for (
size_t k = 0; k <
m_kk; k++) {
290 for (
size_t i = 0; i <
m_kk; i++) {
291 size_t counter = k +
m_kk*i;
296 double dadt = da_dt();
300 for (
size_t k = 0; k <
m_kk; k++) {
308 - 1.0 / (
m_b_current * sqt) * b_vec_Curr_[k] / vpb * fac
314 for (
size_t k = 0; k <
m_kk; k++) {
315 sbar[k] -= -m_partialMolarVolumes[k] *
dpdT_;
325 for (
size_t k = 0; k <
nSpecies(); k++) {
326 ubar[k] -= p * m_partialMolarVolumes[k];
332 for (
size_t k = 0; k <
m_kk; k++) {
334 for (
size_t i = 0; i <
m_kk; i++) {
335 size_t counter = k +
m_kk*i;
339 for (
size_t k = 0; k <
m_kk; k++) {
341 for (
size_t i = 0; i <
m_kk; i++) {
342 size_t counter = k +
m_kk*i;
351 for (
size_t k = 0; k <
m_kk; k++) {
354 - 2.0 *
m_pp[k] / (sqt * vpb)
359 vbar[k] = num / denom;
367 a_vec_Curr_.resize(
m_kk *
m_kk, 0.0);
371 b_vec_Curr_.push_back(NAN);
376 m_partialMolarVolumes.push_back(0.0);
386 std::unordered_map<string, AnyMap*> dbSpecies;
391 if (!isnan(a_coeff_vec(0, k +
m_kk * k))) {
394 bool foundCoeffs =
false;
396 if (data.hasKey(
"equation-of-state") &&
397 data[
"equation-of-state"].hasMapWhere(
"model",
"Redlich-Kwong"))
401 auto eos = data[
"equation-of-state"].getMapWhere(
402 "model",
"Redlich-Kwong");
404 if (eos.hasKey(
"a") && eos.hasKey(
"b")) {
405 double a0 = 0, a1 = 0;
406 if (eos[
"a"].isScalar()) {
407 a0 = eos.convert(
"a",
"Pa*m^6/kmol^2*K^0.5");
409 auto avec = eos[
"a"].asVector<
AnyValue>(2);
410 a0 = eos.units().convert(avec[0],
"Pa*m^6/kmol^2*K^0.5");
411 a1 = eos.units().convert(avec[1],
"Pa*m^6/kmol^2/K^0.5");
413 double b = eos.convert(
"b",
"m^3/kmol");
419 if (eos.hasKey(
"binary-a")) {
422 for (
auto& [name2, coeff] : binary_a) {
423 double a0 = 0, a1 = 0;
424 if (coeff.isScalar()) {
425 a0 = units.
convert(coeff,
"Pa*m^6/kmol^2*K^0.5");
427 auto avec = coeff.asVector<
AnyValue>(2);
428 a0 = units.convert(avec[0],
"Pa*m^6/kmol^2*K^0.5");
429 a1 = units.convert(avec[1],
"Pa*m^6/kmol^2/K^0.5");
441 double Tc = NAN, Pc = NAN;
442 if (data.hasKey(
"critical-parameters")) {
445 auto& critProps = data[
"critical-parameters"].as<
AnyMap>();
446 Tc = critProps.
convert(
"critical-temperature",
"K");
447 Pc = critProps.convert(
"critical-pressure",
"Pa");
451 if (critPropsDb.
empty()) {
453 dbSpecies = critPropsDb[
"species"].asMap(
"name");
457 auto ucName = boost::algorithm::to_upper_copy(
name);
458 if (dbSpecies.count(ucName)) {
459 auto& spec = *dbSpecies.at(ucName);
460 auto& critProps = spec[
"critical-parameters"].as<
AnyMap>();
461 Tc = critProps.
convert(
"critical-temperature",
"K");
462 Pc = critProps.convert(
"critical-pressure",
"Pa");
475 "No critical property or Redlich-Kwong parameters found "
476 "for species {}.",
name);
482 AnyMap& speciesNode)
const
488 auto& eosNode = speciesNode[
"equation-of-state"].getMapWhere(
489 "model",
"Redlich-Kwong",
true);
491 size_t counter = k +
m_kk * k;
492 if (a_coeff_vec(1, counter) != 0.0) {
493 vector<AnyValue> coeffs(2);
494 coeffs[0].setQuantity(a_coeff_vec(0, counter),
"Pa*m^6/kmol^2*K^0.5");
495 coeffs[1].setQuantity(a_coeff_vec(1, counter),
"Pa*m^6/kmol^2/K^0.5");
496 eosNode[
"a"] = std::move(coeffs);
498 eosNode[
"a"].setQuantity(a_coeff_vec(0, counter),
499 "Pa*m^6/kmol^2*K^0.5");
501 eosNode[
"b"].setQuantity(b_vec_Curr_[k],
"m^3/kmol");
503 auto& critProps = speciesNode[
"critical-parameters"];
504 double a = a_coeff_vec(0, k +
m_kk * k);
505 double b = b_vec_Curr_[k];
508 critProps[
"critical-temperature"].setQuantity(Tc,
"K");
509 critProps[
"critical-pressure"].setQuantity(Pc,
"Pa");
513 auto& eosNode = speciesNode[
"equation-of-state"].getMapWhere(
514 "model",
"Redlich-Kwong",
true);
517 if (coeffs.second == 0) {
518 bin_a[name2].setQuantity(coeffs.first,
"Pa*m^6/kmol^2*K^0.5");
520 vector<AnyValue> C(2);
521 C[0].setQuantity(coeffs.first,
"Pa*m^6/kmol^2*K^0.5");
522 C[1].setQuantity(coeffs.second,
"Pa*m^6/kmol^2/K^0.5");
523 bin_a[name2] = std::move(C);
526 eosNode[
"binary-a"] = std::move(bin_a);
535 double molarV = mmw / rho;
538 double dadt = da_dt();
540 double sqT = sqrt(T);
551 double molarV = mmw / rho;
554 double dadt = da_dt();
556 double sqT = sqrt(T);
567 double pres = std::max(
psatEst(TKelvin), presGuess);
569 bool foundLiq =
false;
571 while (m < 100 && !foundLiq) {
572 int nsol =
solveCubic(TKelvin, pres, atmp, btmp, Vroot);
573 if (nsol == 1 || nsol == 2) {
600 if (rhoguess == -1.0) {
601 if (phaseRequested >= FLUID_LIQUID_0) {
603 rhoguess = mmw / lqvol;
611 double volguess = mmw / rhoguess;
614 double molarVolLast = Vroot_[0];
616 if (phaseRequested >= FLUID_LIQUID_0) {
617 molarVolLast = Vroot_[0];
618 }
else if (phaseRequested == FLUID_GAS || phaseRequested == FLUID_SUPERCRIT) {
619 molarVolLast = Vroot_[2];
621 if (volguess > Vroot_[1]) {
622 molarVolLast = Vroot_[2];
624 molarVolLast = Vroot_[0];
627 }
else if (NSolns_ == 1) {
628 if (phaseRequested == FLUID_GAS || phaseRequested == FLUID_SUPERCRIT || phaseRequested == FLUID_UNDEFINED) {
629 molarVolLast = Vroot_[0];
633 }
else if (NSolns_ == -1) {
634 if (phaseRequested >= FLUID_LIQUID_0 || phaseRequested == FLUID_UNDEFINED || phaseRequested == FLUID_SUPERCRIT) {
635 molarVolLast = Vroot_[0];
636 }
else if (TKelvin > tcrit) {
637 molarVolLast = Vroot_[0];
642 molarVolLast = Vroot_[0];
645 return mmw / molarVolLast;
657 auto resid = [
this, T](
double v) {
662 boost::uintmax_t maxiter = 100;
663 auto [lower, upper] = bmt::toms748_solve(
664 resid, Vroot[0], Vroot[1], bmt::eps_tolerance<double>(48), maxiter);
667 return mmw / (0.5 * (lower + upper));
679 auto resid = [
this, T](
double v) {
684 boost::uintmax_t maxiter = 100;
685 auto [lower, upper] = bmt::toms748_solve(
686 resid, Vroot[1], Vroot[2], bmt::eps_tolerance<double>(48), maxiter);
689 return mmw / (0.5 * (lower + upper));
694 double sqt = sqrt(TKelvin);
700 double dpdv = (-
GasConstant * TKelvin / (vmb * vmb)
730 double sqt = sqrt(TKelvin);
733 double dadt = da_dt();
742 for (
size_t i = 0; i <
m_kk; i++) {
743 for (
size_t j = 0; j <
m_kk; j++) {
744 size_t counter = i *
m_kk + j;
745 a_vec_Curr_[counter] = a_coeff_vec(0,counter) + a_coeff_vec(1,counter) * temp;
752 for (
size_t i = 0; i <
m_kk; i++) {
754 for (
size_t j = 0; j <
m_kk; j++) {
760 fmt::memory_buffer b;
761 for (
size_t k = 0; k <
m_kk; k++) {
762 if (isnan(b_vec_Curr_[k])) {
770 throw CanteraError(
"RedlichKwongMFTP::updateMixingExpressions",
771 "Missing Redlich-Kwong coefficients for species: {}", to_string(b));
780 for (
size_t i = 0; i <
m_kk; i++) {
782 for (
size_t j = 0; j <
m_kk; j++) {
783 size_t counter = i *
m_kk + j;
784 double a_vec_Curr = a_coeff_vec(0,counter) + a_coeff_vec(1,counter) * temp;
789 for (
size_t i = 0; i <
m_kk; i++) {
791 for (
size_t j = 0; j <
m_kk; j++) {
792 size_t counter = i *
m_kk + j;
793 double a_vec_Curr = a_coeff_vec(0,counter);
800double RedlichKwongMFTP::da_dt()
const
804 for (
size_t i = 0; i <
m_kk; i++) {
805 for (
size_t j = 0; j <
m_kk; j++) {
806 size_t counter = i *
m_kk + j;
814void RedlichKwongMFTP::calcCriticalConditions(
double& pc,
double& tc,
double& vc)
const
818 for (
size_t i = 0; i <
m_kk; i++) {
819 for (
size_t j = 0; j <
m_kk; j++) {
820 size_t counter = i +
m_kk * j;
844 double sqrttc, f, dfdt, deltatc;
850 for (
int j = 0; j < 10; j++) {
854 deltatc = - f / dfdt;
858 throw CanteraError(
"RedlichKwongMFTP::calcCriticalConditions",
873 double sqt = sqrt(T);
874 double cn = - (
GasConstant * T * b / pres - a/(pres * sqt) + b * b);
875 double dn = - (a * b / (pres * sqt));
879 double tc = pow(tmp, pp);
883 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.
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 getChemPotentials_RT(double *mu) const override
Get the array of non-dimensional species chemical potentials.
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, 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.
void warn_deprecated(const string &source, const AnyBase &node, const string &message)
A deprecation warning for syntax in an input file.
Contains declarations for string manipulation functions within Cantera.
Various templated functions that carry out common vector and polynomial operations (see Templated Arr...