mirror of
https://github.com/Cantera/cantera.git
synced 2025-02-25 18:55:29 -06:00
updates to the radiation classes and starting to add documentation
This commit is contained in:
Christopher Neal
parent
3b8c0ce594
commit
cb1aec990c
@ -15,51 +15,303 @@
|
||||
namespace Cantera
|
||||
{
|
||||
|
||||
/**
|
||||
* Computes the radiative heat loss vector over points jmin to jmax and stores
|
||||
* the data in the qdotRadiation variable.
|
||||
/** 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 `fit-type` of `polynomial` is uses the model described below.
|
||||
* 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.
|
||||
*
|
||||
* The simple radiation model used was established by Liu and Rogg
|
||||
* @cite liu1991. This model considers the radiation of CO2 and H2O.
|
||||
* 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.
|
||||
*
|
||||
* 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 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].
|
||||
* 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 Radiation1D {
|
||||
class TabularPlanckMean : public RadiationPropertyCalculator {
|
||||
public:
|
||||
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);
|
||||
/**
|
||||
* 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);
|
||||
|
||||
// Parse radiation data from YAML input
|
||||
/** 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 radiative heat loss
|
||||
void computeRadiation(double* x, size_t jmin, size_t jmax, vector<double>& qdotRadiation);
|
||||
/** 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);
|
||||
|
||||
//! Set the emissivities for the boundary values
|
||||
/*!
|
||||
* Reads the emissivities for the left and right boundary values in the
|
||||
* radiative term and writes them into the variables, which are used for the
|
||||
* calculation.
|
||||
/** 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);
|
||||
|
||||
@ -74,25 +326,21 @@ public:
|
||||
}
|
||||
|
||||
private:
|
||||
ThermoPhase* m_thermo;
|
||||
double m_press;
|
||||
size_t m_points;
|
||||
ThermoPhase* m_thermo; //!< Pointer to the ThermoPhase object
|
||||
double m_press; //!< Pressure in Pa
|
||||
size_t m_points; //!< Number of grid points
|
||||
|
||||
map<string, size_t> m_absorptionSpecies; //!< Absorbing species
|
||||
AnyMap m_PMAC; //!< Absorption coefficient data for each species
|
||||
//! Property calculator for absorption coefficients
|
||||
std::unique_ptr<RadiationPropertyCalculator> m_props;
|
||||
|
||||
//! Emissivity of the surface to the left of the domain. Used for calculating
|
||||
//! 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;
|
||||
|
||||
//! Emissivity of the surface to the right of the domain. Used for calculating
|
||||
//! radiative heat loss.
|
||||
double m_epsilon_right = 0.0;
|
||||
|
||||
// Helper functions
|
||||
double calculatePolynomial(const vector<double>& coefficients, double temperature);
|
||||
double interpolateTable(const vector<double>& temperatures, const vector<double>& data, double temperature);
|
||||
|
||||
//! Lambda function to get temperature at a given point
|
||||
std::function<double(const double*, size_t)> m_T;
|
||||
|
||||
@ -100,6 +348,70 @@ private:
|
||||
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
|
||||
|
||||
@ -12,38 +12,25 @@
|
||||
namespace Cantera
|
||||
{
|
||||
|
||||
|
||||
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)
|
||||
: m_thermo(thermo), m_press(pressure), m_points(points),
|
||||
m_T(temperatureFunction), m_X(moleFractionFunction)
|
||||
TabularPlanckMean::TabularPlanckMean(ThermoPhase* thermo)
|
||||
: m_thermo(thermo)
|
||||
{
|
||||
parseRadiationData();
|
||||
}
|
||||
|
||||
void Radiation1D::setBoundaryEmissivities(double e_left, double e_right)
|
||||
void TabularPlanckMean::parseRadiationData()
|
||||
{
|
||||
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::parseRadiationData() {
|
||||
AnyMap radiationPropertiesDB;
|
||||
// Search 'crit-properties.yaml' to find Tc and Pc. Load data if needed.
|
||||
if (radiationPropertiesDB.empty()) {
|
||||
|
||||
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.hasKey("PMAC")) {
|
||||
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
|
||||
@ -92,50 +79,29 @@ void Radiation1D::parseRadiationData() {
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// 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;
|
||||
}
|
||||
// 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;
|
||||
}
|
||||
// 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;
|
||||
}
|
||||
}
|
||||
|
||||
void Radiation1D::computeRadiation(double* x, size_t jmin, size_t jmax, std::vector<double>& qdotRadiation) {
|
||||
const double k_P_ref = 1.0 * OneAtm;
|
||||
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++) {
|
||||
double k_P = 0;
|
||||
for (const auto& [sp_name, sp_idx] : m_absorptionSpecies) {
|
||||
const auto& fit_type = m_PMAC[sp_name]["fit-type"].asString();
|
||||
if (fit_type == "table") {
|
||||
k_P += m_press * m_X(x, sp_idx, j) *
|
||||
interpolateTable(m_PMAC[sp_name]["temperatures"].asVector<double>(), m_PMAC[sp_name]["coefficients"].asVector<double>(), m_T(x, j));
|
||||
} else if (fit_type == "polynomial") {
|
||||
k_P += m_press * m_X(x, sp_idx, j) *
|
||||
calculatePolynomial(m_PMAC[sp_name]["coefficients"].asVector<double>(), m_T(x, j));
|
||||
}
|
||||
}
|
||||
qdotRadiation[j] = 2 * k_P * (2 * StefanBoltz * std::pow(m_T(x, j), 4) - boundary_Rad_left - boundary_Rad_right);
|
||||
}
|
||||
}
|
||||
|
||||
double Radiation1D::calculatePolynomial(const std::vector<double>& coefficients, double temperature) {
|
||||
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));
|
||||
@ -143,18 +109,133 @@ double Radiation1D::calculatePolynomial(const std::vector<double>& coefficients,
|
||||
return result / (1.0 * OneAtm);
|
||||
}
|
||||
|
||||
double Radiation1D::interpolateTable(const std::vector<double>& temperatures, const std::vector<double>& data, double temperature) {
|
||||
size_t index = 0;
|
||||
for (size_t i = 1; i < temperatures.size(); ++i) {
|
||||
if (temperature < temperatures[i]) {
|
||||
index = i - 1;
|
||||
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;
|
||||
}
|
||||
}
|
||||
double t1 = temperatures[index], t2 = temperatures[index + 1];
|
||||
double d1 = data[index], d2 = data[index + 1];
|
||||
return d1 + (d2 - d1) * (temperature - t1) / (t2 - t1);
|
||||
|
||||
// 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
|
||||
Loading…
Reference in New Issue
Block a user