diff --git a/include/cantera/oneD/Radiation1D.h b/include/cantera/oneD/Radiation1D.h index 6b7ebd368..c73273911 100644 --- a/include/cantera/oneD/Radiation1D.h +++ b/include/cantera/oneD/Radiation1D.h @@ -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& kabs, + std::vector& 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& kabs, + std::vector& 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 temperatureFunction, - std::function 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& kabs, std::vector& 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& 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& 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& temperatures, + const vector& data, double temperature); + +private: + ThermoPhase* m_thermo; //!< Pointer to the ThermoPhase object + + map 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& kabs, + const std::vector& 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& kabs, + const std::vector& 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 temperatureFunction, + std::function moleFractionFunction, + std::unique_ptr props, + std::unique_ptr 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& 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 m_absorptionSpecies; //!< Absorbing species - AnyMap m_PMAC; //!< Absorption coefficient data for each species + //! Property calculator for absorption coefficients + std::unique_ptr m_props; - //! Emissivity of the surface to the left of the domain. Used for calculating + //! Solver for radiative heat loss + std::unique_ptr 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& coefficients, double temperature); - double interpolateTable(const vector& temperatures, const vector& data, double temperature); - //! Lambda function to get temperature at a given point std::function m_T; @@ -100,6 +348,70 @@ private: std::function 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 createRadiation1D( + const std::string& propertyModel, + const std::string& solverModel, + ThermoPhase* thermo, + double pressure, + size_t points, + std::function Tfunc, + std::function Xfunc, + double e_left, + double e_right +) +{ + // Create the RadiationPropertyCalculator + std::unique_ptr props; + if (propertyModel == "TabularPlanckMean") { + props = std::make_unique(thermo); + } + else if (propertyModel == "RadLibPlanckMean") { + //props = std::make_unique(); + } + else { + throw CanteraError("createRadiation1D", + "Unknown property model: " + propertyModel); + } + + // Create the RadiationSolver + std::unique_ptr solver; + if (solverModel == "OpticallyThin") { + solver = std::make_unique(); + } + else { + throw CanteraError("createRadiation1D", + "Unknown solver model: " + solverModel); + } + + // Build the Radiation1D object + auto rad = std::make_unique( + 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 diff --git a/src/oneD/Radiation1D.cpp b/src/oneD/Radiation1D.cpp index f9a5d2717..eecf07960 100644 --- a/src/oneD/Radiation1D.cpp +++ b/src/oneD/Radiation1D.cpp @@ -12,38 +12,25 @@ namespace Cantera { - -Radiation1D::Radiation1D(ThermoPhase* thermo, double pressure, size_t points, - std::function temperatureFunction, - std::function 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(); // 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 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 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 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 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& 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(), m_PMAC[sp_name]["coefficients"].asVector(), 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(), 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& coefficients, double temperature) { +double TabularPlanckMean::calculatePolynomial(const std::vector& 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(n)); @@ -143,18 +109,133 @@ double Radiation1D::calculatePolynomial(const std::vector& coefficients, return result / (1.0 * OneAtm); } -double Radiation1D::interpolateTable(const std::vector& temperatures, const std::vector& 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& temperatures, + const std::vector& 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& kabs, + std::vector& 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(), + m_PMAC[sp_name]["coefficients"].asVector(), + comp.T); + k_P += comp.P * x_sp * kVal; + } else if (fit_type == "polynomial") { + double kVal = calculatePolynomial( + m_PMAC[sp_name]["coefficients"].asVector(), + 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 temperatureFunction, + std::function moleFractionFunction, + std::unique_ptr props, + std::unique_ptr 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& 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 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 \ No newline at end of file