GreenPolyProducts/cantera
5
0
Fork 0
mirror of https://github.com/Cantera/cantera.git synced 2025-02-25 18:55:29 -06:00
Code Issues Packages Projects Releases Wiki Activity

Merge cb1aec990c into 0720efb02d

Browse Source
This commit is contained in:
Chris Neal 2025-02-13 20:38:13 -06:00 committed by GitHub
commit e64595924f
No known key found for this signature in database
GPG Key ID: B5690EEEBB952194
6 changed files with 1027 additions and 76 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

28
include/cantera/oneD/Flow1D.h
View File

@ -7,11 +7,13 @@
#define CT_FLOW1D_H
#include "Domain1D.h"
#include "Radiation1D.h"
#include "cantera/base/Array.h"
#include "cantera/base/Solution.h"
#include "cantera/thermo/ThermoPhase.h"
#include "cantera/kinetics/Kinetics.h"
namespace Cantera
{
@ -259,12 +261,12 @@ public:
//! Return emissivity at left boundary
double leftEmissivity() const {
return m_epsilon_left;
return m_radiation->leftEmissivity();
}
//! Return emissivity at right boundary
double rightEmissivity() const {
return m_epsilon_right;
return m_radiation->rightEmissivity();
}
//! Specify that the the temperature should be held fixed at point `j`.
@ -461,22 +463,8 @@ protected:
//! to be updated are defined.
virtual void updateProperties(size_t jg, double* x, size_t jmin, size_t jmax);
/**
* Computes the radiative heat loss vector over points jmin to jmax and stores
* the data in the qdotRadiation variable.
*
* The simple radiation model used was established by Liu and Rogg
* @cite liu1991. This model considers the radiation of CO2 and H2O.
*
* This model uses the optically thin limit and the gray-gas approximation to
* simply calculate a volume specified heat flux out of the Planck absorption
* coefficients, the boundary emissivities and the temperature. Polynomial lines
* calculate the species Planck coefficients for H2O and CO2. The data for the
* lines are taken from the RADCAL program @cite RADCAL.
* The coefficients for the polynomials are taken from
* [TNF Workshop](https://tnfworkshop.org/radiation/) material.
*/
void computeRadiation(double* x, size_t jmin, size_t jmax);
//! Compute the radiative heat loss at each grid point
void computeRadiation(double*, size_t, size_t);
//! @}
@ -902,6 +890,9 @@ protected:
//! radiative heat loss.
double m_epsilon_right = 0.0;
//! Radiation object used for calculating radiative heat loss
std::unique_ptr<Radiation1D> m_radiation;
//! Indices within the ThermoPhase of the radiating species. First index is
//! for CO2, second is for H2O.
vector<size_t> m_kRadiating;
@ -1006,5 +997,4 @@ private:
};
}
#endif

417
include/cantera/oneD/Radiation1D.h Normal file
View File

@ -0,0 +1,417 @@
//! @file Radiation1D.h
// 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.
#ifndef RADIATION1D_H
#define RADIATION1D_H
#include "Domain1D.h"
#include "cantera/base/Array.h"
#include "cantera/thermo/ThermoPhase.h"
#include <functional>
namespace Cantera
{
/** Stores the temperature, pressure, an optional soot fraction (fvSoot),
* and a map of species mole fractions. Allows the property calculator
* classes to retrieve the local state information they need .
*
* The temperature is given in Kelvin [K], the pressure in Pascals [Pa],
* and fvSoot is a dimensionless volume fraction for soot. The map 'X'
* holds species names as keys and their mole fractions (unitless) as values.
*/
struct RadComposition {
double T = 0.0; //! Temperature (K)
double P = 0.0; //! Pressure (Pa)
double fvSoot = 0.0; //! Soot volume fraction
Composition X; //! Map of name->mole fraction
};
/** Base class for radiation property calculators.
*
* Responsible for computing the spectral absorption coefficients (kabs) and
* weighting factors (awts) for a given thermodynamic state. Different models
* e.g. polynomial fits, tabular data, or external libraries such as RadLib
* are implemented by deriving from this class.
*
* The data produced by getBandProperties() are used by a RadiationSolver to compute
* the net radiative heat loss at each point in a 1D domain.
*/
class RadiationPropertyCalculator {
public:
virtual ~RadiationPropertyCalculator() = default;
/** Calculate absorption coefficients and weighting factors for each spectral band.
*
* The size of `kabs` and `awts` determine how many "gray gases" or bands are
* used. For a simple Planck mean approach, there may be only one band. For
* multi-band or weighted-sum-of-grey-gases (WSGG) models, there could be multiple.
*
* @param kabs A vector to be filled with absorption coefficients (k_i).
* @param awts A vector to be filled with weighting factors (a_i).
* @param comp A RadComposition containing T, P, composition, etc.
*/
virtual void getBandProperties(std::vector<double>& kabs,
std::vector<double>& awts,
const RadComposition& comp) = 0;
};
/* Commented out for now until RadLib dependency is resolved
class RadLibPlanckMean : public RadiationPropertyCalculator {
public:
RadLibPlanckMean() {
m_rad = new rad_planck_mean();
}
~RadLibPlanckMean() {
delete m_rad;
}
void getBandProperties(std::vector<double>& kabs,
std::vector<double>& awts,
double T, double P, const RadComposition& comp) override
{
m_rad->get_k_a(kabs, awts, T, P, comp.fvSoot, comp.xH2O, comp.xCO2, comp.xCO, comp.xCH4);
}
private:
rad* m_rad; // pointer to rad_planck_mean
};
*/
/* Reads species-specific Planck-mean absorption coefficient data from a YAML file
* (radiation-properties.yaml) if available, falling back to polynomial approximations
* for CO2 and H2O otherwise. A `fit-type` of "table" or "polynomial" can be specified
* in the YAML data.
*
* The table-based data uses interpolation for a discrete set of temperatures,
* whereas the polynomial data uses functional fits.
*
* The `fit-type` of `polynomial` is uses the model described below:
*
* Polynomial lines calculate the species Planck coefficients for H2O and CO2. The
* data for the lines are taken from the RADCAL program @cite RADCAL.
* The coefficients for the polynomials are taken from
* [TNF Workshop](https://tnfworkshop.org/radiation/) material.
*
* The `fit-type` of `table` is uses the model described below.
*
* Spectra for molecules are downloaded with HAPI library from // https://hitran.org/hapi/
* [R.V. Kochanov, I.E. Gordon, L.S. Rothman, P. Wcislo, C. Hill, J.S. Wilzewski,
* HITRAN Application Programming Interface (HAPI): A comprehensive approach
* to working with spectroscopic data, J. Quant. Spectrosc. Radiat. Transfer 177,
* 15-30 (2016), https://doi.org/10.1016/j.jqsrt.2016.03.005].
*
* Planck mean optical path lengths are what are read in from a YAML input file.
*/
class TabularPlanckMean : public RadiationPropertyCalculator {
public:
/**
* The constructor will attempt to parse radiation data from
* "radiation-properties.yaml". If that file doesn't exist, a warning is
* issued and polynomial defaults for CO2 and H2O are used.
*
* @param thermo Pointer to a ThermoPhase object which provides species names
* and other properties. This is needed to match species found in the
* YAML database to the actual species in the simulation.
*/
TabularPlanckMean(ThermoPhase* thermo);
/** Calculate absorption coefficients and weighting factors for each band.
* This method sums absorption contributions from all absorbing species for
* which the table or polynomial data is defined. The final result is stored
* as a single-band coefficient (kabs.size()==1), with awts.size()==1=1.0,
* representing a gray approximation.
*
* @param kabs A vector to store absorption coefficients (k_i).
* @param awts A vector to store weighting factors (a_i).
* @param comp The RadComposition struct with T, P, and species mole fractions.
*/
void getBandProperties(std::vector<double>& kabs, std::vector<double>& awts,
const RadComposition& comp) override;
private:
/** Parse optional YAML data from "radiation-properties.yaml".
* If the file is not found, a warning is issued and default polynomial data
* for H2O and CO2 is used. If it is found, then species listed in the file
* are read into 'm_PMAC' along with their "fit-type" and associated
* coefficients. This might be polynomial or tabulated data.
*
* The method also ensures that H2O and CO2 have some default data even if
* the file does not provide them.
*/
void parseRadiationData();
/** Compute polynomial-based absorption coefficient.
*
* This evaluates a polynomial that has a form given as:
*
* kabs = c0 + c1(1000/T) + c2(1000/T)^2 + c3(1000/T)^3 + c4(1000/T)^4 + c5(1000/T)^5
*
* The value computed is the Plank mean absorption coefficient. This is just one way
* to represent the variation of the absorption coefficient with temperature, and
* it used often for species such as H2O and CO2. The units of the coefficients are
* (m-1 atm-1) and the temperature is in Kelvin( see https://tnfworkshop.org/radiation/).
* This function converts the units to m-1 Pa-1.
*
* @param coefficients The polynomial coefficients for the species.
* @param temperature The local temperature in K.
* @return The Planck-mean absorption coefficient for that species
* at the given temperature in units of m-1 Pa-1.
*/
double calculatePolynomial(const vector<double>& coefficients, double temperature);
/** Compute an absorption coefficient using a log-linear interpolation.
*
* Use log-linear interpolation to compute an absorption coefficient from a
* table of Planck mean optical-path-length values ('data') at discrete
* 'temperatures'. Units of the optical path length are meters and the units of
* temperature are Kelvin.
*
* The tables hold values of the optical path length (OPL) for a gas at
* different temperatures. The absorption coefficient is the inverse of the
* OPL i.e. kabs = 1.0 / OPL
*
* The method uses the following algorithm:
* alpha = 1.0 / data[i]
* ln(alpha) is interpolated linearly vs. T in the bracket
* [temperatures[i-1], temperatures[i]] using the following formula:
*
* ln(alpha) = ln(1/v1) + ( ln(1/v2) - ln(1/v1) ) * (T - t1)/(t2 - t1)
*
* If T is below the lowest or above the highest table entry, the boundary value
* without interpolation is used.
*
* @param temperatures Sorted vector of tabulated temperatures
* @param data Corresponding tabulated OPL data
* (optical path lengths)
* @param temperature Query temperature
* @returns The absorption coefficient at 'temperature'
*
* @param data A vector of corresponding absorption data.
* @return The interpolated absorption coefficient at the given temperature
* in units of m-1 Pa-1.
*/
double interpolateTable(const vector<double>& temperatures,
const vector<double>& data, double temperature);
private:
ThermoPhase* m_thermo; //!< Pointer to the ThermoPhase object
map<string, size_t> m_absorptionSpecies; //!< Absorbing species mapping names to indices
AnyMap m_PMAC; //!< Absorption coefficient data for each species
};
/** Base class for radiation solvers.
*
* Computes the net radiative heat loss (or gain) from
* absorption coefficients, weighting factors, boundary emissivities,
* and local temperature. Different solver implementations should derive from
* this class.
*/
class RadiationSolver {
public:
virtual ~RadiationSolver() = default;
// compute the radiative heat loss given boundary conditions, geometry, etc.
virtual double computeHeatLoss(const std::vector<double>& kabs,
const std::vector<double>& awts,
double T,
double boundaryRadLeft,
double boundaryRadRight) = 0;
};
/*
* The simple radiation model used was established by Liu and Rogg
* @cite liu1991. This model considers the radiation of CO2 and H2O.
*
* This model uses the *optically thin limit* and the *gray-gas approximation*,
* calculating a volumetric heat loss (qdot) using Planck-mean absorption
* coefficients, boundary emissivities, and local temperature.
*
* Typically, qdot for each spectral band i is given by:
*
* \f[
* \dot{q}_i = 2 k_i \left( 2 \sigma T^4 - E_{\text{left}} - E_{\text{right}} \right)
* \f]
*
* Summed over all bands i, weighted by a_i. The 2 factor comes from the assumption
* that radiation can escape in both directions in a 1D domain, ignoring scattering.
*/
class OpticallyThinSolver : public RadiationSolver {
public:
//! Calculate optically thin radiative heat loss for each band, summing the
//! contributions.
double computeHeatLoss(const std::vector<double>& kabs,
const std::vector<double>& awts,
double T,
double boundaryRadLeft,
double boundaryRadRight) override
{
// Sum over each band.
double sum = 0.0;
double sigma = 5.67e-8; // Stefan-Boltzmann constant
for (size_t i=0; i<kabs.size(); ++i) {
sum += awts[i] * 2.0 * kabs[i] *
(2.0 * sigma * std::pow(T,4) - boundaryRadLeft - boundaryRadRight);
}
return sum;
}
};
/**
* The Radiation1D class ties together a RadiationPropertyCalculator and a
* RadiationSolver to compute volumetric radiative heat losses at each grid
* point (j). It fetches the local temperature and species mole fractions
* from the solution vector (using the provided lambda functions) and
* optionally applies boundary emissivities to represent radiative flux
* at the domain edges.
*
*/
class Radiation1D {
public:
/**
* Constructor for a Radiation1D object.
*
* @param thermo A pointer to the ThermoPhase object for species info
* @param pressure The operating pressure [Pa] (assumed uniform)
* @param points Number of grid points in the domain
* @param temperatureFunction A lambda to retrieve T(x, j) from the solution
* @param moleFractionFunction A lambda to retrieve X_k(x, j) from the solution
* @param props A unique pointer to a RadiationPropertyCalculator
* @param solver A unique pointer to a RadiationSolver
*/
Radiation1D(ThermoPhase* thermo, double pressure, size_t points,
std::function<double(const double*, size_t)> temperatureFunction,
std::function<double(const double*, size_t, size_t)> moleFractionFunction,
std::unique_ptr<RadiationPropertyCalculator> props,
std::unique_ptr<RadiationSolver> solver);
/** Compute radiative heat loss from jmin to jmax and fill the qdotRadiation array.
*
* This method extracts T and mole fractions from the solution at each grid point.
* It then uses the RadiationPropertyCalculator to get the band properties, and then
* uses the RadiationSolver to get the heat loss and stores it in qdotRadiation.
*
* @param x Pointer to the solution vector for the 1D domain.
* @param jmin The first grid index to compute.
* @param jmax One past the last grid index to compute.
* @param qdotRadiation A vector of size at least jmax, which will be filled
* with the volumetric radiative heat loss at each grid point.
*/
void computeRadiation(double* x, size_t jmin, size_t jmax,
vector<double>& qdotRadiation);
/**
* Sets the emissivities for the left and right boundary values in the
* radiative term.
*/
void setBoundaryEmissivities(double e_left, double e_right);
//! Return emissivity at left boundary
double leftEmissivity() const {
return m_epsilon_left;
}
//! Return emissivity at right boundary
double rightEmissivity() const {
return m_epsilon_right;
}
private:
ThermoPhase* m_thermo; //!< Pointer to the ThermoPhase object
double m_press; //!< Pressure in Pa
size_t m_points; //!< Number of grid points
//! Property calculator for absorption coefficients
std::unique_ptr<RadiationPropertyCalculator> m_props;
//! Solver for radiative heat loss
std::unique_ptr<RadiationSolver> m_solver;
//! Emissivity of the surface to the left and right of the domain. Used for calculating
//! radiative heat loss.
double m_epsilon_left = 0.0;
double m_epsilon_right = 0.0;
//! Lambda function to get temperature at a given point
std::function<double(const double*, size_t)> m_T;
//! Lambda function to get mole fraction at a given point
std::function<double(const double*, size_t, size_t)> m_X;
};
/**
* Create a Radiation1D instance based on the selected property model and solver.
*
* @param propertyModel String specifying which property calculator to use
* @param solverModel String specifying which solver to use
* @param thermo Pointer to ThermoPhase
* @param pressure Pressure in Pa
* @param points Number of grid points
* @param Tfunc Lambda for temperature
* @param Xfunc Lambda for mole fractions
* @param e_left Emissivity at left boundary
* @param e_right Emissivity at right boundary
*
* @return Radiation1D object
*/
inline std::unique_ptr<Radiation1D> createRadiation1D(
const std::string& propertyModel,
const std::string& solverModel,
ThermoPhase* thermo,
double pressure,
size_t points,
std::function<double(const double*, size_t)> Tfunc,
std::function<double(const double*, size_t, size_t)> Xfunc,
double e_left,
double e_right
)
{
// Create the RadiationPropertyCalculator
std::unique_ptr<RadiationPropertyCalculator> props;
if (propertyModel == "TabularPlanckMean") {
props = std::make_unique<TabularPlanckMean>(thermo);
}
else if (propertyModel == "RadLibPlanckMean") {
//props = std::make_unique<RadLibPlanckMean>();
}
else {
throw CanteraError("createRadiation1D",
"Unknown property model: " + propertyModel);
}
// Create the RadiationSolver
std::unique_ptr<RadiationSolver> solver;
if (solverModel == "OpticallyThin") {
solver = std::make_unique<OpticallyThinSolver>();
}
else {
throw CanteraError("createRadiation1D",
"Unknown solver model: " + solverModel);
}
// Build the Radiation1D object
auto rad = std::make_unique<Radiation1D>(
thermo, pressure, points,
std::move(Tfunc), std::move(Xfunc),
std::move(props), std::move(solver)
);
rad->setBoundaryEmissivities(e_left, e_right);
return rad;
};
} // namespace Cantera
#endif // RADIATION1D_H

86
src/oneD/Flow1D.cpp
View File

@ -10,6 +10,8 @@
#include "cantera/transport/TransportFactory.h"
#include "cantera/numerics/funcs.h"
#include "cantera/base/global.h"
#include "cantera/thermo/Species.h"
using namespace std;
@ -81,10 +83,31 @@ Flow1D::Flow1D(ThermoPhase* ph, size_t nsp, size_t points) :
}
setupGrid(m_points, gr.data());
// Find indices for radiating species
m_kRadiating.resize(2, npos);
m_kRadiating[0] = m_thermo->speciesIndex("CO2");
m_kRadiating[1] = m_thermo->speciesIndex("H2O");
// Initialize the radiation object (hardcoded for now)
std::string propertyModel = "TabularPlanckMean";
std::string solverModel = "OpticallyThin";
// Define lambdas for T(x, j) and X(x, k, j)
auto Tfunc = [this](const double* x, size_t j) {
return this->T(x, j);
};
auto Xfunc = [this](const double* x, size_t k, size_t j) {
return this->X(x, k, j);
};
double emissivityLeft = 0.0;
double emissivityRight = 0.0;
m_radiation = createRadiation1D(
propertyModel,
solverModel,
m_thermo,
m_thermo->pressure(),
m_points,
Tfunc,
Xfunc,
emissivityLeft,
emissivityRight
);
}
Flow1D::Flow1D(shared_ptr<ThermoPhase> th, size_t nsp, size_t points)
@ -466,49 +489,7 @@ void Flow1D::updateDiffFluxes(const double* x, size_t j0, size_t j1)
void Flow1D::computeRadiation(double* x, size_t jmin, size_t jmax)
{
// Variable definitions for the Planck absorption coefficient and the
// radiation calculation:
double k_P_ref = 1.0*OneAtm;
// Polynomial coefficients:
const double c_H2O[6] = {-0.23093, -1.12390, 9.41530, -2.99880,
0.51382, -1.86840e-5};
const double c_CO2[6] = {18.741, -121.310, 273.500, -194.050,
56.310, -5.8169};
// Calculation of the two boundary values
double boundary_Rad_left = m_epsilon_left * StefanBoltz * pow(T(x, 0), 4);
double boundary_Rad_right = m_epsilon_right * StefanBoltz * pow(T(x, m_points - 1), 4);
for (size_t j = jmin; j < jmax; j++) {
// calculation of the mean Planck absorption coefficient
double k_P = 0;
// Absorption coefficient for H2O
if (m_kRadiating[1] != npos) {
double k_P_H2O = 0;
for (size_t n = 0; n <= 5; n++) {
k_P_H2O += c_H2O[n] * pow(1000 / T(x, j), (double) n);
}
k_P_H2O /= k_P_ref;
k_P += m_press * X(x, m_kRadiating[1], j) * k_P_H2O;
}
// Absorption coefficient for CO2
if (m_kRadiating[0] != npos) {
double k_P_CO2 = 0;
for (size_t n = 0; n <= 5; n++) {
k_P_CO2 += c_CO2[n] * pow(1000 / T(x, j), (double) n);
}
k_P_CO2 /= k_P_ref;
k_P += m_press * X(x, m_kRadiating[0], j) * k_P_CO2;
}
// Calculation of the radiative heat loss term
double radiative_heat_loss = 2 * k_P *(2 * StefanBoltz * pow(T(x, j), 4)
- boundary_Rad_left - boundary_Rad_right);
// set the radiative heat loss vector
m_qdotRadiation[j] = radiative_heat_loss;
}
m_radiation->computeRadiation(x, jmin, jmax, m_qdotRadiation);
}
void Flow1D::evalContinuity(double* x, double* rsd, int* diag,
@ -1100,16 +1081,7 @@ bool Flow1D::doElectricField(size_t j) const
void Flow1D::setBoundaryEmissivities(double e_left, double e_right)
{
if (e_left < 0 || e_left > 1) {
throw CanteraError("Flow1D::setBoundaryEmissivities",
"The left boundary emissivity must be between 0.0 and 1.0!");
} else if (e_right < 0 || e_right > 1) {
throw CanteraError("Flow1D::setBoundaryEmissivities",
"The right boundary emissivity must be between 0.0 and 1.0!");
} else {
m_epsilon_left = e_left;
m_epsilon_right = e_right;
}
m_radiation->setBoundaryEmissivities(e_left, e_right);
}
void Flow1D::fixTemperature(size_t j)

241
src/oneD/Radiation1D.cpp Normal file
View File

@ -0,0 +1,241 @@
//! @file Radiation1D.cpp
// 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/oneD/Radiation1D.h"
#include "cantera/thermo/Species.h"
#include "cantera/base/global.h"
namespace Cantera
{
TabularPlanckMean::TabularPlanckMean(ThermoPhase* thermo)
: m_thermo(thermo)
{
parseRadiationData();
}
void TabularPlanckMean::parseRadiationData()
{
AnyMap radiationPropertiesDB;
try {
radiationPropertiesDB = AnyMap::fromYamlFile("radiation-properties.yaml");
} catch (CanteraError& err) {
warn_user("TabularPlanckMean::parseRadiationData",
"Failed to load 'radiation-properties.yaml':\n{}"
"\nFalling back to default polynomial data for CO2, H2O.", err.what());
}
if(!radiationPropertiesDB.empty() && radiationPropertiesDB.hasKey("PMAC")) {
auto& data = radiationPropertiesDB["PMAC"].as<AnyMap>();
// Needs to loop over only the species that are in the input yaml data
for (const auto& name : m_thermo->speciesNames()) {
if (data.hasKey("radiation")) {
std::cout << "Radiation data found for species " << name << std::endl;
m_absorptionSpecies.insert({name, m_thermo->speciesIndex(name)});
if (data["radiation"].hasKey("fit-type")) {
m_PMAC[name]["fit-type"] = data["radiation"]["fit-type"].asString();
} else {
throw InputFileError("Flow1D::Flow1D", data,
"No 'fit-type' entry found for species '{}'", name);
}
// This is the direct tabulation of the optical path length
if (data["radiation"]["fit-type"] == "table") {
if (data["radiation"].hasKey("temperatures")) {
std::cout << "Storing temperatures for species " << name << std::endl;
// Each species may have a specific set of temperatures that are used
m_PMAC[name]["temperatures"] = data["radiation"]["temperatures"].asVector<double>();
} else {
throw InputFileError("Flow1D::Flow1D", data,
"No 'temperatures' entry found for species '{}'", name);
}
if (data["radiation"].hasKey("data")) {
std::cout << "Storing data for species " << name << std::endl;
// This data is the Plank mean absorption coefficient
m_PMAC[name]["coefficients"] = data["radiation"]["data"].asVector<double>();
} else {
throw InputFileError("Flow1D::Flow1D", data,
"No 'data' entry found for species '{}'", name);
}
} else if (data["radiation"]["fit-type"] == "polynomial") {
std::cout << "Polynomial fit found for species " << name << std::endl;
if (data["radiation"].hasKey("data")) {
std::cout << "Storing data for species " << name << std::endl;
m_PMAC[name]["coefficients"] = data["radiation"]["data"].asVector<double>();
} else {
throw InputFileError("Flow1D::Flow1D", data,
"No 'data' entry found for species '{}'", name);
}
} else {
throw InputFileError("Flow1D::Flow1D", data,
"Invalid 'fit-type' entry found for species '{}'", name);
}
}
}
}
// Polynomial coefficients for CO2 and H2O (backwards compatibility)
// Check if "CO2" is already in the map, if not, add the polynomial fit data
if (!m_PMAC.hasKey("CO2")) {
const std::vector<double> c_CO2 = {18.741, -121.310, 273.500, -194.050, 56.310,
-5.8169};
m_PMAC["CO2"]["fit-type"] = "polynomial";
m_PMAC["CO2"]["coefficients"] = c_CO2;
}
// Check if "H2O" is already in the map, if not, add the polynomial fit data
if (!m_PMAC.hasKey("H2O")) {
const std::vector<double> c_H2O = {-0.23093, -1.12390, 9.41530, -2.99880,
0.51382, -1.86840e-5};
m_PMAC["H2O"]["fit-type"] = "polynomial";
m_PMAC["H2O"]["coefficients"] = c_H2O;
}
}
double TabularPlanckMean::calculatePolynomial(const std::vector<double>& coefficients,
double temperature)
{
double result = 0.0;
for (size_t n = 0; n < coefficients.size(); ++n) {
result += coefficients[n] * std::pow(1000 / temperature, static_cast<double>(n));
}
return result / (1.0 * OneAtm);
}
double TabularPlanckMean::interpolateTable(const std::vector<double>& temperatures,
const std::vector<double>& data,
double temperature)
{
// Handle edge cases first
if (temperature <= temperatures.front()) {
// alpha = 1.0 / data[0]
return 1.0 / data.front();
} else if (temperature >= temperatures.back()) {
// alpha = 1.0 / data[last]
return 1.0 / data.back();
}
// Find the interval [t1, t2] where t1 <= T < t2
// so that temperatures[i-1] <= T < temperatures[i]
size_t idx = 1;
for (; idx < temperatures.size(); ++idx) {
if (temperature < temperatures[idx]) {
break;
}
}
// Perform linear interpolation
double t1 = temperatures[idx - 1];
double t2 = temperatures[idx];
double v1 = data[idx - 1];
double v2 = data[idx];
// ln(alpha) = ln(1/v1) + ( ln(1/v2) - ln(1/v1) ) * (T - t1)/(t2 - t1)
double frac = (temperature - t1) / (t2 - t1);
double lnAlpha = log(1.0 / v1) + (log(1.0 / v2) - log(1.0 / v1)) * frac;
return exp(lnAlpha) / (1.0 * OneAtm);
}
void TabularPlanckMean::getBandProperties(std::vector<double>& kabs,
std::vector<double>& awts,
const RadComposition& comp)
{
double k_P = 0;
// Loop over absorbing species
for (const auto& [sp_name, sp_idx] : m_absorptionSpecies) {
const auto& fit_type = m_PMAC[sp_name]["fit-type"].asString();
// Get the species mole fraction from the Composition
// If the species doesn't exist in comp.X, error out
double x_sp = 0;
if (comp.X.find(sp_name) == comp.X.end()) {
throw CanteraError("TabularPlanckMean::getBandProperties",
"Species '{}' not found in composition data", sp_name);
} else {
x_sp = comp.X.at(sp_name);
}
if (fit_type == "table") {
double kVal = interpolateTable(
m_PMAC[sp_name]["temperatures"].asVector<double>(),
m_PMAC[sp_name]["coefficients"].asVector<double>(),
comp.T);
k_P += comp.P * x_sp * kVal;
} else if (fit_type == "polynomial") {
double kVal = calculatePolynomial(
m_PMAC[sp_name]["coefficients"].asVector<double>(),
comp.T);
k_P += comp.P * x_sp * kVal;
}
}
// Store the single-band result
kabs.resize(1);
awts.resize(1);
kabs[0] = k_P;
awts[0] = 1.0; // single “band” weighting
}
Radiation1D::Radiation1D(ThermoPhase* thermo, double pressure, size_t points,
std::function<double(const double*, size_t)> temperatureFunction,
std::function<double(const double*, size_t, size_t)> moleFractionFunction,
std::unique_ptr<RadiationPropertyCalculator> props,
std::unique_ptr<RadiationSolver> solver)
: m_thermo(thermo), m_press(pressure), m_points(points),
m_T(temperatureFunction), m_X(moleFractionFunction),
m_props(std::move(props)), m_solver(std::move(solver))
{
}
void Radiation1D::setBoundaryEmissivities(double e_left, double e_right)
{
if (e_left < 0 || e_left > 1) {
throw CanteraError("Radiation1D::setBoundaryEmissivities",
"The left boundary emissivity must be between 0.0 and 1.0!");
} else if (e_right < 0 || e_right > 1) {
throw CanteraError("Radiation1D::setBoundaryEmissivities",
"The right boundary emissivity must be between 0.0 and 1.0!");
} else {
m_epsilon_left = e_left;
m_epsilon_right = e_right;
}
}
void Radiation1D::computeRadiation(double* x, size_t jmin, size_t jmax,
std::vector<double>& qdotRadiation) {
const double StefanBoltz = 5.67e-8;
double boundary_Rad_left = m_epsilon_left * StefanBoltz * std::pow(m_T(x, 0), 4);
double boundary_Rad_right = m_epsilon_right * StefanBoltz * std::pow(m_T(x, m_points - 1), 4);
for (size_t j = jmin; j < jmax; j++) {
RadComposition comp;
comp.T = m_T(x, j);
comp.P = m_press;
comp.X = m_thermo->getMoleFractionsByName();
// Get the band absorption coefficients and weighting factors
std::vector<double> kabs, awts;
m_props->getBandProperties(kabs, awts, comp);
// Solve for radiative heat loss
qdotRadiation[j] = m_solver->computeHeatLoss(kabs, awts, comp.T,
boundary_Rad_left,
boundary_Rad_right);
}
}
} // namespace Cantera

276
test/data/air_radiation.yaml Normal file
View File

@ -0,0 +1,276 @@
description: |-
Ideal gas properties of air. Includes several reactions among
the included species.
generator: ck2yaml
input-files: [air.inp, gri30_tran.dat]
cantera-version: 2.5.0
date: Wed, 11 Dec 2019 16:59:03 -0500
units: {length: cm, time: s, quantity: mol, activation-energy: cal/mol}
phases:
- name: air
thermo: ideal-gas
elements: [O, N, Ar]
species: [O, O2, N, NO, NO2, N2O, N2, AR, H2O CO2]
kinetics: gas
transport: mixture-averaged
state: {T: 300.0, P: 1 atm, X: {O2: 0.21, N2: 0.78, AR: 0.01}}
radiation:
- radiation-parameters.yaml
species:
- name: O
composition: {O: 1}
thermo:
model: NASA7
temperature-ranges: [200.0, 1000.0, 3500.0]
data:
- [3.1682671, -3.27931884e-03, 6.64306396e-06, -6.12806624e-09, 2.11265971e-12,
2.91222592e+04, 2.05193346]
- [2.56942078, -8.59741137e-05, 4.19484589e-08, -1.00177799e-11, 1.22833691e-15,
2.92175791e+04, 4.78433864]
note: L1/90
transport:
model: gas
geometry: atom
well-depth: 80.0
diameter: 2.75
- name: O2
composition: {O: 2}
thermo:
model: NASA7
temperature-ranges: [200.0, 1000.0, 3500.0]
data:
- [3.78245636, -2.99673416e-03, 9.84730201e-06, -9.68129509e-09, 3.24372837e-12,
-1063.94356, 3.65767573]
- [3.28253784, 1.48308754e-03, -7.57966669e-07, 2.09470555e-10, -2.16717794e-14,
-1088.45772, 5.45323129]
note: TPIS89
transport:
model: gas
geometry: linear
well-depth: 107.4
diameter: 3.458
polarizability: 1.6
rotational-relaxation: 3.8
- name: N
composition: {N: 1}
thermo:
model: NASA7
temperature-ranges: [200.0, 1000.0, 6000.0]
data:
- [2.5, 0.0, 0.0, 0.0, 0.0, 5.6104637e+04, 4.1939087]
- [2.4159429, 1.7489065e-04, -1.1902369e-07, 3.0226245e-11, -2.0360982e-15,
5.6133773e+04, 4.6496096]
note: L6/88
transport:
model: gas
geometry: atom
well-depth: 71.4
diameter: 3.298
note: '*'
- name: NO
composition: {N: 1, O: 1}
thermo:
model: NASA7
temperature-ranges: [200.0, 1000.0, 6000.0]
data:
- [4.2184763, -4.638976e-03, 1.1041022e-05, -9.3361354e-09, 2.803577e-12,
9844.623, 2.2808464]
- [3.2606056, 1.1911043e-03, -4.2917048e-07, 6.9457669e-11, -4.0336099e-15,
9920.9746, 6.3693027]
note: RUS78
transport:
model: gas
geometry: linear
well-depth: 97.53
diameter: 3.621
polarizability: 1.76
rotational-relaxation: 4.0
radiation:
model: PMAC
fit-type: table
temperatures: [200.0, 300.0, 400.0, 500.0, 600.0,
700.0, 800.0, 900.0, 1000.0, 1100.0, 1200.0, 1300.0, 1400.0, 1500.0, 1600.0,
1700.0, 1800.0, 1900.0, 2000.0, 2100.0, 2200.0, 2300.0, 2400.0, 2500.0, 2600.0,
2700.0, 2800.0, 2900.0, 3000.0, 3100.0, 3200.0, 3300.0, 3400.0, 3500.0]
data: [12.49908285125374, 1.3459019533042522,
0.6011825839900713, 0.4748113234371179, 0.47742862256124935, 0.5370442220062525,
0.6376193053372756, 0.7764870585324466, 0.9555723091275765, 1.1788559185177765,
1.4514629725987682, 1.7793270399233352, 2.1689573586466553, 2.627448906362707,
3.1623084968864976, 3.781539187660205, 4.49349909489309, 5.307064481132861,
6.231329265582485, 7.275885774773834, 8.450611521315984, 9.765632286526694,
11.231430006248385, 12.858500531681175, 14.657751390671171, 16.639987005055776,
18.816075851493107, 21.19670280326081, 23.793015918675906, 26.61491214156301,
29.673133893532462, 32.97754031606446, 36.53746836667566, 40.3622966223194]
- name: NO2
composition: {N: 1, O: 2}
thermo:
model: NASA7
temperature-ranges: [200.0, 1000.0, 6000.0]
data:
- [3.9440312, -1.585429e-03, 1.6657812e-05, -2.0475426e-08, 7.8350564e-12,
2896.6179, 6.3119917]
- [4.8847542, 2.1723956e-03, -8.2806906e-07, 1.574751e-10, -1.0510895e-14,
2316.4983, -0.11741695]
note: L7/88
transport:
model: gas
geometry: nonlinear
well-depth: 200.0
diameter: 3.5
rotational-relaxation: 1.0
note: '*'
radiation:
model: PMAC
fit-type: table
temperatures: [200.0, 300.0, 400.0, 500.0, 600.0,
700.0, 800.0, 900.0, 1000.0, 1100.0, 1200.0, 1300.0, 1400.0, 1500.0, 1600.0,
1700.0, 1800.0, 1900.0, 2000.0, 2100.0, 2200.0, 2300.0, 2400.0, 2500.0, 2600.0,
2700.0, 2800.0, 2900.0, 3000.0, 3100.0, 3200.0, 3300.0, 3400.0, 3500.0]
data: [0.18967185355199417, 0.04428846239231481,
0.028416106886409966, 0.028275414773864766, 0.03422377575827348, 0.04560041986270607,
0.06381585795362235, 0.09150569222640514, 0.13253556724406976, 0.1922144503571911,
0.27758651879037544, 0.39778260471688276, 0.564395391544332, 0.791908677828003,
1.0981494235689724, 1.5047895936423843, 2.037869062204293, 2.728368349256181,
3.6128356924408402, 4.7340233488511325, 6.141594179781417, 7.892863536806141,
10.053585199881873, 12.698770218856987, 15.91363338637731, 19.79442325583844,
24.449499989156376, 30.000309241474096, 36.58245963378199, 44.34696368197642,
53.461298650892346, 64.1107214629735, 76.49956856745997, 90.85268289351858]
- name: N2O
composition: {N: 2, O: 1}
thermo:
model: NASA7
temperature-ranges: [200.0, 1000.0, 6000.0]
data:
- [2.2571502, 0.011304728, -1.3671319e-05, 9.6819806e-09, -2.9307182e-12,
8741.7744, 10.757992]
- [4.8230729, 2.6270251e-03, -9.5850874e-07, 1.6000712e-10, -9.7752303e-15,
8073.4048, -2.2017207]
note: L7/88
transport:
model: gas
geometry: linear
well-depth: 232.4
diameter: 3.828
rotational-relaxation: 1.0
note: '*'
radiation:
model: PMAC
fit-type: table
temperatures: [200.0, 300.0, 400.0, 500.0, 600.0,
700.0, 800.0, 900.0, 1000.0, 1100.0, 1200.0, 1300.0, 1400.0, 1500.0, 1600.0,
1700.0, 1800.0, 1900.0, 2000.0, 2100.0, 2200.0, 2300.0, 2400.0, 2500.0, 2600.0,
2700.0, 2800.0, 2900.0, 3000.0, 3100.0, 3200.0, 3300.0, 3400.0, 3500.0]
data: [0.05604186150530222, 0.04991781435183458,
0.04001827500393808, 0.036150506950346904, 0.03710719640348422, 0.04180788677648103,
0.05018622510260605, 0.06285657185631298, 0.08094410959148458, 0.10604926499576005,
0.1402734389982863, 0.18627028306484095, 0.24730302890664432, 0.32730217346266893,
0.430908689506439, 0.5635079200362725, 0.7312449464709093, 0.9410221482388516,
1.2004777818389483, 1.5179552157303737, 1.9024578243125905, 2.363587182232489,
2.911482794626652, 3.5567676931747707, 4.310436165387461, 5.18386757038981,
6.188658945204697, 7.336641004525838, 8.639751162834212, 10.110065918784443,
11.759665166149045, 13.600597157732668, 15.644893691074058, 17.9045229031322]
- name: N2
composition: {N: 2}
thermo:
model: NASA7
temperature-ranges: [300.0, 1000.0, 5000.0]
data:
- [3.298677, 1.4082404e-03, -3.963222e-06, 5.641515e-09, -2.444854e-12,
-1020.8999, 3.950372]
- [2.92664, 1.4879768e-03, -5.68476e-07, 1.0097038e-10, -6.753351e-15,
-922.7977, 5.980528]
note: '121286'
transport:
model: gas
geometry: linear
well-depth: 97.53
diameter: 3.621
polarizability: 1.76
rotational-relaxation: 4.0
- name: AR
composition: {Ar: 1}
thermo:
model: NASA7
temperature-ranges: [300.0, 1000.0, 5000.0]
data:
- [2.5, 0.0, 0.0, 0.0, 0.0, -745.375, 4.366]
- [2.5, 0.0, 0.0, 0.0, 0.0, -745.375, 4.366]
note: '120186'
transport:
model: gas
geometry: atom
well-depth: 136.5
diameter: 3.33
- name: H2O
composition: {H: 2, O: 1}
thermo:
model: NASA7
temperature-ranges: [200.0, 1000.0, 3500.0]
data:
- [4.19864056, -2.0364341e-03, 6.52040211e-06, -5.48797062e-09, 1.77197817e-12,
-3.02937267e+04, -0.849032208]
- [3.03399249, 2.17691804e-03, -1.64072518e-07, -9.7041987e-11, 1.68200992e-14,
-3.00042971e+04, 4.9667701]
note: L8/89
transport:
model: gas
geometry: nonlinear
well-depth: 572.4
diameter: 2.605
dipole: 1.844
rotational-relaxation: 4.0
radiation:
model: PMAC
fit-type: polynomial
data: [-0.23093, -1.12390, 9.41530, -2.99880, 0.51382, -1.86840e-5]
- name: CO2
composition: {C: 1, O: 2}
thermo:
model: NASA7
temperature-ranges: [200.0, 1000.0, 3500.0]
data:
- [2.35677352, 8.98459677e-03, -7.12356269e-06, 2.45919022e-09, -1.43699548e-13,
-4.83719697e+04, 9.90105222]
- [3.85746029, 4.41437026e-03, -2.21481404e-06, 5.23490188e-10, -4.72084164e-14,
-4.8759166e+04, 2.27163806]
note: L7/88
transport:
model: gas
geometry: linear
well-depth: 244.0
diameter: 3.763
polarizability: 2.65
rotational-relaxation: 2.1
radiation:
model: PMAC
fit-type: polynomial
data: [18.741, -121.310, 273.500, -194.050, 56.310, -5.8169]
reactions:
- equation: 2 O + M <=> O2 + M # Reaction 1
type: three-body
rate-constant: {A: 1.2e+17, b: -1.0, Ea: 0.0}
efficiencies: {AR: 0.83}
- equation: N + NO <=> N2 + O # Reaction 2
rate-constant: {A: 2.7e+13, b: 0.0, Ea: 355.0}
- equation: N + O2 <=> NO + O # Reaction 3
rate-constant: {A: 9.0e+09, b: 1.0, Ea: 6500.0}
- equation: N2O + O <=> N2 + O2 # Reaction 4
rate-constant: {A: 1.4e+12, b: 0.0, Ea: 1.081e+04}
- equation: N2O + O <=> 2 NO # Reaction 5
rate-constant: {A: 2.9e+13, b: 0.0, Ea: 2.315e+04}
- equation: N2O (+M) <=> N2 + O (+M) # Reaction 6
type: falloff
low-P-rate-constant: {A: 6.37e+14, b: 0.0, Ea: 5.664e+04}
high-P-rate-constant: {A: 7.91e+10, b: 0.0, Ea: 5.602e+04}
efficiencies: {AR: 0.625}
- equation: NO + O + M <=> NO2 + M # Reaction 7
type: three-body
rate-constant: {A: 1.06e+20, b: -1.41, Ea: 0.0}
efficiencies: {AR: 0.7}
- equation: NO2 + O <=> NO + O2 # Reaction 8
rate-constant: {A: 3.9e+12, b: 0.0, Ea: -240.0}

55
test/data/radiation-parameters.yaml Normal file
View File

@ -0,0 +1,55 @@
# radiation-parameters.yaml
PMAC:
- NO:
fit-type: table
temperatures: [200.0, 300.0, 400.0, 500.0, 600.0,
700.0, 800.0, 900.0, 1000.0, 1100.0, 1200.0, 1300.0, 1400.0, 1500.0, 1600.0,
1700.0, 1800.0, 1900.0, 2000.0, 2100.0, 2200.0, 2300.0, 2400.0, 2500.0, 2600.0,
2700.0, 2800.0, 2900.0, 3000.0, 3100.0, 3200.0, 3300.0, 3400.0, 3500.0]
data: [12.49908285125374, 1.3459019533042522,
0.6011825839900713, 0.4748113234371179, 0.47742862256124935, 0.5370442220062525,
0.6376193053372756, 0.7764870585324466, 0.9555723091275765, 1.1788559185177765,
1.4514629725987682, 1.7793270399233352, 2.1689573586466553, 2.627448906362707,
3.1623084968864976, 3.781539187660205, 4.49349909489309, 5.307064481132861,
6.231329265582485, 7.275885774773834, 8.450611521315984, 9.765632286526694,
11.231430006248385, 12.858500531681175, 14.657751390671171, 16.639987005055776,
18.816075851493107, 21.19670280326081, 23.793015918675906, 26.61491214156301,
29.673133893532462, 32.97754031606446, 36.53746836667566, 40.3622966223194]
- NO2:
fit-type: table
temperatures: [200.0, 300.0, 400.0, 500.0, 600.0,
700.0, 800.0, 900.0, 1000.0, 1100.0, 1200.0, 1300.0, 1400.0, 1500.0, 1600.0,
1700.0, 1800.0, 1900.0, 2000.0, 2100.0, 2200.0, 2300.0, 2400.0, 2500.0, 2600.0,
2700.0, 2800.0, 2900.0, 3000.0, 3100.0, 3200.0, 3300.0, 3400.0, 3500.0]
data: [0.18967185355199417, 0.04428846239231481,
0.028416106886409966, 0.028275414773864766, 0.03422377575827348, 0.04560041986270607,
0.06381585795362235, 0.09150569222640514, 0.13253556724406976, 0.1922144503571911,
0.27758651879037544, 0.39778260471688276, 0.564395391544332, 0.791908677828003,
1.0981494235689724, 1.5047895936423843, 2.037869062204293, 2.728368349256181,
3.6128356924408402, 4.7340233488511325, 6.141594179781417, 7.892863536806141,
10.053585199881873, 12.698770218856987, 15.91363338637731, 19.79442325583844,
24.449499989156376, 30.000309241474096, 36.58245963378199, 44.34696368197642,
53.461298650892346, 64.1107214629735, 76.49956856745997, 90.85268289351858]
- N2O:
fit-type: table
temperatures: [200.0, 300.0, 400.0, 500.0, 600.0,
700.0, 800.0, 900.0, 1000.0, 1100.0, 1200.0, 1300.0, 1400.0, 1500.0, 1600.0,
1700.0, 1800.0, 1900.0, 2000.0, 2100.0, 2200.0, 2300.0, 2400.0, 2500.0, 2600.0,
2700.0, 2800.0, 2900.0, 3000.0, 3100.0, 3200.0, 3300.0, 3400.0, 3500.0]
data: [0.05604186150530222, 0.04991781435183458,
0.04001827500393808, 0.036150506950346904, 0.03710719640348422, 0.04180788677648103,
0.05018622510260605, 0.06285657185631298, 0.08094410959148458, 0.10604926499576005,
0.1402734389982863, 0.18627028306484095, 0.24730302890664432, 0.32730217346266893,
0.430908689506439, 0.5635079200362725, 0.7312449464709093, 0.9410221482388516,
1.2004777818389483, 1.5179552157303737, 1.9024578243125905, 2.363587182232489,
2.911482794626652, 3.5567676931747707, 4.310436165387461, 5.18386757038981,
6.188658945204697, 7.336641004525838, 8.639751162834212, 10.110065918784443,
11.759665166149045, 13.600597157732668, 15.644893691074058, 17.9045229031322]
- H2O:
fit-type: polynomial
data: [-0.23093, -1.12390, 9.41530, -2.99880, 0.51382, -1.86840e-5]
- CO2:
fit-type: polynomial
data: [18.741, -121.310, 273.500, -194.050, 56.310, -5.8169]