Skip to content

Commit

Permalink
initial work on Radiation1D class
Browse files Browse the repository at this point in the history
  • Loading branch information
@wandadars
wandadars committed Jan 13, 2025
1 parent 7cccced commit f4aabb8
Show file tree
Hide file tree
Showing 6 changed files with 337 additions and 184 deletions.
54 changes: 9 additions & 45 deletions include/cantera/oneD/Flow1D.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
{

Expand Down Expand Up @@ -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`.
Expand Down Expand Up @@ -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);

//! @}

Expand Down Expand Up @@ -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<Radiation1D> m_radiation;

//! Indices within the ThermoPhase of the radiating species. First index is
//! for CO2, second is for H2O.
vector<size_t> m_kRadiating;
Expand Down Expand Up @@ -968,16 +943,6 @@ 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
//! Fixed values of the temperature at each grid point that are used when solving
//! with the energy equation disabled.
Expand Down Expand Up @@ -1032,5 +997,4 @@ class Flow1D : public Domain1D
};

}

#endif
105 changes: 105 additions & 0 deletions include/cantera/oneD/Radiation1D.h
View file
Original file line number Diff line number Diff line change
@@ -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 <functional>

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<double(const double*, size_t)> temperatureFunction,
std::function<double(const double*, size_t, size_t)> moleFractionFunction);

// Parse radiation data from YAML input
void parseRadiationData();

// Compute radiative heat loss
void computeRadiation(double* x, size_t jmin, size_t jmax, vector<double>& 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<string, size_t> 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<double>& coefficients, double temperature);
double interpolateTable(const vector<double>& temperatures, const vector<double>& data, double temperature);

//! Lambda function to get temperature at a given point
std::function<double(const double*, size_t)> m_T;

//! Lambda function to get mole fraction at a given point
std::function<double(const double*, size_t, size_t)> m_X;
};

} // namespace Cantera

#endif // RADIATION1D_H
145 changes: 6 additions & 139 deletions src/oneD/Flow1D.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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<Radiation1D>(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;

Expand Down Expand Up @@ -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<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 @@ -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<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;
}
}

// 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,
Expand Down Expand Up @@ -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)
Expand Down
Loading

0 comments on commit f4aabb8

Please sign in to comment.