Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Revive 1D radiation improvements of PR965 #1799

Open
wants to merge 2 commits into
base: main
Choose a base branch
from
Open

Revive 1D radiation improvements of PR965 #1799

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
a9bf1e5
[oneD] reproducing PR965 with yaml file addition
wandadars Oct 21, 2024
336a341
[oneD] added default PMAC values to maintain backwards compatibility
wandadars Oct 22, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
34 changes: 26 additions & 8 deletions include/cantera/oneD/Flow1D.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -447,6 +447,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.
*
Expand All @@ -457,6 +459,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);

Expand Down Expand Up @@ -853,14 +869,6 @@ class Flow1D : public Domain1D
Kinetics* m_kin = nullptr;
Transport* m_trans = nullptr;

// boundary emissivities for the radiation calculations
double m_epsilon_left = 0.0;
double m_epsilon_right = 0.0;

//! Indices within the ThermoPhase of the radiating species. First index is
//! for CO2, second is for H2O.
vector<size_t> m_kRadiating;

// flags
vector<bool> m_do_energy;
bool m_do_soret = false;
Expand All @@ -874,6 +882,16 @@ class Flow1D : public Domain1D
//! radiative heat loss vector
vector<double> m_qdotRadiation;

// boundary emissivities for the radiation calculations
double m_epsilon_left = 0.0;
double m_epsilon_right = 0.0;

std::map<std::string, int> 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<size_t> m_kRadiating;

// fixed T and Y values
vector<double> m_fixedtemp;
vector<double> m_zfix;
Expand Down
134 changes: 108 additions & 26 deletions src/oneD/Flow1D.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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;

Expand Down Expand Up @@ -84,10 +86,69 @@ 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<double>();
} 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<double>();
} 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<double>();
} else {
throw InputFileError("Flow1D::Flow1D", data,
"No 'data' entry found for species '{}'", name);
}
} else {
throw InputFileError("Flow1D::Flow1D", data,
"Invalid 'fit-type' entry found for species '{}'", name);
}
}
}

// Polynomial coefficients for CO2 and H2O (backwards compatibility)
// Check if "CO2" is already in the map, if not, add the polynomial fit data
if (!m_PMAC.hasKey("CO2")) {
const std::vector<double> c_CO2 = {18.741, -121.310, 273.500, -194.050, 56.310,
-5.8169};
m_PMAC["CO2"]["fit-type"] = "polynomial";
m_PMAC["CO2"]["coefficients"] = c_CO2;
}

// Check if "H2O" is already in the map, if not, add the polynomial fit data
if (!m_PMAC.hasKey("H2O")) {
const std::vector<double> c_H2O = {-0.23093, -1.12390, 9.41530, -2.99880,
0.51382, -1.86840e-5};
m_PMAC["H2O"]["fit-type"] = "polynomial";
m_PMAC["H2O"]["coefficients"] = c_H2O;
}
}

Flow1D::Flow1D(shared_ptr<ThermoPhase> th, size_t nsp, size_t points)
Expand Down Expand Up @@ -468,36 +529,57 @@ 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);

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<double>().size();
for (int k = 0; k < OPL_table_size; k++) {
if (T(x, j) < m_PMAC[sp_name]["temperatures"].asVector<double>()[k]) {
if (T(x, j) < m_PMAC[sp_name]["temperatures"].asVector<double>()[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<double>()[T_index]);
}
else {
coef=log(1.0/m_PMAC[sp_name]["coefficients"].asVector<double>()[T_index-1])+
(log(1.0/m_PMAC[sp_name]["coefficients"].asVector<double>()[T_index])-log(1.0/m_PMAC[sp_name]["coefficients"].asVector<double>()[T_index-1]))*
(T(x, j)-m_PMAC[sp_name]["temperatures"].asVector<double>()[T_index-1])/(m_PMAC[sp_name]["temperatures"].asVector<double>()[T_index]-m_PMAC[sp_name]["temperatures"].asVector<double>()[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<double>()[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
Expand Down
Loading
Loading