From 9e98aef969d83f142d07fc4e594ceaa41da2ef15 Mon Sep 17 00:00:00 2001 From: Ray Speth Date: Sat, 2 Feb 2019 22:36:44 -0500 Subject: [PATCH] [Input] Create RedlichKwongMFTP objects from YAML definitions --- include/cantera/base/AnyMap.h | 1 + include/cantera/thermo/RedlichKwongMFTP.h | 1 + src/base/AnyMap.cpp | 4 + src/thermo/RedlichKwongMFTP.cpp | 147 +++++++++++++++------- src/thermo/ThermoFactory.cpp | 1 + test/data/thermo-models.yaml | 55 ++++++++ test/thermo/thermoFromYaml.cpp | 13 ++ 7 files changed, 174 insertions(+), 48 deletions(-) diff --git a/include/cantera/base/AnyMap.h b/include/cantera/base/AnyMap.h index 7cf5f84d4..8a809dfb2 100644 --- a/include/cantera/base/AnyMap.h +++ b/include/cantera/base/AnyMap.h @@ -64,6 +64,7 @@ public: template bool is() const; + bool isScalar() const; explicit AnyValue(const std::string& value); explicit AnyValue(const char* value); diff --git a/include/cantera/thermo/RedlichKwongMFTP.h b/include/cantera/thermo/RedlichKwongMFTP.h index 0fa3448e0..f6a0faac8 100644 --- a/include/cantera/thermo/RedlichKwongMFTP.h +++ b/include/cantera/thermo/RedlichKwongMFTP.h @@ -177,6 +177,7 @@ public: virtual void setParametersFromXML(const XML_Node& thermoNode); virtual void setToEquilState(const doublereal* lambda_RT); virtual void initThermoXML(XML_Node& phaseNode, const std::string& id); + virtual void initThermo(); //! Retrieve a and b coefficients by looking up tabulated critical parameters /*! diff --git a/src/base/AnyMap.cpp b/src/base/AnyMap.cpp index 859ca4953..e22bb9b4a 100644 --- a/src/base/AnyMap.cpp +++ b/src/base/AnyMap.cpp @@ -304,6 +304,10 @@ std::string AnyValue::type_str() const { return demangle(type()); } +bool AnyValue::isScalar() const { + return is() || is() || is() || is(); +} + AnyValue::AnyValue(const std::string& value) : m_value(new boost::any{value}) {} AnyValue::AnyValue(const char* value) : m_value(new boost::any{std::string(value)}) {} diff --git a/src/thermo/RedlichKwongMFTP.cpp b/src/thermo/RedlichKwongMFTP.cpp index 2ce92c0e9..2d4b426aa 100644 --- a/src/thermo/RedlichKwongMFTP.cpp +++ b/src/thermo/RedlichKwongMFTP.cpp @@ -577,74 +577,125 @@ void RedlichKwongMFTP::initThermoXML(XML_Node& phaseNode, const std::string& id) "Unknown thermo model : " + model); } + // Reset any coefficients which may have been set using values from + // 'critProperties.xml' as part of non-XML initialization, so that + // off-diagonal elements can be correctly initialized + a_coeff_vec.data().assign(a_coeff_vec.data().size(), NAN); + // Go get all of the coefficients and factors in the // activityCoefficients XML block if (thermoNode.hasChild("activityCoefficients")) { XML_Node& acNode = thermoNode.child("activityCoefficients"); - // Count the number of species with parameters provided in the - // input file: - size_t nParams = 0; - // Loop through the children and read out fluid parameters. Process // all the pureFluidParameters, first: - for (size_t i = 0; i < acNode.nChildren(); i++) { - XML_Node& xmlACChild = acNode.child(i); - if (caseInsensitiveEquals(xmlACChild.name(), "purefluidparameters")) { - readXMLPureFluid(xmlACChild); - nParams += 1; - } - } - - // If any species exist which have undefined pureFluidParameters, - // search the database in 'critProperties.xml' to find critical - // temperature and pressure to calculate a and b. - - // Loop through all species in the CTI file - size_t iSpecies = 0; - - for (size_t i = 0; i < m_kk; i++) { - string iName = speciesName(i); - - // Get the index of the species - iSpecies = speciesIndex(iName); - - // Check if a and b are already populated (only the diagonal elements of a). - size_t counter = iSpecies + m_kk * iSpecies; - - // If not, then search the database: - if (isnan(a_coeff_vec(0, counter))) { - - vector coeffArray; - - // Search the database for the species name and calculate - // coefficients a and b, from critical properties: - // coeffArray[0] = a0, coeffArray[1] = b; - coeffArray = getCoeff(iName); - - // Check if species was found in the database of critical properties, - // and assign the results - if (!isnan(coeffArray[0])) { - //Assuming no temperature dependence (i,e a1 = 0) - setSpeciesCoeffs(iName, coeffArray[0], 0.0, coeffArray[1]); - } - } - } - // Loop back through the "activityCoefficients" children and process the // crossFluidParameters in the XML tree: for (size_t i = 0; i < acNode.nChildren(); i++) { XML_Node& xmlACChild = acNode.child(i); - if (caseInsensitiveEquals(xmlACChild.name(), "crossfluidparameters")) { + if (caseInsensitiveEquals(xmlACChild.name(), "purefluidparameters")) { + readXMLPureFluid(xmlACChild); + } else if (caseInsensitiveEquals(xmlACChild.name(), "crossfluidparameters")) { readXMLCrossFluid(xmlACChild); } } } + // If any species exist which have undefined pureFluidParameters, + // search the database in 'critProperties.xml' to find critical + // temperature and pressure to calculate a and b. + + // Loop through all species in the CTI file + size_t iSpecies = 0; + + for (size_t i = 0; i < m_kk; i++) { + string iName = speciesName(i); + + // Get the index of the species + iSpecies = speciesIndex(iName); + + // Check if a and b are already populated (only the diagonal elements of a). + size_t counter = iSpecies + m_kk * iSpecies; + + // If not, then search the database: + if (isnan(a_coeff_vec(0, counter))) { + + vector coeffArray; + + // Search the database for the species name and calculate + // coefficients a and b, from critical properties: + // coeffArray[0] = a0, coeffArray[1] = b; + coeffArray = getCoeff(iName); + + // Check if species was found in the database of critical properties, + // and assign the results + if (!isnan(coeffArray[0])) { + //Assuming no temperature dependence (i,e a1 = 0) + setSpeciesCoeffs(iName, coeffArray[0], 0.0, coeffArray[1]); + } + } + } } MixtureFugacityTP::initThermoXML(phaseNode, id); } +void RedlichKwongMFTP::initThermo() +{ + for (auto& item : m_species) { + // Read a and b coefficients from species 'input' information (i.e. as + // specified in a YAML input file) + if (item.second->input.hasKey("equation-of-state")) { + auto eos = item.second->input["equation-of-state"].as(); + if (eos.getString("model", "") != "Redlich-Kwong") { + throw CanteraError("RedlichKwongMFTP::initThermo", + "Expected species equation of state to be 'Redlich-Kwong', " + "but got '{}' instead", eos.getString("model", "")); + } + double a0 = 0, a1 = 0; + if (eos["a"].isScalar()) { + a0 = eos.convert("a", "Pa*m^6/kmol^2*K^0.5"); + } else { + auto avec = eos["a"].asVector(2); + a0 = eos.units().convert(avec[0], "Pa*m^6/kmol^2*K^0.5"); + a1 = eos.units().convert(avec[1], "Pa*m^6/kmol^2/K^0.5"); + } + double b = eos.convert("b", "m^3/kmol"); + setSpeciesCoeffs(item.first, a0, a1, b); + if (eos.hasKey("binary-a")) { + AnyMap& binary_a = eos["binary-a"].as(); + const UnitSystem& units = binary_a.units(); + for (auto& item2 : binary_a) { + double a0 = 0, a1 = 0; + if (item2.second.isScalar()) { + a0 = units.convert(item2.second, "Pa*m^6/kmol^2*K^0.5"); + } else { + auto avec = item2.second.asVector(2); + a0 = units.convert(avec[0], "Pa*m^6/kmol^2*K^0.5"); + a1 = units.convert(avec[1], "Pa*m^6/kmol^2/K^0.5"); + } + setBinaryCoeffs(item.first, item2.first, a0, a1); + } + } + } else { + // Check if a and b are already populated for this species (only the + // diagonal elements of a). If not, then search 'critProperties.xml' + // to find critical temperature and pressure to calculate a and b. + size_t k = speciesIndex(item.first); + if (isnan(a_coeff_vec(0, k + m_kk * k))) { + // coeffs[0] = a0, coeffs[1] = b; + vector coeffs = getCoeff(item.first); + + // Check if species was found in the database of critical + // properties, and assign the results + if (!isnan(coeffs[0])) { + // Assuming no temperature dependence (i.e. a1 = 0) + setSpeciesCoeffs(item.first, coeffs[0], 0.0, coeffs[1]); + } + } + } + } +} + vector RedlichKwongMFTP::getCoeff(const std::string& iName) { vector_fp spCoeff{NAN, NAN}; diff --git a/src/thermo/ThermoFactory.cpp b/src/thermo/ThermoFactory.cpp index 4934143b0..fdb359b72 100644 --- a/src/thermo/ThermoFactory.cpp +++ b/src/thermo/ThermoFactory.cpp @@ -81,6 +81,7 @@ ThermoFactory::ThermoFactory() reg("Redlich-Kister", []() { return new RedlichKisterVPSSTP(); }); reg("RedlichKwong", []() { return new RedlichKwongMFTP(); }); m_synonyms["RedlichKwongMFTP"] = "RedlichKwong"; + m_synonyms["Redlich-Kwong"] = "RedlichKwong"; reg("MaskellSolidSolnPhase", []() { return new MaskellSolidSolnPhase(); }); m_synonyms["Maskell-solid-solution"] = "MaskellSolidSolnPhase"; reg("PureLiquidWater", []() { return new WaterSSTP(); }); diff --git a/test/data/thermo-models.yaml b/test/data/thermo-models.yaml index 6ec10db3d..792a8f82c 100644 --- a/test/data/thermo-models.yaml +++ b/test/data/thermo-models.yaml @@ -169,6 +169,12 @@ phases: A_Debye: variable interactions: *hmw-ineractions +- name: CO2-RK + species: [{rk-species: [CO2, H2O, H2]}] + thermo: Redlich-Kwong + state: {T: 300, P: 200 atm, mole-fractions: {CO2: 0.9998, H2O: 0.0002}} + + species: - name: NaCl(s) composition: {Na: 1, Cl: 1} @@ -428,6 +434,55 @@ gas-species: - [3.28253784, 0.00148308754, -7.57966669e-07, 2.09470555e-10, -2.16717794e-14, -1088.45772, 5.45323129] +rk-species: +- name: H2 + composition: {H: 2} + thermo: + model: NASA7 + temperature-ranges: [200.0, 1000.0, 3500.0] + data: + - [2.34433112, 0.00798052075, -1.9478151e-05, 2.01572094e-08, -7.37611761e-12, + -917.935173, 0.683010238] + - [3.3372792, -4.94024731e-05, 4.99456778e-07, -1.79566394e-10, 2.00255376e-14, + -950.158922, -3.20502331] + equation-of-state: + model: Redlich-Kwong + units: {length: cm, pressure: bar, quantity: mol} + a: [3.0e+08, -3.30e+06] + b: 31.0 +- name: H2O + composition: {H: 2, O: 1} + thermo: + model: NASA7 + temperature-ranges: [200.0, 1000.0, 3500.0] + data: + - [4.19864056, -0.0020364341, 6.52040211e-06, -5.48797062e-09, 1.77197817e-12, + -30293.7267, -0.849032208] + - [3.03399249, 0.00217691804, -1.64072518e-07, -9.7041987e-11, 1.68200992e-14, + -30004.2971, 4.9667701] + equation-of-state: + model: Redlich-Kwong + a: [1.7458e+08 bar*cm^6/mol^2*K^0.5, -8.0e+04 bar*cm^6/mol^2*K^-0.5] + b: 18.18 cm^3/mol + binary-a: + H2: [4 bar*cm^6/mol^2*K^0.5, 40 bar*cm^6/mol^2*K^-0.5] + CO2: 7.897e7 bar*cm^6/mol^2*K^0.5 +- name: CO2 + composition: {C: 1, O: 2} + thermo: + model: NASA7 + temperature-ranges: [200.0, 1000.0, 3500.0] + data: + - [2.35677352, 0.00898459677, -7.12356269e-06, 2.45919022e-09, -1.43699548e-13, + -48371.9697, 9.90105222] + - [3.85746029, 0.00441437026, -2.21481404e-06, 5.23490188e-10, -4.72084164e-14, + -48759.166, 2.27163806] + equation-of-state: + model: Redlich-Kwong + units: {length: cm, pressure: bar, quantity: mol} + a: [7.54e7, -4.13e4] + b: 27.80 + HMW-species: - name: H2O(L) diff --git a/test/thermo/thermoFromYaml.cpp b/test/thermo/thermoFromYaml.cpp index ad410d588..b5667313f 100644 --- a/test/thermo/thermoFromYaml.cpp +++ b/test/thermo/thermoFromYaml.cpp @@ -311,3 +311,16 @@ TEST(ThermoFromYaml, HMWSoln_HKFT) EXPECT_NEAR(acoeff[k], acoeffRef[k], 2e-8); } } + +TEST(ThermoFromYaml, RedlichKwong_CO2) +{ + auto thermo = newThermo("thermo-models.yaml", "CO2-RK"); + EXPECT_NEAR(thermo->density(), 892.420938853, 1e-8); + EXPECT_NEAR(thermo->enthalpy_mass(), -9199743.7500511, 1e-6); + EXPECT_NEAR(thermo->cp_mass(), 2219.899777820, 1e-8); + + thermo->setState_TPX(350, 180*OneAtm, "CO2:0.6, H2O:0.02, H2:0.38"); + EXPECT_NEAR(thermo->density(), 181.567887542, 1e-8); + EXPECT_NEAR(thermo->enthalpy_mass(), -8872890.9496462, 1e-6); + EXPECT_NEAR(thermo->cp_mass(), 3358.439021094, 1e-8); +}