[Thermo] Allow instantiation of RedlichKwongMFTP without XML

This also adds the first test which instantiates a RedlichKwongMFTP object, and
removes some unused member variables and private methods from the class.

include/cantera/thermo/RedlichKwongMFTP.h

@@ -191,6 +191,40 @@ public:
     virtual void setToEquilState(const doublereal* lambda_RT);
     virtual void initThermoXML(XML_Node& phaseNode, const std::string& id);
 
+    //! Set the pure fluid interaction parameters for a species
+    /*!
+     *  The "a" parameter for species *i* in the Redlich-Kwong model is assumed
+     *  to be a linear function of temperature:
+     *  \f[ a = a_0 + a_1 T \f]
+     *
+     *  @param species   Name of the species
+     *  @param a0        constant term in the expression for the "a" parameter
+     *      of the specified species [Pa-m^6/kmol^2]
+     *  @param a1        temperature-proportional term in the expression for the
+     *      "a" parameter of the specified species [Pa-m^6/kmol^2/K]
+     *  @param b         "b" parameter in the Redlich-Kwong model [m^3/kmol]
+     */
+    void setSpeciesCoeffs(const std::string& species, double a0, double a1,
+                              double b);
+
+    //! Set values for the interaction parameter between two species
+    /*!
+     *  The "a" parameter for interactions between species *i* and *j* is
+     *  assumed by default to be computed as:
+     *  \f[ a_{ij} = \sqrt(a_{i,0} a_{j,0}) + \sqrt(a_{i,1} a_{j,1}) T \f]
+     *
+     *  This function overrides the defaults with the specified parameters:
+     *  \f[ a_{ij} = a_{ij,0} + a_{ij,1} T \f]
+     *
+     *  @param species_i   Name of one species
+     *  @param species_j   Name of the other species
+     *  @param a0          constant term in the "a" expression [Pa-m^6/kmol^2]
+     *  @param a1          temperature-proportional term in the "a" expression
+     *      [Pa-m^6/kmol^2/K]
+     */
+    void setBinaryCoeffs(const std::string& species_i,
+                         const std::string& species_j, double a0, double a1);
+
 private:
     //! Read the pure species RedlichKwong input parameters
     /*!
@@ -198,9 +232,6 @@ private:
      */
     void readXMLPureFluid(XML_Node& pureFluidParam);
 
-    //! Apply mixing rules for a coefficients
-    void applyStandardMixingRules();
-
     //! Read the cross species RedlichKwong input parameters
     /*!
      *  @param pureFluidParam   XML_Node for the cross fluid parameters
@@ -271,13 +302,6 @@ public:
                      doublereal Vroot[3]) const;
 
 protected:
-    //! boolean indicating whether standard mixing rules are applied
-    /*!
-     *  - 1 = Yes, there are standard cross terms in the a coefficient matrices.
-     *  - 0 = No, there are nonstandard cross terms in the a coefficient matrices.
-     */
-    int m_standardMixingRules;
-
     //! Form of the temperature parameterization
     /*!
      *  - 0 = There is no temperature parameterization of a or b
@@ -302,10 +326,6 @@ protected:
 
     Array2D a_coeff_vec;
 
-    vector_fp m_pc_Species;
-    vector_fp m_tc_Species;
-    vector_fp m_vc_Species;
-
     int NSolns_;
 
     doublereal Vroot_[3];

src/thermo/RedlichKwongMFTP.cpp

@@ -22,7 +22,6 @@ const doublereal RedlichKwongMFTP::omega_b = 8.66403499650E-02;
 const doublereal RedlichKwongMFTP::omega_vc = 3.33333333333333E-01;
 
 RedlichKwongMFTP::RedlichKwongMFTP() :
-    m_standardMixingRules(0),
     m_formTempParam(0),
     m_b_current(0.0),
     m_a_current(0.0),
@@ -34,7 +33,6 @@ RedlichKwongMFTP::RedlichKwongMFTP() :
 }
 
 RedlichKwongMFTP::RedlichKwongMFTP(const std::string& infile, const std::string& id_) :
-    m_standardMixingRules(0),
     m_formTempParam(0),
     m_b_current(0.0),
     m_a_current(0.0),
@@ -47,7 +45,6 @@ RedlichKwongMFTP::RedlichKwongMFTP(const std::string& infile, const std::string&
 }
 
 RedlichKwongMFTP::RedlichKwongMFTP(XML_Node& phaseRefRoot, const std::string& id_) :
-    m_standardMixingRules(0),
     m_formTempParam(0),
     m_b_current(0.0),
     m_a_current(0.0),
@@ -60,7 +57,6 @@ RedlichKwongMFTP::RedlichKwongMFTP(XML_Node& phaseRefRoot, const std::string& id
 }
 
 RedlichKwongMFTP::RedlichKwongMFTP(const RedlichKwongMFTP& b) :
-    m_standardMixingRules(0),
     m_formTempParam(0),
     m_b_current(0.0),
     m_a_current(0.0),
@@ -79,7 +75,6 @@ RedlichKwongMFTP& RedlichKwongMFTP::operator=(const RedlichKwongMFTP& b)
         MixtureFugacityTP::operator=(b);
 
         // However, we have to handle data that we own.
-        m_standardMixingRules = b.m_standardMixingRules;
         m_formTempParam = b.m_formTempParam;
         m_b_current = b.m_b_current;
         m_a_current = b.m_a_current;
@@ -87,9 +82,6 @@ RedlichKwongMFTP& RedlichKwongMFTP::operator=(const RedlichKwongMFTP& b)
         b_vec_Curr_ = b.b_vec_Curr_;
         a_coeff_vec = b.a_coeff_vec;
 
-        m_pc_Species = b.m_pc_Species;
-        m_tc_Species = b.m_tc_Species;
-        m_vc_Species = b.m_vc_Species;
         NSolns_ = b.NSolns_;
         Vroot_[0] = b.Vroot_[0];
         Vroot_[1] = b.Vroot_[1];
@@ -116,6 +108,63 @@ int RedlichKwongMFTP::eosType() const
     return cRedlichKwongMFTP;
 }
 
+void RedlichKwongMFTP::setSpeciesCoeffs(const std::string& species,
+                                        double a0, double a1, double b)
+{
+    size_t k = speciesIndex(species);
+    if (k == npos) {
+        throw CanteraError("RedlichKwongMFTP::setSpeciesCoeffs",
+            "Unknown species '{}'.", species);
+    }
+
+    if (a1 != 0.0) {
+        m_formTempParam = 1; // expression is temperature-dependent
+    }
+
+    size_t counter = k + m_kk * k;
+    a_coeff_vec(0, counter) = a0;
+    a_coeff_vec(1, counter) = a1;
+
+    // standard mixing rule for cross-species interaction term
+    for (size_t j = 0; j < m_kk; j++) {
+        if (k == j) {
+            continue;
+        }
+        double a0kj = sqrt(a_coeff_vec(0, j + m_kk * j) * a0);
+        double a1kj = sqrt(a_coeff_vec(1, j + m_kk * j) * a1);
+        if  (a_coeff_vec(0, j + m_kk * k) == 0) {
+            a_coeff_vec(0, j + m_kk * k) = a0kj;
+            a_coeff_vec(1, j + m_kk * k) = a1kj;
+            a_coeff_vec(0, k + m_kk * j) = a0kj;
+            a_coeff_vec(1, k + m_kk * j) = a1kj;
+        }
+    }
+    b_vec_Curr_[k] = b;
+}
+
+void RedlichKwongMFTP::setBinaryCoeffs(const std::string& species_i,
+        const std::string& species_j, double a0, double a1)
+{
+    size_t ki = speciesIndex(species_i);
+    if (ki == npos) {
+        throw CanteraError("RedlichKwongMFTP::setBinaryCoeffs",
+            "Unknown species '{}'.", species_i);
+    }
+    size_t kj = speciesIndex(species_j);
+    if (kj == npos) {
+        throw CanteraError("RedlichKwongMFTP::setBinaryCoeffs",
+            "Unknown species '{}'.", species_j);
+    }
+
+    if (a1 != 0.0) {
+        m_formTempParam = 1; // expression is temperature-dependent
+    }
+    size_t counter1 = ki + m_kk * kj;
+    size_t counter2 = kj + m_kk * ki;
+    a_coeff_vec(0, counter1) = a_coeff_vec(0, counter2) = a0;
+    a_coeff_vec(1, counter1) = a_coeff_vec(1, counter2) = a1;
+}
+
 // ------------Molar Thermodynamic Properties -------------------------
 
 doublereal RedlichKwongMFTP::enthalpy_mole() const
@@ -554,10 +603,6 @@ bool RedlichKwongMFTP::addSpecies(shared_ptr<Species> spec)
 
         a_coeff_vec.resize(2, m_kk * m_kk, 0.0);
 
-        m_pc_Species.push_back(0.0);
-        m_tc_Species.push_back(0.0);
-        m_vc_Species.push_back(0.0);
-
         m_pp.push_back(0.0);
         m_tmpV.push_back(0.0);
         m_partialMolarVolumes.push_back(0.0);
@@ -568,201 +613,101 @@ bool RedlichKwongMFTP::addSpecies(shared_ptr<Species> spec)
 
 void RedlichKwongMFTP::initThermoXML(XML_Node& phaseNode, const std::string& id)
 {
-    // Check the model parameter for the Redlich-Kwong equation of state
-    // two are allowed
-    //       RedlichKwong        mixture of species, each of which are RK fluids
-    //       RedlichKwongMFTP    mixture of species with cross term coefficients
     if (phaseNode.hasChild("thermo")) {
         XML_Node& thermoNode = phaseNode.child("thermo");
         std::string model = thermoNode["model"];
-        if (model == "RedlichKwong") {
+        if (model != "RedlichKwong" && model != "RedlichKwongMFTP") {
-            m_standardMixingRules = 1;
-        } else if (model == "RedlichKwongMFTP") {
-            m_standardMixingRules = 0;
-        } else {
             throw CanteraError("RedlichKwongMFTP::initThermoXML",
                                "Unknown thermo model : " + model);
         }
 
         // Go get all of the coefficients and factors in the
         // activityCoefficients XML block
-        XML_Node* acNodePtr = 0;
         if (thermoNode.hasChild("activityCoefficients")) {
             XML_Node& acNode = thermoNode.child("activityCoefficients");
-            acNodePtr = &acNode;
-            size_t nC = acNode.nChildren();
 
             // Loop through the children getting multiple instances of
             // parameters
-            for (size_t i = 0; i < nC; i++) {
+            for (size_t i = 0; i < acNode.nChildren(); i++) {
-                XML_Node& xmlACChild = acNodePtr->child(i);
+                XML_Node& xmlACChild = acNode.child(i);
                 if (ba::iequals(xmlACChild.name(), "purefluidparameters")) {
                     readXMLPureFluid(xmlACChild);
-                }
+                } else if (ba::iequals(xmlACChild.name(), "crossfluidparameters")) {
-            }
-            if (m_standardMixingRules == 1) {
-                applyStandardMixingRules();
-            }
-
-            // Loop through the children getting multiple instances of
-            // parameters
-            for (size_t i = 0; i < nC; i++) {
-                XML_Node& xmlACChild = acNodePtr->child(i);
-                if (ba::iequals(xmlACChild.name(), "crossfluidparameters")) {
                     readXMLCrossFluid(xmlACChild);
                 }
             }
         }
     }
 
-    for (size_t i = 0; i < m_kk; i++) {
-        double a0coeff = a_coeff_vec(0, i*m_kk + i);
-        double aTcoeff = a_coeff_vec(1, i*m_kk + i);
-        double ai = a0coeff + aTcoeff * 500.;
-        double bi = b_vec_Curr_[i];
-        calcCriticalConditions(ai, bi, a0coeff, aTcoeff, m_pc_Species[i], m_tc_Species[i], m_vc_Species[i]);
-    }
-
     MixtureFugacityTP::initThermoXML(phaseNode, id);
 }
 
 void RedlichKwongMFTP::readXMLPureFluid(XML_Node& pureFluidParam)
 {
-    vector_fp vParams;
     string xname = pureFluidParam.name();
     if (xname != "pureFluidParameters") {
         throw CanteraError("RedlichKwongMFTP::readXMLPureFluid",
                            "Incorrect name for processing this routine: " + xname);
     }
 
-    // Read the species. Find the index of the species in the current phase.
+    double a0 = 0.0;
-    // It's not an error to not find the species
+    double a1 = 0.0;
-    string iName = pureFluidParam.attrib("species");
+    double b = 0.0;
-    if (iName == "") {
+    for (size_t iChild = 0; iChild < pureFluidParam.nChildren(); iChild++) {
-        throw CanteraError("RedlichKwongMFTP::readXMLPureFluid", "no species attribute");
-    }
-    size_t iSpecies = speciesIndex(iName);
-    if (iSpecies == npos) {
-        return;
-    }
-    size_t counter = iSpecies + m_kk * iSpecies;
-    size_t nParamsExpected, nParamsFound;
-    size_t num = pureFluidParam.nChildren();
-    for (size_t iChild = 0; iChild < num; iChild++) {
         XML_Node& xmlChild = pureFluidParam.child(iChild);
         string nodeName = ba::to_lower_copy(xmlChild.name());
 
         if (nodeName == "a_coeff") {
+            vector_fp vParams;
             string iModel = ba::to_lower_copy(xmlChild.attrib("model"));
-            if (iModel == "constant") {
-                nParamsExpected = 1;
-            } else if (iModel == "linear_a") {
-                nParamsExpected = 2;
-                if (m_formTempParam == 0) {
-                    m_formTempParam = 1;
-                }
-            } else {
-                throw CanteraError("RedlichKwongMFTP::readXMLPureFluid", "unknown model");
-            }
-
             getFloatArray(xmlChild, vParams, true, "Pascal-m6/kmol2", "a_coeff");
-            nParamsFound = vParams.size();
+
-            if (nParamsFound != nParamsExpected) {
+            if (iModel == "constant" && vParams.size() == 1) {
-                throw CanteraError("RedlichKwongMFTP::readXMLPureFluid(for a_coeff" + iName + ")",
+                a0 = vParams[0];
-                                   "wrong number of params found");
+                a1 = 0;
+            } else if (iModel == "linear_a" && vParams.size() == 2) {
+                a0 = vParams[0];
+                a1 = vParams[1];
+            } else {
+                throw CanteraError("RedlichKwongMFTP::readXMLPureFluid",
+                    "unknown model or incorrect number of parameters");
             }
 
-            for (size_t i = 0; i < nParamsFound; i++) {
-                a_coeff_vec(i, counter) = vParams[i];
-            }
         } else if (nodeName == "b_coeff") {
-            getFloatArray(xmlChild, vParams, true, "m3/kmol", "b_coeff");
+            b = getFloatCurrent(xmlChild, "toSI");
-            nParamsFound = vParams.size();
-            if (nParamsFound != 1) {
-                throw CanteraError("RedlichKwongMFTP::readXMLPureFluid(for b_coeff" + iName + ")",
-                                   "wrong number of params found");
-            }
-            b_vec_Curr_[iSpecies] = vParams[0];
-        }
-    }
-}
-
-void RedlichKwongMFTP::applyStandardMixingRules()
-{
-    int nParam = 2;
-    for (size_t i = 0; i < m_kk; i++) {
-        size_t icounter = i + m_kk * i;
-        for (size_t j = 0; j < m_kk; j++) {
-            if (i != j) {
-                size_t counter = i + m_kk * j;
-                size_t jcounter = j + m_kk * j;
-                for (int n = 0; n < nParam; n++) {
-                    a_coeff_vec(n, counter) = sqrt(a_coeff_vec(n, icounter) * a_coeff_vec(n, jcounter));
-                }
-            }
         }
     }
+    setSpeciesCoeffs(pureFluidParam.attrib("species"), a0, a1, b);
 }
 
 void RedlichKwongMFTP::readXMLCrossFluid(XML_Node& CrossFluidParam)
 {
-    vector_fp vParams;
     string xname = CrossFluidParam.name();
     if (xname != "crossFluidParameters") {
         throw CanteraError("RedlichKwongMFTP::readXMLCrossFluid",
                            "Incorrect name for processing this routine: " + xname);
     }
 
-    // Read the species. Find the index of the species in the current phase.
-    // It's not an error to not find the species
     string iName = CrossFluidParam.attrib("species1");
-    if (iName == "") {
-        throw CanteraError("RedlichKwongMFTP::readXMLCrossFluid", "no species1 attribute");
-    }
-    size_t iSpecies = speciesIndex(iName);
-    if (iSpecies == npos) {
-        return;
-    }
     string jName = CrossFluidParam.attrib("species2");
-    if (iName == "") {
-        throw CanteraError("RedlichKwongMFTP::readXMLCrossFluid", "no species2 attribute");
-    }
-    size_t jSpecies = speciesIndex(jName);
-    if (jSpecies == npos) {
-        return;
-    }
 
-    size_t counter = iSpecies + m_kk * jSpecies;
-    size_t counter0 = jSpecies + m_kk * iSpecies;
-    size_t nParamsExpected, nParamsFound;
     size_t num = CrossFluidParam.nChildren();
     for (size_t iChild = 0; iChild < num; iChild++) {
         XML_Node& xmlChild = CrossFluidParam.child(iChild);
         string nodeName = ba::to_lower_copy(xmlChild.name());
 
         if (nodeName == "a_coeff") {
-            string iModel = ba::to_lower_copy(xmlChild.attrib("model"));
+            vector_fp vParams;
-            if (iModel == "constant") {
-                nParamsExpected = 1;
-            } else if (iModel == "linear_a") {
-                nParamsExpected = 2;
-                if (m_formTempParam == 0) {
-                    m_formTempParam = 1;
-                }
-            } else {
-                throw CanteraError("RedlichKwongMFTP::readXMLCrossFluid", "unknown model");
-            }
-
             getFloatArray(xmlChild, vParams, true, "Pascal-m6/kmol2", "a_coeff");
-            nParamsFound = vParams.size();
+            string iModel = ba::to_lower_copy(xmlChild.attrib("model"));
-            if (nParamsFound != nParamsExpected) {
+            if (iModel == "constant" && vParams.size() == 1) {
-                throw CanteraError("RedlichKwongMFTP::readXMLCrossFluid(for a_coeff" + iName + ")",
+                setBinaryCoeffs(iName, jName, vParams[0], 0.0);
-                                   "wrong number of params found");
+            } else if (iModel == "linear_a") {
-            }
+                setBinaryCoeffs(iName, jName, vParams[0], vParams[1]);
-
+            } else {
-            for (size_t i = 0; i < nParamsFound; i++) {
+                throw CanteraError("RedlichKwongMFTP::readXMLCrossFluid",
-                a_coeff_vec(i, counter) = vParams[i];
+                    "unknown model ({}) or wrong number of parameters ({})",
-                a_coeff_vec(i, counter0) = vParams[i];
+                    iModel, vParams.size());
             }
         }
     }

test/thermo/phaseConstructors.cpp

@@ -3,6 +3,7 @@
 #include "cantera/thermo/FixedChemPotSSTP.h"
 #include "cantera/thermo/PureFluidPhase.h"
 #include "cantera/thermo/WaterSSTP.h"
+#include "cantera/thermo/RedlichKwongMFTP.h"
 #include "cantera/thermo/NasaPoly2.h"
 #include "cantera/thermo/ShomatePoly.h"
 #include "cantera/thermo/IdealGasPhase.h"
@@ -114,15 +115,15 @@ public:
         sO2->thermo.reset(new NasaPoly2(200, 3500, 101325, o2_nasa_coeffs));
         sOH->thermo.reset(new NasaPoly2(200, 3500, 101325, oh_nasa_coeffs));
         sCO->thermo.reset(new NasaPoly2(200, 3500, 101325, o2_nasa_coeffs));
-        sCO2->thermo.reset(new NasaPoly2(200, 3500, 101325, h2o_nasa_coeffs));
+        sCO2->thermo.reset(new ShomatePoly2(200, 3500, 101325, co2_shomate_coeffs));
     }
 
-    IdealGasPhase p;
     shared_ptr<Species> sH2O, sH2, sO2, sOH, sCO, sCO2;
 };
 
 TEST_F(ConstructFromScratch, AddElements)
 {
+    IdealGasPhase p;
     p.addElement("H");
     p.addElement("O");
     ASSERT_EQ((size_t) 2, p.nElements());
@@ -132,6 +133,7 @@ TEST_F(ConstructFromScratch, AddElements)
 
 TEST_F(ConstructFromScratch, AddSpeciesDefaultBehavior)
 {
+    IdealGasPhase p;
     p.addElement("H");
     p.addElement("O");
     p.addSpecies(sH2O);
@@ -151,6 +153,7 @@ TEST_F(ConstructFromScratch, AddSpeciesDefaultBehavior)
 
 TEST_F(ConstructFromScratch, ignoreUndefinedElements)
 {
+    IdealGasPhase p;
     p.addElement("H");
     p.addElement("O");
     p.ignoreUndefinedElements();
@@ -168,6 +171,7 @@ TEST_F(ConstructFromScratch, ignoreUndefinedElements)
 
 TEST_F(ConstructFromScratch, addUndefinedElements)
 {
+    IdealGasPhase p;
     p.addElement("H");
     p.addElement("O");
     p.addUndefinedElements();
@@ -187,6 +191,29 @@ TEST_F(ConstructFromScratch, addUndefinedElements)
     ASSERT_DOUBLE_EQ(0.5, p.massFraction("CO2"));
 }
 
+TEST_F(ConstructFromScratch, RedlichKwongMFTP)
+{
+    RedlichKwongMFTP p;
+    p.addUndefinedElements();
+    p.addSpecies(sCO2);
+    p.addSpecies(sH2O);
+    p.addSpecies(sH2);
+    double fa = toSI("bar-cm6/mol2");
+    double fb = toSI("cm3/mol");
+    p.setBinaryCoeffs("H2", "H2O", 4 * fa, 40 * fa);
+    p.setSpeciesCoeffs("CO2", 7.54e7 * fa, -4.13e4 * fa, 27.80 * fb);
+    p.setBinaryCoeffs("CO2", "H2O", 7.897e7 * fa, 0.0);
+    p.setSpeciesCoeffs("H2O", 1.7458e8 * fa, -8e4 * fa, 18.18 * fb);
+    p.setSpeciesCoeffs("H2", 30e7 * fa, -330e4 * fa, 31 * fb);
+    p.initThermo();
+    p.setMoleFractionsByName("CO2:0.9998, H2O:0.0002");
+    p.setState_TP(300, 200 * OneAtm);
+    EXPECT_NEAR(p.pressure(), 200 * OneAtm, 1e-5);
+    // Arbitrary regression test values
+    EXPECT_NEAR(p.density(), 892.421, 2e-3);
+    EXPECT_NEAR(p.enthalpy_mole(), -404848642.3797, 1e-3);
+}
+
 TEST(PureFluidFromScratch, CarbonDioxide)
 {
     PureFluidPhase p;