IdealMolalSoln.cpp

Go to the documentation of this file.

/**
 *  @file IdealMolalSoln.cpp
 *   ThermoPhase object for the ideal molal equation of
 * state (see @ref thermoprops
 * and class @link Cantera::IdealMolalSoln IdealMolalSoln@endlink).
 *
 * Definition file for a derived class of ThermoPhase that handles variable
 * pressure standard state methods for calculating thermodynamic properties that
 * are further based upon activities on the molality scale. The Ideal molal
 * solution assumes that all molality-based activity coefficients are equal to
 * one. This turns out, actually, to be highly nonlinear when the solvent
 * densities get low.
 */
 
// 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.
 
#include "cantera/thermo/IdealMolalSoln.h"
#include "cantera/thermo/ThermoFactory.h"
#include "cantera/thermo/PDSS.h"
#include "cantera/base/stringUtils.h"
 
#include <iostream>
 
namespace {
double X_o_cutoff_default = 0.20;
double gamma_o_min_default = 0.00001;
double gamma_k_min_default = 10.0;
double slopefCut_default = 0.6;
double slopegCut_default = 0.0;
double cCut_default = .05;
}
 
namespace Cantera
{
 
IdealMolalSoln::IdealMolalSoln(const string& inputFile, const string& id_) :
    IMS_X_o_cutoff_(X_o_cutoff_default),
    IMS_gamma_o_min_(gamma_o_min_default),
    IMS_gamma_k_min_(gamma_k_min_default),
    IMS_slopefCut_(slopefCut_default),
    IMS_slopegCut_(slopegCut_default),
    IMS_cCut_(cCut_default)
{
    initThermoFile(inputFile, id_);
}
 
double IdealMolalSoln::intEnergy_mole() const
{
    getPartialMolarIntEnergies(m_workS.data());
    return mean_X(m_workS);
}
 
// ------- Mechanical Equation of State Properties ------------------------
 
double IdealMolalSoln::isothermalCompressibility() const
{
    return 0.0;
}
 
double IdealMolalSoln::thermalExpansionCoeff() const
{
    return 0.0;
}
 
// ------- Activities and Activity Concentrations
 
Units IdealMolalSoln::standardConcentrationUnits() const
{
    if (m_formGC == 0) {
        return Units(1.0); // dimensionless
    } else {
        // kmol/m^3 for bulk phases
        return Units(1.0, 0, -static_cast<double>(nDim()), 0, 0, 0, 1);
    }
}
 
void IdealMolalSoln::getActivityConcentrations(double* c) const
{
    if (m_formGC != 1) {
        double c_solvent = standardConcentration();
        getActivities(c);
        for (size_t k = 0; k < m_kk; k++) {
            c[k] *= c_solvent;
        }
    } else {
        getActivities(c);
        for (size_t k = 0; k < m_kk; k++) {
            double c0 = standardConcentration(k);
            c[k] *= c0;
        }
    }
}
 
double IdealMolalSoln::standardConcentration(size_t k) const
{
    switch (m_formGC) {
    case 0:
        return 1.0;
    case 1:
        return 1.0 / m_speciesMolarVolume[k];
    case 2:
        return 1.0 / m_speciesMolarVolume[0];
    default:
        throw CanteraError("IdealMolalSoln::standardConcentration",
                       "m_formGC is set to an incorrect value. \
                       Allowed values are 0, 1, and 2");
    }
}
 
void IdealMolalSoln::getActivities(double* ac) const
{
    _updateStandardStateThermo();
 
    // Update the molality array, m_molalities(). This requires an update due to
    // mole fractions
    if (IMS_typeCutoff_ == 0) {
        calcMolalities();
        for (size_t k = 0; k < m_kk; k++) {
            ac[k] = m_molalities[k];
        }
        double xmolSolvent = moleFraction(0);
        // Limit the activity coefficient to be finite as the solvent mole
        // fraction goes to zero.
        xmolSolvent = std::max(m_xmolSolventMIN, xmolSolvent);
        ac[0] = exp((xmolSolvent - 1.0)/xmolSolvent);
    } else {
 
        s_updateIMS_lnMolalityActCoeff();
 
        // Now calculate the array of activities.
        for (size_t k = 1; k < m_kk; k++) {
            ac[k] = m_molalities[k] * exp(IMS_lnActCoeffMolal_[k]);
        }
        double xmolSolvent = moleFraction(0);
        ac[0] = exp(IMS_lnActCoeffMolal_[0]) * xmolSolvent;
    }
}
 
void IdealMolalSoln::getMolalityActivityCoefficients(double* acMolality) const
{
    if (IMS_typeCutoff_ == 0) {
        for (size_t k = 0; k < m_kk; k++) {
            acMolality[k] = 1.0;
        }
        double xmolSolvent = moleFraction(0);
        // Limit the activity coefficient to be finite as the solvent mole
        // fraction goes to zero.
        xmolSolvent = std::max(m_xmolSolventMIN, xmolSolvent);
        acMolality[0] = exp((xmolSolvent - 1.0)/xmolSolvent) / xmolSolvent;
    } else {
        s_updateIMS_lnMolalityActCoeff();
        std::copy(IMS_lnActCoeffMolal_.begin(), IMS_lnActCoeffMolal_.end(), acMolality);
        for (size_t k = 0; k < m_kk; k++) {
            acMolality[k] = exp(acMolality[k]);
        }
    }
}
 
// ------ Partial Molar Properties of the Solution -----------------
 
void IdealMolalSoln::getChemPotentials(double* mu) const
{
    // First get the standard chemical potentials. This requires updates of
    // standard state as a function of T and P These are defined at unit
    // molality.
    getStandardChemPotentials(mu);
 
    // Update the molality array, m_molalities(). This requires an update due to
    // mole fractions
    calcMolalities();
 
    // get the solvent mole fraction
    double xmolSolvent = moleFraction(0);
 
    if (IMS_typeCutoff_ == 0 || xmolSolvent > 3.* IMS_X_o_cutoff_/2.0) {
        for (size_t k = 1; k < m_kk; k++) {
            double xx = std::max(m_molalities[k], SmallNumber);
            mu[k] += RT() * log(xx);
        }
 
        // Do the solvent
        //  -> see my notes
        double xx = std::max(xmolSolvent, SmallNumber);
        mu[0] += (RT() * (xmolSolvent - 1.0) / xx);
    } else {
        // Update the activity coefficients. This also updates the internal
        // molality array.
        s_updateIMS_lnMolalityActCoeff();
 
        for (size_t k = 1; k < m_kk; k++) {
            double xx = std::max(m_molalities[k], SmallNumber);
            mu[k] += RT() * (log(xx) + IMS_lnActCoeffMolal_[k]);
        }
        double xx = std::max(xmolSolvent, SmallNumber);
        mu[0] += RT() * (log(xx) + IMS_lnActCoeffMolal_[0]);
    }
}
 
void IdealMolalSoln::getPartialMolarEnthalpies(double* hbar) const
{
    getEnthalpy_RT(hbar);
    for (size_t k = 0; k < m_kk; k++) {
        hbar[k] *= RT();
    }
}
 
void IdealMolalSoln::getPartialMolarIntEnergies(double* ubar) const
{
    getIntEnergy_RT(ubar);
    for (size_t k = 0; k < m_kk; k++) {
        ubar[k] *= RT();
    }
}
 
void IdealMolalSoln::getPartialMolarEntropies(double* sbar) const
{
    getEntropy_R(sbar);
    calcMolalities();
    if (IMS_typeCutoff_ == 0) {
        for (size_t k = 1; k < m_kk; k++) {
            double mm = std::max(SmallNumber, m_molalities[k]);
            sbar[k] -= GasConstant * log(mm);
        }
        double xmolSolvent = moleFraction(0);
        sbar[0] -= (GasConstant * (xmolSolvent - 1.0) / xmolSolvent);
    } else {
        // Update the activity coefficients, This also update the internally
        // stored molalities.
        s_updateIMS_lnMolalityActCoeff();
 
        // First we will add in the obvious dependence on the T term out front
        // of the log activity term
        double mm;
        for (size_t k = 1; k < m_kk; k++) {
            mm = std::max(SmallNumber, m_molalities[k]);
            sbar[k] -= GasConstant * (log(mm) + IMS_lnActCoeffMolal_[k]);
        }
        double xmolSolvent = moleFraction(0);
        mm = std::max(SmallNumber, xmolSolvent);
        sbar[0] -= GasConstant *(log(mm) + IMS_lnActCoeffMolal_[0]);
    }
}
 
void IdealMolalSoln::getPartialMolarVolumes(double* vbar) const
{
    getStandardVolumes(vbar);
}
 
void IdealMolalSoln::getPartialMolarCp(double* cpbar) const
{
    // Get the nondimensional Gibbs standard state of the species at the T and P
    // of the solution.
    getCp_R(cpbar);
    for (size_t k = 0; k < m_kk; k++) {
        cpbar[k] *= GasConstant;
    }
}
 
// -------------- Utilities -------------------------------
 
bool IdealMolalSoln::addSpecies(shared_ptr<Species> spec)
{
    bool added = MolalityVPSSTP::addSpecies(spec);
    if (added) {
        m_speciesMolarVolume.push_back(0.0);
        IMS_lnActCoeffMolal_.push_back(0.0);
    }
    return added;
}
 
void IdealMolalSoln::initThermo()
{
    MolalityVPSSTP::initThermo();
 
    if (m_input.hasKey("standard-concentration-basis")) {
        setStandardConcentrationModel(m_input["standard-concentration-basis"].asString());
    }
    if (m_input.hasKey("cutoff")) {
        auto& cutoff = m_input["cutoff"].as<AnyMap>();
        setCutoffModel(cutoff.getString("model", "none"));
        IMS_gamma_o_min_ = cutoff.getDouble("gamma_o", gamma_o_min_default);
        IMS_gamma_k_min_ = cutoff.getDouble("gamma_k", gamma_k_min_default);
        IMS_X_o_cutoff_ = cutoff.getDouble("X_o", X_o_cutoff_default);
        IMS_cCut_ = cutoff.getDouble("c_0", cCut_default);
        IMS_slopefCut_ = cutoff.getDouble("slope_f", slopefCut_default);
        IMS_slopegCut_ = cutoff.getDouble("slope_g", slopegCut_default);
    }
 
    for (size_t k = 0; k < nSpecies(); k++) {
        m_speciesMolarVolume[k] = providePDSS(k)->molarVolume();
    }
    if (IMS_typeCutoff_ == 2) {
        calcIMSCutoffParams_();
    }
    setMoleFSolventMin(1.0E-5);
}
 
void IdealMolalSoln::getParameters(AnyMap& phaseNode) const
{
    MolalityVPSSTP::getParameters(phaseNode);
 
    // "solvent-molar-volume" (m_formGC == 2) is the default, and can be omitted
    if (m_formGC == 0) {
        phaseNode["standard-concentration-basis"] = "unity";
    } else if (m_formGC == 1) {
        phaseNode["standard-concentration-basis"] = "species-molar-volume";
    }
 
    AnyMap cutoff;
    if (IMS_typeCutoff_ == 1) {
        cutoff["model"] = "poly";
    } else if (IMS_typeCutoff_ == 2) {
        cutoff["model"] = "polyexp";
    }
 
    if (IMS_gamma_o_min_ != gamma_o_min_default) {
        cutoff["gamma_o"] = IMS_gamma_o_min_;
    }
    if (IMS_gamma_k_min_ != gamma_k_min_default) {
        cutoff["gamma_k"] = IMS_gamma_k_min_;
    }
    if (IMS_X_o_cutoff_ != X_o_cutoff_default) {
        cutoff["X_o"] = IMS_X_o_cutoff_;
    }
    if (IMS_cCut_ != cCut_default) {
        cutoff["c_0"] = IMS_cCut_;
    }
    if (IMS_slopefCut_ != slopefCut_default) {
        cutoff["slope_f"] = IMS_slopefCut_;
    }
    if (IMS_slopegCut_ != slopegCut_default) {
        cutoff["slope_g"] = IMS_slopegCut_;
    }
 
    if (cutoff.size()) {
        phaseNode["cutoff"] = std::move(cutoff);
    }
}
 
void IdealMolalSoln::setStandardConcentrationModel(const string& model)
{
    if (caseInsensitiveEquals(model, "unity")) {
        m_formGC = 0;
    } else if (caseInsensitiveEquals(model, "species-molar-volume")
               || caseInsensitiveEquals(model, "molar_volume")) {
        m_formGC = 1;
    } else if (caseInsensitiveEquals(model, "solvent-molar-volume")
               || caseInsensitiveEquals(model, "solvent_volume")) {
        m_formGC = 2;
    } else {
        throw CanteraError("IdealMolalSoln::setStandardConcentrationModel",
                           "Unknown standard concentration model '{}'", model);
    }
}
 
void IdealMolalSoln::setCutoffModel(const string& model)
{
    if (caseInsensitiveEquals(model, "none")) {
        IMS_typeCutoff_ = 0;
    } else if (caseInsensitiveEquals(model, "poly")) {
        IMS_typeCutoff_ = 1;
    } else if (caseInsensitiveEquals(model, "polyexp")) {
        IMS_typeCutoff_ = 2;
    } else {
        throw CanteraError("IdealMolalSoln::setCutoffModel",
                           "Unknown cutoff model '{}'", model);
    }
}
 
// ------------ Private and Restricted Functions ------------------
 
void IdealMolalSoln::s_updateIMS_lnMolalityActCoeff() const
{
    // Calculate the molalities. Currently, the molalities may not be current
    // with respect to the contents of the State objects' data.
    calcMolalities();
 
    double xmolSolvent = moleFraction(0);
    double xx = std::max(m_xmolSolventMIN, xmolSolvent);
 
    if (IMS_typeCutoff_ == 0) {
        for (size_t k = 1; k < m_kk; k++) {
            IMS_lnActCoeffMolal_[k]= 0.0;
        }
        IMS_lnActCoeffMolal_[0] = - log(xx) + (xx - 1.0)/xx;
        return;
    } else if (IMS_typeCutoff_ == 1) {
        if (xmolSolvent > 3.0 * IMS_X_o_cutoff_/2.0) {
            for (size_t k = 1; k < m_kk; k++) {
                IMS_lnActCoeffMolal_[k]= 0.0;
            }
            IMS_lnActCoeffMolal_[0] = - log(xx) + (xx - 1.0)/xx;
            return;
        } else if (xmolSolvent < IMS_X_o_cutoff_/2.0) {
            double tmp = log(xx * IMS_gamma_k_min_);
            for (size_t k = 1; k < m_kk; k++) {
                IMS_lnActCoeffMolal_[k]= tmp;
            }
            IMS_lnActCoeffMolal_[0] = log(IMS_gamma_o_min_);
            return;
        } else {
            // If we are in the middle region, calculate the connecting polynomials
            double xminus = xmolSolvent - IMS_X_o_cutoff_/2.0;
            double xminus2 = xminus * xminus;
            double xminus3 = xminus2 * xminus;
            double x_o_cut2 = IMS_X_o_cutoff_ * IMS_X_o_cutoff_;
            double x_o_cut3 = x_o_cut2 * IMS_X_o_cutoff_;
 
            double h2 = 3.5 * xminus2 / IMS_X_o_cutoff_ - 2.0 * xminus3 / x_o_cut2;
            double h2_prime = 7.0 * xminus / IMS_X_o_cutoff_ - 6.0 * xminus2 / x_o_cut2;
 
            double h1 = (1.0 - 3.0 * xminus2 / x_o_cut2 + 2.0 * xminus3/ x_o_cut3);
            double h1_prime = (- 6.0 * xminus / x_o_cut2 + 6.0 * xminus2/ x_o_cut3);
 
            double h1_g = h1 / IMS_gamma_o_min_;
            double h1_g_prime = h1_prime / IMS_gamma_o_min_;
 
            double alpha = 1.0 / (exp(1.0) * IMS_gamma_k_min_);
            double h1_f = h1 * alpha;
            double h1_f_prime = h1_prime * alpha;
 
            double f = h2 + h1_f;
            double f_prime = h2_prime + h1_f_prime;
 
            double g = h2 + h1_g;
            double g_prime = h2_prime + h1_g_prime;
 
            double tmp = (xmolSolvent/ g * g_prime + (1.0-xmolSolvent) / f * f_prime);
            double lngammak = -1.0 - log(f) + tmp * xmolSolvent;
            double lngammao =-log(g) - tmp * (1.0-xmolSolvent);
 
            tmp = log(xmolSolvent) + lngammak;
            for (size_t k = 1; k < m_kk; k++) {
                IMS_lnActCoeffMolal_[k]= tmp;
            }
            IMS_lnActCoeffMolal_[0] = lngammao;
        }
    } else if (IMS_typeCutoff_ == 2) {
        // Exponentials - trial 2
        if (xmolSolvent > IMS_X_o_cutoff_) {
            for (size_t k = 1; k < m_kk; k++) {
                IMS_lnActCoeffMolal_[k]= 0.0;
            }
            IMS_lnActCoeffMolal_[0] = - log(xx) + (xx - 1.0)/xx;
            return;
        } else {
            double xoverc = xmolSolvent/IMS_cCut_;
            double eterm = std::exp(-xoverc);
 
            double fptmp = IMS_bfCut_ - IMS_afCut_ / IMS_cCut_ - IMS_bfCut_*xoverc
                           + 2.0*IMS_dfCut_*xmolSolvent - IMS_dfCut_*xmolSolvent*xoverc;
            double f_prime = 1.0 + eterm*fptmp;
            double f = xmolSolvent + IMS_efCut_ + eterm * (IMS_afCut_ + xmolSolvent * (IMS_bfCut_ + IMS_dfCut_*xmolSolvent));
 
            double gptmp = IMS_bgCut_ - IMS_agCut_ / IMS_cCut_ - IMS_bgCut_*xoverc
                           + 2.0*IMS_dgCut_*xmolSolvent - IMS_dgCut_*xmolSolvent*xoverc;
            double g_prime = 1.0 + eterm*gptmp;
            double g = xmolSolvent + IMS_egCut_ + eterm * (IMS_agCut_ + xmolSolvent * (IMS_bgCut_ + IMS_dgCut_*xmolSolvent));
 
            double tmp = (xmolSolvent / g * g_prime + (1.0 - xmolSolvent) / f * f_prime);
            double lngammak = -1.0 - log(f) + tmp * xmolSolvent;
            double lngammao =-log(g) - tmp * (1.0-xmolSolvent);
 
            tmp = log(xx) + lngammak;
            for (size_t k = 1; k < m_kk; k++) {
                IMS_lnActCoeffMolal_[k]= tmp;
            }
            IMS_lnActCoeffMolal_[0] = lngammao;
        }
    }
}
 
void IdealMolalSoln::calcIMSCutoffParams_()
{
    IMS_afCut_ = 1.0 / (std::exp(1.0) * IMS_gamma_k_min_);
    IMS_efCut_ = 0.0;
    bool converged = false;
    for (int its = 0; its < 100 && !converged; its++) {
        double oldV = IMS_efCut_;
        IMS_afCut_ = 1.0 / (std::exp(1.0) * IMS_gamma_k_min_) - IMS_efCut_;
        IMS_bfCut_ = IMS_afCut_ / IMS_cCut_ + IMS_slopefCut_ - 1.0;
        IMS_dfCut_ = ((- IMS_afCut_/IMS_cCut_ + IMS_bfCut_ - IMS_bfCut_*IMS_X_o_cutoff_/IMS_cCut_)
                      /
                      (IMS_X_o_cutoff_*IMS_X_o_cutoff_/IMS_cCut_ - 2.0 * IMS_X_o_cutoff_));
        double tmp = IMS_afCut_ + IMS_X_o_cutoff_*(IMS_bfCut_ + IMS_dfCut_ * IMS_X_o_cutoff_);
        double eterm = std::exp(-IMS_X_o_cutoff_/IMS_cCut_);
        IMS_efCut_ = - eterm * (tmp);
        if (fabs(IMS_efCut_ - oldV) < 1.0E-14) {
            converged = true;
        }
    }
    if (!converged) {
        throw CanteraError("IdealMolalSoln::calcCutoffParams_",
                           "failed to converge on the f polynomial");
    }
    converged = false;
    double f_0 = IMS_afCut_ + IMS_efCut_;
    double f_prime_0 = 1.0 - IMS_afCut_ / IMS_cCut_ + IMS_bfCut_;
    IMS_egCut_ = 0.0;
    for (int its = 0; its < 100 && !converged; its++) {
        double oldV = IMS_egCut_;
        double lng_0 = -log(IMS_gamma_o_min_) - f_prime_0 / f_0;
        IMS_agCut_ = exp(lng_0) - IMS_egCut_;
        IMS_bgCut_ = IMS_agCut_ / IMS_cCut_ + IMS_slopegCut_ - 1.0;
        IMS_dgCut_ = ((- IMS_agCut_/IMS_cCut_ + IMS_bgCut_ - IMS_bgCut_*IMS_X_o_cutoff_/IMS_cCut_)
                      /
                      (IMS_X_o_cutoff_*IMS_X_o_cutoff_/IMS_cCut_ - 2.0 * IMS_X_o_cutoff_));
        double tmp = IMS_agCut_ + IMS_X_o_cutoff_*(IMS_bgCut_ + IMS_dgCut_ *IMS_X_o_cutoff_);
        double eterm = std::exp(-IMS_X_o_cutoff_/IMS_cCut_);
        IMS_egCut_ = - eterm * (tmp);
        if (fabs(IMS_egCut_ - oldV) < 1.0E-14) {
            converged = true;
        }
    }
    if (!converged) {
        throw CanteraError("IdealMolalSoln::calcCutoffParams_",
                           "failed to converge on the g polynomial");
    }
}
 
}