471 lines
14 KiB
C
471 lines
14 KiB
C
/*---------------------------------------------------------------------------*\
|
|
========= |
|
|
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
|
|
\\ / O peration |
|
|
\\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation
|
|
\\/ M anipulation |
|
|
-------------------------------------------------------------------------------
|
|
License
|
|
This file is part of OpenFOAM.
|
|
|
|
OpenFOAM is free software: you can redistribute it and/or modify it
|
|
under the terms of the GNU General Public License as published by
|
|
the Free Software Foundation, either version 3 of the License, or
|
|
(at your option) any later version.
|
|
|
|
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
|
|
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
|
|
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
|
|
for more details.
|
|
|
|
You should have received a copy of the GNU General Public License
|
|
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
|
|
|
|
Application
|
|
chemFoam
|
|
|
|
Description
|
|
Solver for chemistry problems, designed for use on single cell cases to
|
|
provide comparison against other chemistry solvers, that uses a single cell
|
|
mesh, and fields created from the initial conditions.
|
|
|
|
\*---------------------------------------------------------------------------*/
|
|
|
|
#include <cantera/transport.h>
|
|
#include <cantera/IdealGasMix.h>
|
|
#include <cantera/thermo/speciesThermoTypes.h>
|
|
|
|
#include <fstream>
|
|
|
|
#include <cmath>
|
|
|
|
#include "cxxopts/cxxopts.hpp"
|
|
|
|
#include "ofFormats.h"
|
|
|
|
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
|
|
|
|
std::string stringNASACoefs (doublereal *coef)
|
|
{
|
|
std::string arr;
|
|
|
|
for (int j = 0; j < 7; j++)
|
|
{
|
|
arr += fmt::format(" {:15.10E} ", coef[j]);
|
|
}
|
|
|
|
return arr;
|
|
}
|
|
|
|
void calculateSutherland (Cantera::IdealGasMix *gas_, doublereal As[], doublereal C[])
|
|
{
|
|
std::shared_ptr<Cantera::Transport> tr_ ( Cantera::newTransportMgr("Mix", gas_) );
|
|
|
|
doublereal mu0[gas_->nSpecies()];
|
|
doublereal muRatio[gas_->nSpecies()];
|
|
|
|
const doublereal T0 = gas_->minTemp();
|
|
gas_->setState_TP(T0, Cantera::OneAtm);
|
|
tr_->getSpeciesViscosities(mu0);
|
|
|
|
size_t nInterval = 100;
|
|
doublereal T[nInterval];
|
|
doublereal TRatioPower[nInterval];
|
|
doublereal dTemp = (gas_->maxTemp() - gas_->minTemp()) / doublereal(nInterval);
|
|
|
|
for (size_t i = 0; i < nInterval; i++)
|
|
{
|
|
T[i] = gas_->minTemp() + (i+1) * dTemp;
|
|
}
|
|
|
|
for (size_t i = 0; i < nInterval; i++)
|
|
{
|
|
TRatioPower[i] = pow(T[i]/T0, 3./2.);
|
|
}
|
|
|
|
for (size_t j = 0; j < gas_->nSpecies(); j++)
|
|
{
|
|
C[j] = 0.0;
|
|
}
|
|
|
|
for (size_t i = 0; i < nInterval; i++)
|
|
{
|
|
doublereal Ti = T[i];
|
|
|
|
gas_->setTemperature(Ti);
|
|
|
|
tr_->getSpeciesViscosities(muRatio);
|
|
for (size_t j = 0; j < gas_->nSpecies(); j++)
|
|
{
|
|
muRatio[j] /= mu0[j];
|
|
}
|
|
|
|
for (size_t j = 0; j < gas_->nSpecies(); j++)
|
|
{
|
|
C[j] += (muRatio[j] * Ti - TRatioPower[i] * T0) / (TRatioPower[i] - muRatio[j]);
|
|
}
|
|
}
|
|
|
|
for (size_t j = 0; j < gas_->nSpecies(); j++)
|
|
{
|
|
C[j] /= doublereal(nInterval);
|
|
}
|
|
|
|
|
|
for (size_t j = 0; j < gas_->nSpecies(); j++)
|
|
{
|
|
As[j] = mu0[j] * (T0 + C[j]) / pow(T0, 3./2.);
|
|
}
|
|
|
|
return ;
|
|
}
|
|
|
|
|
|
std::string ofReactionString (const std::shared_ptr<Cantera::Reaction> r)
|
|
{
|
|
|
|
std::ostringstream reaction;
|
|
|
|
|
|
for (auto iter = r->reactants.begin(); iter != r->reactants.end(); ++iter)
|
|
{
|
|
if (iter != r->reactants.begin())
|
|
{
|
|
reaction << " + ";
|
|
}
|
|
if (iter->second != 1.0)
|
|
{
|
|
reaction << iter->second;
|
|
}
|
|
reaction << iter->first;
|
|
if (iter->second != r->orders[iter->first] && r->orders[iter->first] != 0.0 )
|
|
{
|
|
reaction << "^" << r->orders[iter->first];
|
|
}
|
|
}
|
|
|
|
reaction << " = ";
|
|
|
|
for (Cantera::Composition::iterator iter = r->products.begin(); iter != r->products.end(); ++iter)
|
|
{
|
|
if (iter != r->products.begin())
|
|
{
|
|
reaction << " + ";
|
|
}
|
|
if (iter->second != 1.0)
|
|
{
|
|
reaction << iter->second;
|
|
}
|
|
reaction << iter->first;
|
|
if (iter->second != r->orders[iter->first] && r->orders[iter->first] != 0.0)
|
|
{
|
|
reaction << "^" << r->orders[iter->first];
|
|
}
|
|
}
|
|
|
|
return reaction.str();
|
|
}
|
|
|
|
|
|
bool is_file_exist(const std::string &fileName)
|
|
{
|
|
std::ifstream infile(fileName);
|
|
return infile.good();
|
|
}
|
|
|
|
|
|
int main(int argc, char *argv[])
|
|
{
|
|
|
|
cxxopts::Options options("ct2foam", "Utility Converts Cantera formatted data in OpenFOAM format");
|
|
|
|
std::string scti;
|
|
std::string sgas("");
|
|
std::string sthermo;
|
|
std::string srxn;
|
|
|
|
std::ofstream fthermo;
|
|
std::ofstream frxn;
|
|
|
|
options.add_options()
|
|
("name", "Name of ideal gas mix in CTI file", cxxopts::value(sgas))
|
|
("w,overwrite", "Overwrite exsisting files")
|
|
("h,help", "print help message")
|
|
;
|
|
|
|
options.add_options("internal")
|
|
("cti", "Input CTI File name", cxxopts::value(scti))
|
|
("reaction", "Output OpenFOAM Reaction File name", cxxopts::value(srxn)->default_value("reactions"))
|
|
("thermo", "Output OpenFOAM Thermo File name", cxxopts::value(sthermo)->default_value("thermos"))
|
|
("rest", "Input CTI File name", cxxopts::value<std::vector<std::string>>())
|
|
;
|
|
|
|
options.parse_positional({"cti", "reaction", "thermo", "rest"});
|
|
|
|
options.positional_help("CTI [REACTION [THERMO]]");
|
|
|
|
cxxopts::ParseResult result = options.parse(argc, argv);
|
|
|
|
if (result.count("help"))
|
|
{
|
|
std::cout << options.help();
|
|
|
|
return 0;
|
|
}
|
|
|
|
if (!result["cti"].count())
|
|
{
|
|
std::cerr << "CTI file name must be provided" << std::endl << std::endl;
|
|
|
|
std::cout << options.help();
|
|
|
|
return -1;
|
|
}
|
|
|
|
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
|
|
|
|
Cantera::IdealGasMix gas_(scti, sgas);
|
|
|
|
bool overwrite = result["overwrite"].as<bool>();
|
|
|
|
if (is_file_exist(sthermo) && (!overwrite))
|
|
{
|
|
throw std::system_error(EEXIST, std::system_category(), sthermo);
|
|
}
|
|
else
|
|
{
|
|
fthermo.open(sthermo);
|
|
}
|
|
|
|
if (is_file_exist(srxn) && (!overwrite))
|
|
{
|
|
throw std::system_error(EEXIST, std::system_category(), srxn);
|
|
}
|
|
else
|
|
{
|
|
frxn.open(srxn);
|
|
}
|
|
|
|
|
|
// Thermo Part =============================================================
|
|
|
|
// Calculate Sutherland Coefficients
|
|
doublereal As[gas_.nSpecies()] ;
|
|
doublereal C[gas_.nSpecies()] ;
|
|
|
|
calculateSutherland (&gas_, As, C);
|
|
|
|
|
|
// these constants define the location of coefficient "a6" in the
|
|
// cofficient array c. The c array contains Tmid in the first
|
|
// location, followed by the 7 low-temperature coefficients, then
|
|
// the seven high-temperature ones.
|
|
for (size_t n = 0; n < gas_.nSpecies(); n++)
|
|
{
|
|
int type;
|
|
doublereal c[15];
|
|
doublereal minTemp, maxTemp, refPressure;
|
|
|
|
// get the NASA coefficients in array c
|
|
switch(gas_.speciesThermo().reportType(n))
|
|
{
|
|
case NASA2:
|
|
gas_.speciesThermo().reportParams(n, type, c, minTemp, maxTemp, refPressure);
|
|
break;
|
|
default:
|
|
throw std::domain_error(gas_.speciesName(n) + " has a non-NASA2-type thermodynamics interpolator");
|
|
}
|
|
|
|
fthermo<<(
|
|
fmt::format(
|
|
thermoFormat,
|
|
gas_.speciesName(n),
|
|
gas_.molecularWeight(n),
|
|
gas_.charge(n),
|
|
minTemp,
|
|
maxTemp,
|
|
c[0],
|
|
stringNASACoefs(&c[1]),
|
|
stringNASACoefs(&c[8]),
|
|
As[n],
|
|
C[n]
|
|
)
|
|
);
|
|
}
|
|
|
|
fthermo.close();
|
|
|
|
// End of Thermo Part ======================================================
|
|
|
|
|
|
// Reaction Part
|
|
|
|
// Write species name list
|
|
frxn<<("species\n");
|
|
frxn<<fmt::format("{}\n", gas_.nSpecies());
|
|
frxn<<("(\n");
|
|
for (size_t n = 0; n < gas_.nSpecies(); n++)
|
|
{
|
|
frxn<<fmt::format("{}\n", gas_.speciesName(n));
|
|
}
|
|
frxn<<(")\n");
|
|
frxn<<(";\n\n");
|
|
|
|
// Write reaction list
|
|
frxn<<( "reactions\n{\n"); // begin "reactions" dictionary
|
|
|
|
for (size_t k = 0; k < gas_.nReactions(); k++)
|
|
{
|
|
std::shared_ptr<Cantera::Reaction> r(gas_.reaction(k));
|
|
|
|
std::string irn( r->reversible ? "reversible" : "irreversible" );
|
|
|
|
std::string rate("Arrhenius");
|
|
|
|
std::string rxn("Reaction");
|
|
|
|
std::string ffn;
|
|
|
|
double c[5] = {0};
|
|
|
|
// reaction name
|
|
frxn<<( fmt::format( " un-named-reaction-{}\n", k)); // begin a reaction dictionary
|
|
frxn<<( " {\n" );
|
|
|
|
switch (gas_.reactionType(k))
|
|
{
|
|
case Cantera::ELEMENTARY_RXN:
|
|
|
|
frxn<<(
|
|
fmt::format(
|
|
rxnTypeFormat + rxnEqnFormat + arrheniusFormat,
|
|
irn + rate + rxn,
|
|
ofReactionString(r),
|
|
std::dynamic_pointer_cast<Cantera::ElementaryReaction>(r)->rate.preExponentialFactor(),
|
|
std::dynamic_pointer_cast<Cantera::ElementaryReaction>(r)->rate.temperatureExponent(),
|
|
std::dynamic_pointer_cast<Cantera::ElementaryReaction>(r)->rate.activationEnergy_R()
|
|
)
|
|
);
|
|
break;
|
|
|
|
case Cantera::THREE_BODY_RXN:
|
|
|
|
rate = "thirdBody" + rate;
|
|
|
|
frxn<<(
|
|
fmt::format(
|
|
rxnTypeFormat + rxnEqnFormat + arrheniusFormat,
|
|
irn + rate + rxn,
|
|
ofReactionString(r),
|
|
std::dynamic_pointer_cast<Cantera::ElementaryReaction>(r)->rate.preExponentialFactor(),
|
|
std::dynamic_pointer_cast<Cantera::ElementaryReaction>(r)->rate.temperatureExponent(),
|
|
std::dynamic_pointer_cast<Cantera::ElementaryReaction>(r)->rate.activationEnergy_R()
|
|
)
|
|
);
|
|
|
|
frxn<<(" coeffs \n");
|
|
|
|
frxn<<(fmt::format("{}\n(\n", gas_.nSpecies()));
|
|
|
|
for (size_t l = 0; l < gas_.nSpecies(); l++)
|
|
{
|
|
frxn<<(
|
|
fmt::format(
|
|
"({} {})\n",
|
|
gas_.speciesName(l),
|
|
std::dynamic_pointer_cast<Cantera::ThreeBodyReaction>(r)->third_body.efficiency(gas_.speciesName(l))
|
|
)
|
|
);
|
|
}
|
|
|
|
frxn<<(")\n;\n");
|
|
|
|
break;
|
|
|
|
case Cantera::FALLOFF_RXN:
|
|
case Cantera::CHEMACT_RXN:
|
|
|
|
rxn = (r->reaction_type == Cantera::FALLOFF_RXN ? "FallOff" : "ChemicallyActivated") + rxn;
|
|
|
|
switch (std::dynamic_pointer_cast<Cantera::FalloffReaction>(r)->falloff->getType())
|
|
{
|
|
case Cantera::SIMPLE_FALLOFF:
|
|
rxn = "Lindemann" + rxn;
|
|
ffn = lindemannFormat;
|
|
break;
|
|
|
|
case Cantera::TROE_FALLOFF:
|
|
rxn = "Troe" + rxn;
|
|
std::dynamic_pointer_cast<Cantera::FalloffReaction>(r)->falloff->getParameters(c);
|
|
ffn = fmt::format(troeFormat, c[0], c[1], c[2], c[3]);
|
|
break;
|
|
|
|
case Cantera::SRI_FALLOFF:
|
|
rxn = "SRI" + rxn;
|
|
std::dynamic_pointer_cast<Cantera::FalloffReaction>(r)->falloff->getParameters(c);
|
|
ffn = fmt::format(sriFormat, c[0], c[1], c[2], c[3], c[4]);
|
|
break;
|
|
}
|
|
|
|
frxn<<(
|
|
fmt::format(
|
|
rxnTypeFormat + rxnEqnFormat + falloffFormat,
|
|
irn + rate + rxn,
|
|
ofReactionString(r),
|
|
std::dynamic_pointer_cast<Cantera::FalloffReaction>(r)->low_rate.preExponentialFactor(),
|
|
std::dynamic_pointer_cast<Cantera::FalloffReaction>(r)->low_rate.temperatureExponent(),
|
|
std::dynamic_pointer_cast<Cantera::FalloffReaction>(r)->low_rate.activationEnergy_R(),
|
|
std::dynamic_pointer_cast<Cantera::FalloffReaction>(r)->high_rate.preExponentialFactor(),
|
|
std::dynamic_pointer_cast<Cantera::FalloffReaction>(r)->high_rate.temperatureExponent(),
|
|
std::dynamic_pointer_cast<Cantera::FalloffReaction>(r)->high_rate.activationEnergy_R()
|
|
)
|
|
);
|
|
|
|
frxn<<( ffn );
|
|
|
|
frxn<<( std::string() +
|
|
" thirdBodyEfficiencies\n" +
|
|
" {\n" +
|
|
" coeffs \n"
|
|
);
|
|
|
|
frxn<<(fmt::format("{}\n(\n", gas_.nSpecies()));
|
|
|
|
for (size_t l = 0; l < gas_.nSpecies(); l++)
|
|
{
|
|
frxn<<(
|
|
fmt::format(
|
|
"({} {})\n",
|
|
gas_.speciesName(l),
|
|
std::dynamic_pointer_cast<Cantera::FalloffReaction>(r)->third_body.efficiency(gas_.speciesName(l))
|
|
)
|
|
);
|
|
}
|
|
|
|
frxn<<(")\n;\n");
|
|
frxn<<(" }\n");
|
|
|
|
break;
|
|
|
|
case Cantera::PLOG_RXN:
|
|
break;
|
|
|
|
case Cantera::CHEBYSHEV_RXN:
|
|
break;
|
|
|
|
default:
|
|
break;
|
|
}
|
|
|
|
frxn<<( " }\n" ); // end of a reaction dictionary
|
|
}
|
|
|
|
frxn<<( "}\n"); // end of "reactions" dictionary
|
|
|
|
frxn.close();
|
|
|
|
return 0;
|
|
}
|
|
|
|
|
|
// ************************************************************************* //
|