From 59cc717837147c6458cc81f916415f1aec2d39cc Mon Sep 17 00:00:00 2001 From: Christopher Neal Date: Mon, 21 Oct 2024 03:07:05 -0400 Subject: [PATCH 1/5] [oneD] reproducing PR965 with yaml file addition --- include/cantera/oneD/Flow1D.h | 26 ++++ src/oneD/Flow1D.cpp | 115 +++++++++++--- test/data/air_radiation.yaml | 274 ++++++++++++++++++++++++++++++++++ 3 files changed, 393 insertions(+), 22 deletions(-) create mode 100644 test/data/air_radiation.yaml diff --git a/include/cantera/oneD/Flow1D.h b/include/cantera/oneD/Flow1D.h index 17c158f9df..dfc1e00773 100644 --- a/include/cantera/oneD/Flow1D.h +++ b/include/cantera/oneD/Flow1D.h @@ -465,6 +465,8 @@ class Flow1D : public Domain1D * Computes the radiative heat loss vector over points jmin to jmax and stores * the data in the qdotRadiation variable. * + * The `fit-type` of `polynomial` is uses the model described below. + * * The simple radiation model used was established by Liu and Rogg * @cite liu1991. This model considers the radiation of CO2 and H2O. * @@ -475,6 +477,20 @@ class Flow1D : public Domain1D * 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. + * + * + * */ void computeRadiation(double* x, size_t jmin, size_t jmax); @@ -952,6 +968,16 @@ class Flow1D : public Domain1D //! radiative heat loss vector vector m_qdotRadiation; + // boundary emissivities for the radiation calculations + double m_epsilon_left = 0.0; + double m_epsilon_right = 0.0; + + std::map m_absorptionSpecies; //!< Absorbing species + AnyMap m_PMAC; //!< Absorption coefficient data for each species + + // Old radiation variable that can not be deleted for some reason + std::vector m_kRadiating; + // fixed T and Y values //! Fixed values of the temperature at each grid point that are used when solving //! with the energy equation disabled. diff --git a/src/oneD/Flow1D.cpp b/src/oneD/Flow1D.cpp index 4f07a6986c..f1609a64fa 100644 --- a/src/oneD/Flow1D.cpp +++ b/src/oneD/Flow1D.cpp @@ -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,52 @@ 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"); + // Parse radiation data from the YAML file + for (auto& name : m_thermo->speciesNames()) { + auto& data = m_thermo->species(name)->input; + if (data.hasKey("radiation")) { + cout << "Radiation data found for species " << name << 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")) { + cout << "Storing temperatures for species " << name << endl; + // Each species may have a specific set of temperatures that are used + m_PMAC[name]["temperatures"] = data["radiation"]["temperatures"].asVector(); + } else { + throw InputFileError("Flow1D::Flow1D", data, + "No 'temperatures' entry found for species '{}'", name); + } + if (data["radiation"].hasKey("data")) { + cout << "Storing data for species " << name << endl; + // This data is the Plank mean absorption coefficient + m_PMAC[name]["coefficients"] = data["radiation"]["data"].asVector(); + } else { + throw InputFileError("Flow1D::Flow1D", data, + "No 'data' entry found for species '{}'", name); + } + } else if (data["radiation"]["fit-type"] == "polynomial") { + cout << "Polynomial fit found for species " << name << endl; + if (data["radiation"].hasKey("data")) { + cout << "Storing data for species " << name << endl; + m_PMAC[name]["coefficients"] = data["radiation"]["data"].asVector(); + } 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); + } + } + } } Flow1D::Flow1D(shared_ptr th, size_t nsp, size_t points) @@ -467,34 +511,61 @@ void Flow1D::computeRadiation(double* x, size_t jmin, size_t jmax) // Polynomial coefficients: const double c_H2O[6] = {-0.23093, -1.12390, 9.41530, -2.99880, - 0.51382, -1.86840e-5}; + 0.51382, -1.86840e-5}; const double c_CO2[6] = {18.741, -121.310, 273.500, -194.050, - 56.310, -5.8169}; + 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); + double coef = 0.0; 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); + + for(const auto& [sp_name, sp_idx] : m_absorptionSpecies) { + if (m_PMAC[sp_name]["fit-type"].asString() == "table") { + // temperature table interval search + int T_index = 0; + const int OPL_table_size = m_PMAC[sp_name]["temperatures"].asVector().size(); + for (int k = 0; k < OPL_table_size; k++) { + if (T(x, j) < m_PMAC[sp_name]["temperatures"].asVector()[k]) { + if (T(x, j) < m_PMAC[sp_name]["temperatures"].asVector()[0]) { + T_index = 0; //lower table limit + } + else { + T_index = k; + } + break; + } + else { + T_index=OPL_table_size-1; //upper table limit + } + } + + // absorption coefficient for specie + double k_P_specie = 0.0; + if ((T_index == 0) || (T_index == OPL_table_size-1)) { + coef=log(1.0/m_PMAC[sp_name]["coefficients"].asVector()[T_index]); + } + else { + coef=log(1.0/m_PMAC[sp_name]["coefficients"].asVector()[T_index-1])+ + (log(1.0/m_PMAC[sp_name]["coefficients"].asVector()[T_index])-log(1.0/m_PMAC[sp_name]["coefficients"].asVector()[T_index-1]))* + (T(x, j)-m_PMAC[sp_name]["temperatures"].asVector()[T_index-1])/(m_PMAC[sp_name]["temperatures"].asVector()[T_index]-m_PMAC[sp_name]["temperatures"].asVector()[T_index-1]); + } + k_P_specie = exp(coef); + + k_P_specie /= k_P_ref; + k_P += m_press * X(x, m_absorptionSpecies[sp_name], j) * k_P_specie; + } else if (m_PMAC[sp_name]["fit-type"].asString() == "polynomial") { + double k_P_specie = 0.0; + for (size_t n = 0; n <= 5; n++) { + k_P_specie += m_PMAC[sp_name]["coefficients"].asVector()[n] * pow(1000 / T(x, j), (double) n); + } + k_P_specie /= k_P_ref; + k_P += m_press * X(x, m_absorptionSpecies[sp_name], j) * k_P_specie; } - 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 diff --git a/test/data/air_radiation.yaml b/test/data/air_radiation.yaml new file mode 100644 index 0000000000..6be8f62be6 --- /dev/null +++ b/test/data/air_radiation.yaml @@ -0,0 +1,274 @@ +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}} + +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} From 7cccceda77d6d546f2f85140bd8b83b5f4a21daa Mon Sep 17 00:00:00 2001 From: Christopher Neal Date: Tue, 22 Oct 2024 00:09:17 -0400 Subject: [PATCH 2/5] [oneD] added default PMAC values to maintain backwards compatibility --- src/oneD/Flow1D.cpp | 23 +++++++++++++++++------ 1 file changed, 17 insertions(+), 6 deletions(-) diff --git a/src/oneD/Flow1D.cpp b/src/oneD/Flow1D.cpp index f1609a64fa..869aece6b2 100644 --- a/src/oneD/Flow1D.cpp +++ b/src/oneD/Flow1D.cpp @@ -129,6 +129,23 @@ Flow1D::Flow1D(ThermoPhase* ph, size_t nsp, size_t points) : } } } + + // 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; + } } Flow1D::Flow1D(shared_ptr th, size_t nsp, size_t points) @@ -509,12 +526,6 @@ void Flow1D::computeRadiation(double* x, size_t jmin, size_t jmax) // 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); From f4aabb8cc2595c775f41d0c28f931a7ca5e8ddcb Mon Sep 17 00:00:00 2001 From: Christopher Neal Date: Fri, 10 Jan 2025 01:29:07 -0500 Subject: [PATCH 3/5] initial work on Radiation1D class --- include/cantera/oneD/Flow1D.h | 54 ++-------- include/cantera/oneD/Radiation1D.h | 105 ++++++++++++++++++ src/oneD/Flow1D.cpp | 145 ++----------------------- src/oneD/Radiation1D.cpp | 160 ++++++++++++++++++++++++++++ test/data/air_radiation.yaml | 2 + test/data/radiation-parameters.yaml | 55 ++++++++++ 6 files changed, 337 insertions(+), 184 deletions(-) create mode 100644 include/cantera/oneD/Radiation1D.h create mode 100644 src/oneD/Radiation1D.cpp create mode 100644 test/data/radiation-parameters.yaml diff --git a/include/cantera/oneD/Flow1D.h b/include/cantera/oneD/Flow1D.h index dfc1e00773..c1d2316e8a 100644 --- a/include/cantera/oneD/Flow1D.h +++ b/include/cantera/oneD/Flow1D.h @@ -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 @@ class Flow1D : public Domain1D //! 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,38 +463,8 @@ class Flow1D : public Domain1D //! 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 `fit-type` of `polynomial` is uses the model described below. - * - * 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. - * - * - * 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. - * - * - * - */ - 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); //! @} @@ -918,6 +890,9 @@ class Flow1D : public Domain1D //! radiative heat loss. double m_epsilon_right = 0.0; + //! Radiation object used for calculating radiative heat loss + std::unique_ptr m_radiation; + //! Indices within the ThermoPhase of the radiating species. First index is //! for CO2, second is for H2O. vector m_kRadiating; @@ -968,16 +943,6 @@ class Flow1D : public Domain1D //! radiative heat loss vector vector m_qdotRadiation; - // boundary emissivities for the radiation calculations - double m_epsilon_left = 0.0; - double m_epsilon_right = 0.0; - - std::map m_absorptionSpecies; //!< Absorbing species - AnyMap m_PMAC; //!< Absorption coefficient data for each species - - // Old radiation variable that can not be deleted for some reason - std::vector m_kRadiating; - // fixed T and Y values //! Fixed values of the temperature at each grid point that are used when solving //! with the energy equation disabled. @@ -1032,5 +997,4 @@ class Flow1D : public Domain1D }; } - #endif diff --git a/include/cantera/oneD/Radiation1D.h b/include/cantera/oneD/Radiation1D.h new file mode 100644 index 0000000000..6b7ebd368f --- /dev/null +++ b/include/cantera/oneD/Radiation1D.h @@ -0,0 +1,105 @@ +//! @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 + +namespace Cantera +{ + +/** + * Computes the radiative heat loss vector over points jmin to jmax and stores + * the data in the qdotRadiation variable. + * + * The `fit-type` of `polynomial` is uses the model described below. + * + * 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. + * + * + * 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 Radiation1D { +public: + Radiation1D(ThermoPhase* thermo, double pressure, size_t points, + std::function temperatureFunction, + std::function moleFractionFunction); + + // Parse radiation data from YAML input + void parseRadiationData(); + + // Compute radiative heat loss + void computeRadiation(double* x, size_t jmin, size_t jmax, vector& qdotRadiation); + + //! 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. + */ + 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; + double m_press; + size_t m_points; + + map m_absorptionSpecies; //!< Absorbing species + AnyMap m_PMAC; //!< Absorption coefficient data for each species + + //! Emissivity of the surface to the left 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; + + //! Lambda function to get mole fraction at a given point + std::function m_X; +}; + +} // namespace Cantera + +#endif // RADIATION1D_H diff --git a/src/oneD/Flow1D.cpp b/src/oneD/Flow1D.cpp index 869aece6b2..abdab89dcb 100644 --- a/src/oneD/Flow1D.cpp +++ b/src/oneD/Flow1D.cpp @@ -20,7 +20,10 @@ namespace Cantera Flow1D::Flow1D(ThermoPhase* ph, size_t nsp, size_t points) : Domain1D(nsp+c_offset_Y, points), - m_nsp(nsp) + m_nsp(nsp), + m_radiation(make_unique(ph, ph->pressure(), points, + [this](const double* x, size_t j) {return this->T(x,j);}, + [this](const double* x, size_t k, size_t j) {return this->X(x,k,j);})) { m_points = points; @@ -82,70 +85,6 @@ Flow1D::Flow1D(ThermoPhase* ph, size_t nsp, size_t points) : gr.push_back(1.0*ng/m_points); } setupGrid(m_points, gr.data()); - - // Parse radiation data from the YAML file - for (auto& name : m_thermo->speciesNames()) { - auto& data = m_thermo->species(name)->input; - if (data.hasKey("radiation")) { - cout << "Radiation data found for species " << name << 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")) { - cout << "Storing temperatures for species " << name << endl; - // Each species may have a specific set of temperatures that are used - m_PMAC[name]["temperatures"] = data["radiation"]["temperatures"].asVector(); - } else { - throw InputFileError("Flow1D::Flow1D", data, - "No 'temperatures' entry found for species '{}'", name); - } - if (data["radiation"].hasKey("data")) { - cout << "Storing data for species " << name << endl; - // This data is the Plank mean absorption coefficient - m_PMAC[name]["coefficients"] = data["radiation"]["data"].asVector(); - } else { - throw InputFileError("Flow1D::Flow1D", data, - "No 'data' entry found for species '{}'", name); - } - } else if (data["radiation"]["fit-type"] == "polynomial") { - cout << "Polynomial fit found for species " << name << endl; - if (data["radiation"].hasKey("data")) { - cout << "Storing data for species " << name << endl; - m_PMAC[name]["coefficients"] = data["radiation"]["data"].asVector(); - } 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 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; - } } Flow1D::Flow1D(shared_ptr th, size_t nsp, size_t points) @@ -522,70 +461,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; - - // 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); - - double coef = 0.0; - for (size_t j = jmin; j < jmax; j++) { - // calculation of the mean Planck absorption coefficient - double k_P = 0; - - for(const auto& [sp_name, sp_idx] : m_absorptionSpecies) { - if (m_PMAC[sp_name]["fit-type"].asString() == "table") { - // temperature table interval search - int T_index = 0; - const int OPL_table_size = m_PMAC[sp_name]["temperatures"].asVector().size(); - for (int k = 0; k < OPL_table_size; k++) { - if (T(x, j) < m_PMAC[sp_name]["temperatures"].asVector()[k]) { - if (T(x, j) < m_PMAC[sp_name]["temperatures"].asVector()[0]) { - T_index = 0; //lower table limit - } - else { - T_index = k; - } - break; - } - else { - T_index=OPL_table_size-1; //upper table limit - } - } - - // absorption coefficient for specie - double k_P_specie = 0.0; - if ((T_index == 0) || (T_index == OPL_table_size-1)) { - coef=log(1.0/m_PMAC[sp_name]["coefficients"].asVector()[T_index]); - } - else { - coef=log(1.0/m_PMAC[sp_name]["coefficients"].asVector()[T_index-1])+ - (log(1.0/m_PMAC[sp_name]["coefficients"].asVector()[T_index])-log(1.0/m_PMAC[sp_name]["coefficients"].asVector()[T_index-1]))* - (T(x, j)-m_PMAC[sp_name]["temperatures"].asVector()[T_index-1])/(m_PMAC[sp_name]["temperatures"].asVector()[T_index]-m_PMAC[sp_name]["temperatures"].asVector()[T_index-1]); - } - k_P_specie = exp(coef); - - k_P_specie /= k_P_ref; - k_P += m_press * X(x, m_absorptionSpecies[sp_name], j) * k_P_specie; - } else if (m_PMAC[sp_name]["fit-type"].asString() == "polynomial") { - double k_P_specie = 0.0; - for (size_t n = 0; n <= 5; n++) { - k_P_specie += m_PMAC[sp_name]["coefficients"].asVector()[n] * pow(1000 / T(x, j), (double) n); - } - k_P_specie /= k_P_ref; - k_P += m_press * X(x, m_absorptionSpecies[sp_name], j) * k_P_specie; - } - } - - // 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, @@ -1175,16 +1051,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) diff --git a/src/oneD/Radiation1D.cpp b/src/oneD/Radiation1D.cpp new file mode 100644 index 0000000000..f9a5d27177 --- /dev/null +++ b/src/oneD/Radiation1D.cpp @@ -0,0 +1,160 @@ +//! @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 +{ + + +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) +{ + parseRadiationData(); +} + +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::parseRadiationData() { + AnyMap radiationPropertiesDB; + // Search 'crit-properties.yaml' to find Tc and Pc. Load data if needed. + if (radiationPropertiesDB.empty()) { + radiationPropertiesDB = AnyMap::fromYamlFile("radiation-properties.yaml"); + } + + if( radiationPropertiesDB.hasKey("PMAC")) { + auto& data = radiationPropertiesDB["PMAC"].as(); + + // 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(); + } 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(); + } 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(); + } 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 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; + } + } +} + +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 result = 0.0; + for (size_t n = 0; n < coefficients.size(); ++n) { + result += coefficients[n] * std::pow(1000 / temperature, static_cast(n)); + } + 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; + 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); +} + + +} // namespace Cantera \ No newline at end of file diff --git a/test/data/air_radiation.yaml b/test/data/air_radiation.yaml index 6be8f62be6..e408c08971 100644 --- a/test/data/air_radiation.yaml +++ b/test/data/air_radiation.yaml @@ -17,6 +17,8 @@ phases: 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 diff --git a/test/data/radiation-parameters.yaml b/test/data/radiation-parameters.yaml new file mode 100644 index 0000000000..56ba6e1523 --- /dev/null +++ b/test/data/radiation-parameters.yaml @@ -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] + From 3b8c0ce59448f5481b7d86a21f07c1082c5b586b Mon Sep 17 00:00:00 2001 From: Christopher Neal Date: Tue, 28 Jan 2025 02:31:16 -0500 Subject: [PATCH 4/5] [oneD] use a factory-like initialization for radiation object help by Flow1D --- src/oneD/Flow1D.cpp | 31 +++++++++++++++++++++++++++---- 1 file changed, 27 insertions(+), 4 deletions(-) diff --git a/src/oneD/Flow1D.cpp b/src/oneD/Flow1D.cpp index abdab89dcb..2aa0f938d0 100644 --- a/src/oneD/Flow1D.cpp +++ b/src/oneD/Flow1D.cpp @@ -20,10 +20,7 @@ namespace Cantera Flow1D::Flow1D(ThermoPhase* ph, size_t nsp, size_t points) : Domain1D(nsp+c_offset_Y, points), - m_nsp(nsp), - m_radiation(make_unique(ph, ph->pressure(), points, - [this](const double* x, size_t j) {return this->T(x,j);}, - [this](const double* x, size_t k, size_t j) {return this->X(x,k,j);})) + m_nsp(nsp) { m_points = points; @@ -85,6 +82,32 @@ Flow1D::Flow1D(ThermoPhase* ph, size_t nsp, size_t points) : gr.push_back(1.0*ng/m_points); } setupGrid(m_points, gr.data()); + + // 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 th, size_t nsp, size_t points) From cb1aec990cb6fadb66d0bae8c986f31300d7716c Mon Sep 17 00:00:00 2001 From: Christopher Neal Date: Tue, 28 Jan 2025 02:32:14 -0500 Subject: [PATCH 5/5] updates to the radiation classes and starting to add documentation --- include/cantera/oneD/Radiation1D.h | 388 ++++++++++++++++++++++++++--- src/oneD/Radiation1D.cpp | 219 +++++++++++----- 2 files changed, 500 insertions(+), 107 deletions(-) diff --git a/include/cantera/oneD/Radiation1D.h b/include/cantera/oneD/Radiation1D.h index 6b7ebd368f..c732739116 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. + * + * 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. * - * 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 `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 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& 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 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); + + /** 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::function moleFractionFunction, + std::unique_ptr props, + std::unique_ptr solver); - // Parse radiation data from YAML input - void parseRadiationData(); - - // Compute radiative heat loss - void computeRadiation(double* x, size_t jmin, size_t jmax, vector& qdotRadiation); + /** 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); - //! 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. + /** + * 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 @@ class Radiation1D { } 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 - //! radiative heat loss. - double m_epsilon_left = 0.0; + //! Solver for radiative heat loss + std::unique_ptr m_solver; - //! Emissivity of the surface to the right of the domain. Used for calculating + //! 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; - // 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 @@ class Radiation1D { 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 f9a5d27177..eecf079607 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; - } - - // 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); + // 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; + } - 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); + // 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; } } -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