Skip to content

Radiation model update in StFlow.cpp #965

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

Closed
lavrenyukiv wants to merge 7 commits into from
Closed

Radiation model update in StFlow.cpp #965

Changes from 4 commits
Commits
Show all changes
7 commits
Select commit Hold shift + click to select a range
05c7925
Update StFlow.cpp
lavrenyukiv Jan 25, 2021
6201d6a
Update StFlow.cpp
lavrenyukiv Jan 25, 2021
234fa4b
Update StFlow.cpp
lavrenyukiv Jan 25, 2021
5db79c0
Update StFlow.cpp
lavrenyukiv Jan 26, 2021
0cfb26f
Update StFlow.cpp
lavrenyukiv Feb 1, 2021
ea5d074
Update StFlow.h
lavrenyukiv Feb 1, 2021
89cc361
Update StFlow.cpp
lavrenyukiv Feb 1, 2021
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
168 changes: 153 additions & 15 deletions src/oneD/StFlow.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -94,9 +94,12 @@ StFlow::StFlow(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.resize(5, npos);

@ischoegl ischoegl Jan 26, 2021 •
edited
Loading

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

There has been some work to make the identification of species more flexible (which would mean that lowercase species names would be recognized). Rather than hard coding names in the subsequent lines, specifying the elemental composition will in most cases work (there usually is only one isomer). The function name is findIsomers

@lavrenyukiv lavrenyukiv Feb 3, 2021

Copy link
Copy Markdown
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I changed this part.

It would be great if the user could select the components involved in the radiation himself.

It would probably be better to load the absorption coefficient data from a separate file, so that the user can add his own dependencies.

Unfortunately, my experience with С++ is not enough.

@ischoegl ischoegl Feb 3, 2021

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I agree that it would be great to select species / import parameters instead of hard coding them. There are likely multiple ways of getting this done. I can probably give you pointers, but it’s up to @bryanwweber and @speth ...

m_kRadiating[0] = m_thermo->speciesIndex("CO2");
m_kRadiating[1] = m_thermo->speciesIndex("H2O");
m_kRadiating[2] = m_thermo->speciesIndex("CO");
m_kRadiating[3] = m_thermo->speciesIndex("CH4");
m_kRadiating[4] = m_thermo->speciesIndex("OH");
}

void StFlow::resize(size_t ncomponents, size_t points)
Expand Down Expand Up @@ -313,18 +316,56 @@ void StFlow::evalResidual(double* x, double* rsd, int* diag,
if (m_do_radiation) {
// variable definitions for the Planck absorption coefficient and the
// radiation calculation:
doublereal k_P_ref = 1.0*OneAtm;

// polynomial coefficients:
const doublereal c_H2O[6] = {-0.23093, -1.12390, 9.41530, -2.99880,
0.51382, -1.86840e-5};
const doublereal c_CO2[6] = {18.741, -121.310, 273.500, -194.050,
56.310, -5.8169};

double k_P_ref = 1.0*OneAtm;

// Temperatures for Planck optical path length evaluation, K
const int OPL_table_size = 28;
const double TemperatureOPL[OPL_table_size] = {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};

// Planck optical path length, m
const double PlanckOPL_CO2[OPL_table_size] = {0.0385, 0.039, 0.033, 0.0303,
0.0307, 0.0335, 0.0381, 0.0447, 0.0533, 0.0645, 0.0787, 0.0966, 0.119,
0.147, 0.183, 0.227, 0.281, 0.348, 0.431, 0.532, 0.655, 0.803, 0.983,
1.198, 1.455, 1.761, 2.123, 2.55};
const double PlanckOPL_H2O[OPL_table_size] = {0.0194, 0.0349, 0.0523, 0.0728,
0.097, 0.125, 0.158, 0.196, 0.24, 0.29, 0.352, 0.424, 0.51, 0.614, 0.737,
0.885, 1.06, 1.27, 1.52, 1.82, 2.18, 2.6, 3.1, 3.69, 4.38, 5.2, 6.15, 7.28};
const double PlanckOPL_CH4[OPL_table_size] = {0.00936, 0.0128, 0.0189, 0.029,
0.0458, 0.0729, 0.1166, 0.186, 0.296, 0.467, 0.732, 1.136, 1.746, 2.655,
3.999, 5.96, 8.8, 12.87, 18.64, 26.75, 38.07, 53.71, 75.17, 104.37,
143.83, 196.8, 267.3, 360.8};
const double PlanckOPL_CO[OPL_table_size] = {0.0139, 0.024, 0.0439, 0.0785,
0.134, 0.22, 0.346, 0.53, 0.799, 1.189, 1.75, 2.57, 3.728, 5.372, 7.68,
10.87, 15.24, 21.15, 29.05, 39.46, 53.03, 70.49, 92.71, 120.67, 155.48,
198.39, 250.77, 314.1};
const double PlanckOPL_OH[OPL_table_size] = {0.00794, 0.0113, 0.01624, 0.024,
0.03635, 0.05585, 0.08625, 0.133, 0.2042, 0.3108, 0.4687, 0.6994, 1.0324,
1.507, 2.177, 3.111, 4.401, 6.163, 8.55, 11.75, 16.0, 21.62, 28.98, 38.54,
50.88, 66.69, 86.85, 112.37};

// natural logarithms of reversed optical path lengths for linear interpolation in logarithm scale
double OPL_log_CO2[OPL_table_size];
double OPL_log_H2O[OPL_table_size];
double OPL_log_CH4[OPL_table_size];
double OPL_log_CO[OPL_table_size];
double OPL_log_OH[OPL_table_size];

for (int j = 0; j < OPL_table_size; j++) {
OPL_log_CO2[j] = log(1.0/PlanckOPL_CO2[j]);
OPL_log_H2O[j] = log(1.0/PlanckOPL_H2O[j]);
OPL_log_CH4[j] = log(1.0/PlanckOPL_CH4[j]);
OPL_log_CO[j] = log(1.0/PlanckOPL_CO[j]);
OPL_log_OH[j] = log(1.0/PlanckOPL_OH[j]);
}

// 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;
// loop over all grid points
for (size_t j = jmin; j < jmax; j++) {
// helping variable for the calculation
Expand All @@ -335,22 +376,119 @@ void StFlow::evalResidual(double* x, double* rsd, int* diag,
// 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);
for (int k = 0; k < OPL_table_size; k++) {
if (T(x, j) < TemperatureOPL[k]) {

@ischoegl ischoegl Jan 27, 2021

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This implementation could be improved in terms of computational efficiency as it appears that the same temperature lookup loop is called for each of the species. Beyond, It should be possible to consolidate the species-specific code sections, which - as far as I can tell - are mostly identical except for the species index.

@ischoegl ischoegl Jan 28, 2021 •
edited
Loading

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

PS: I am also not convinced that a linear interpolation within a lookup table is an improvement over a polynomial fit - unless the approximation is poor for the newly added curves. (Btw, how do new H20/CO2 curves compare to the previous implementation?)

@lavrenyukiv lavrenyukiv Jan 28, 2021 •
edited
Loading

Copy link
Copy Markdown
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I will make a comparison between old and new curves for H2O and CO2.
Linear interpolation is not an improvement, but rather a more general assumption. For other species we have not polynomial type assumption and polynomial coefficients. But we have spectra for it. Integration of spectra is a very expensive operation in terms of calculations. So tabulated data is accurate and fast.
The use of a polynomial can lead to negative values of the absorption coefficients.
I found that CO, CH4 and OH molecules spectra integrated before and presented by me are wrong, so I will change this data.

@lavrenyukiv lavrenyukiv Jan 28, 2021 •
edited
Loading

Copy link
Copy Markdown
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Plots in logarithmic coordinates are at the left part, the same plots in linear coordinates are at the right part.
Results for H2O look similar. Most likely the polynomial approximation was done on a linear scale. So there is a large difference in the logarithmic scale at high temperatures. So linear interpolation in logarithmic scale looks preferable.
The results for CO2 are very different. I will try to find the data from the original source.
models_check

If we compute S * a * (T^4-300^4), where S is Stefan - Boltzmann constant, we can make an estimate of the absorbed energy(W/m3/atm). Here the differences between approximations will be huge.
models_check_2

@ischoegl ischoegl Jan 28, 2021 •
edited
Loading

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@lavrenyukiv ... thank you - this looks interesting. I was not involved in the original model, so it's hard to tell. I agree that the differences for energy are substantial. The original polynomial coefficients appear to be taken directly from the TNF link you had posted (where it states that "These curve fits were generated for temperatures between 300K and 2500K and may be very inaccurate outside this range."; the fit itself is for 1000/T so deviations may not look as dramatic?). I am not familiar with RADCAL and assumptions used therein, so some legitimate differences may exist.

PS: There appears to be an updated version of the radcal code on GitHub.

PPS: Incidentally, there was a post on the User Group about unphysically low temperatures but it's hard to tell whether these are related as no code was posted.

@lavrenyukiv lavrenyukiv Jan 28, 2021

Copy link
Copy Markdown
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes, it is a negative absorption coefficient at low temperatures.
models_check_3

@lavrenyukiv lavrenyukiv Jan 30, 2021 •
edited
Loading

Copy link
Copy Markdown
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I am trying to compare spectra.
Results for CO2 from https://github.com/firemodels/radcal/blob/master/Documentation/Figures/CO2_300K.png

image

and spectra from hapi

full range
co2_full

the RADCAL range
co2_crop

detailed view
co2_closest

Spectra look similar in the specified wavenumber range. The hapi has a greater resolution.

if (T(x, j) < TemperatureOPL[0]) {
coef=OPL_log_H2O[0];
}
else {
coef=OPL_log_H2O[k-1]+(OPL_log_H2O[k]-OPL_log_H2O[k-1])*
(T(x, j)-TemperatureOPL[k-1])/
(TemperatureOPL[k]-TemperatureOPL[k-1]);
}
break;
}
else {
coef=OPL_log_H2O[OPL_table_size-1];
}
}
k_P_H2O = exp(coef);
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 (int k = 0; k < OPL_table_size; k++) {
if (T(x, j) < TemperatureOPL[k]) {
if (T(x, j) < TemperatureOPL[0]) {
coef=OPL_log_CO2[0];
}
else {
coef=OPL_log_CO2[k-1]+(OPL_log_CO2[k]-OPL_log_CO2[k-1])*
(T(x, j)-TemperatureOPL[k-1])/
(TemperatureOPL[k]-TemperatureOPL[k-1]);
}
break;
}
else {
coef=OPL_log_CO2[OPL_table_size-1];
}
}
k_P_CO2 = exp(coef);
k_P_CO2 /= k_P_ref;
k_P += m_press * X(x, m_kRadiating[0], j) * k_P_CO2;
}

// absorption coefficient for CO
if (m_kRadiating[2] != npos) {
double k_P_CO = 0;
for (int k = 0; k < OPL_table_size; k++) {
if (T(x, j) < TemperatureOPL[k]) {
if (T(x, j) < TemperatureOPL[0]) {
coef=OPL_log_CO[0];
}
else {
coef=OPL_log_CO[k-1]+(OPL_log_CO[k]-OPL_log_CO[k-1])*
(T(x, j)-TemperatureOPL[k-1])/
(TemperatureOPL[k]-TemperatureOPL[k-1]);
}
break;
}
else {
coef=OPL_log_CO[OPL_table_size-1];
}
}
k_P_CO = exp(coef);
k_P_CO /= k_P_ref;
k_P += m_press * X(x, m_kRadiating[2], j) * k_P_CO;
}
// absorption coefficient for CH4
if (m_kRadiating[3] != npos) {
double k_P_CH4 = 0;
for (int k = 0; k < OPL_table_size; k++) {
if (T(x, j) < TemperatureOPL[k]) {
if (T(x, j) < TemperatureOPL[0]) {
coef=OPL_log_CH4[0];
}
else {
coef=OPL_log_CH4[k-1]+(OPL_log_CH4[k]-OPL_log_CH4[k-1])*
(T(x, j)-TemperatureOPL[k-1])/
(TemperatureOPL[k]-TemperatureOPL[k-1]);
}
break;
}
else {
coef=OPL_log_CH4[OPL_table_size-1];
}
}
k_P_CH4 = exp(coef);
k_P_CH4 /= k_P_ref;
k_P += m_press * X(x, m_kRadiating[3], j) * k_P_CH4;
}
// absorption coefficient for OH
if (m_kRadiating[4] != npos) {
double k_P_OH = 0;
for (int k = 0; k < OPL_table_size; k++) {
if (T(x, j) < TemperatureOPL[k]) {
if (T(x, j) < TemperatureOPL[0]) {
coef=OPL_log_OH[0];
}
else {
coef=OPL_log_OH[k-1]+(OPL_log_OH[k]-OPL_log_OH[k-1])*
(T(x, j)-TemperatureOPL[k-1])/
(TemperatureOPL[k]-TemperatureOPL[k-1]);
}
break;
}
else {
coef=OPL_log_OH[OPL_table_size-1];
}
}
k_P_OH = exp(coef);
k_P_OH /= k_P_ref;
k_P += m_press * X(x, m_kRadiating[4], j) * k_P_OH;
}
// calculation of the radiative heat loss term
radiative_heat_loss = 2 * k_P *(2 * StefanBoltz * pow(T(x, j), 4)
- boundary_Rad_left - boundary_Rad_right);
Expand Down