From 05c7925df7914b5c013346e38f9aad3f408b5120 Mon Sep 17 00:00:00 2001 From: lavrenyukiv Date: Mon, 25 Jan 2021 22:36:27 +0500 Subject: [PATCH 1/7] Update StFlow.cpp --- src/oneD/StFlow.cpp | 143 ++++++++++++++++++++++++++++++++++++++++---- 1 file changed, 130 insertions(+), 13 deletions(-) diff --git a/src/oneD/StFlow.cpp b/src/oneD/StFlow.cpp index df9b1e7d65..615c02c09c 100644 --- a/src/oneD/StFlow.cpp +++ b/src/oneD/StFlow.cpp @@ -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); 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) @@ -315,12 +318,32 @@ void StFlow::evalResidual(double* x, double* rsd, int* diag, // 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}; - + // Temperatures for Planck optical path length evaluation, K + const int OPL_table_size = 28; + const doublereal 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 doublereal 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 doublereal 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 doublereal 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 doublereal 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 doublereal PlanckOPL_OH[OPL_table_size] = {0.00794, 0.0113, .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 + doublereal OPL_log_CO2[OPL_table_size]; + doublereal OPL_log_H2O[OPL_table_size]; + doublereal OPL_log_CH4[OPL_table_size]; + doublereal OPL_log_CO[OPL_table_size]; + doublereal 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); @@ -335,21 +358,115 @@ 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]) { + 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[3] != 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) From 6201d6a2b1610f4d72522fed93453c43c05b03dc Mon Sep 17 00:00:00 2001 From: lavrenyukiv Date: Mon, 25 Jan 2021 23:13:30 +0500 Subject: [PATCH 2/7] Update StFlow.cpp index error fixed. --- src/oneD/StFlow.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/oneD/StFlow.cpp b/src/oneD/StFlow.cpp index 615c02c09c..32c7c2dd63 100644 --- a/src/oneD/StFlow.cpp +++ b/src/oneD/StFlow.cpp @@ -446,7 +446,7 @@ void StFlow::evalResidual(double* x, double* rsd, int* diag, } // absorption coefficient for OH - if (m_kRadiating[3] != npos) { + 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]) { From 234fa4b64bf3ac778d4655da872ca3793e7bb08a Mon Sep 17 00:00:00 2001 From: lavrenyukiv Date: Mon, 25 Jan 2021 23:26:36 +0500 Subject: [PATCH 3/7] Update StFlow.cpp error fixed. --- src/oneD/StFlow.cpp | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/oneD/StFlow.cpp b/src/oneD/StFlow.cpp index 32c7c2dd63..51c357a94f 100644 --- a/src/oneD/StFlow.cpp +++ b/src/oneD/StFlow.cpp @@ -347,8 +347,10 @@ void StFlow::evalResidual(double* x, double* rsd, int* diag, // 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 + // loop over all grid points for (size_t j = jmin; j < jmax; j++) { // helping variable for the calculation double radiative_heat_loss = 0; From 5db79c0d54cd447541c9a2d4cb38567e261c962e Mon Sep 17 00:00:00 2001 From: lavrenyukiv Date: Tue, 26 Jan 2021 23:30:55 +0500 Subject: [PATCH 4/7] Update StFlow.cpp - replace tabs with spaces - replace doublereal types with double --- src/oneD/StFlow.cpp | 217 ++++++++++++++++++++++++-------------------- 1 file changed, 118 insertions(+), 99 deletions(-) diff --git a/src/oneD/StFlow.cpp b/src/oneD/StFlow.cpp index 51c357a94f..0bc0ab3ab0 100644 --- a/src/oneD/StFlow.cpp +++ b/src/oneD/StFlow.cpp @@ -316,41 +316,57 @@ 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; + double k_P_ref = 1.0*OneAtm; // Temperatures for Planck optical path length evaluation, K const int OPL_table_size = 28; - const doublereal 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}; + 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 doublereal 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 doublereal 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 doublereal 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 doublereal 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 doublereal PlanckOPL_OH[OPL_table_size] = {0.00794, 0.0113, .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}; + 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 - doublereal OPL_log_CO2[OPL_table_size]; - doublereal OPL_log_H2O[OPL_table_size]; - doublereal OPL_log_CH4[OPL_table_size]; - doublereal OPL_log_CO[OPL_table_size]; - doublereal OPL_log_OH[OPL_table_size]; + 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]); - } + 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 + double coef; + // loop over all grid points for (size_t j = jmin; j < jmax; j++) { // helping variable for the calculation double radiative_heat_loss = 0; @@ -361,21 +377,22 @@ void StFlow::evalResidual(double* x, double* rsd, int* diag, if (m_kRadiating[1] != npos) { double k_P_H2O = 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_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); - + if (T(x, j) < TemperatureOPL[k]) { + 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; } @@ -384,21 +401,22 @@ void StFlow::evalResidual(double* x, double* rsd, int* diag, double k_P_CO2 = 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_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); - + 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; } @@ -406,21 +424,22 @@ void StFlow::evalResidual(double* x, double* rsd, int* diag, 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); - + 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; } @@ -428,48 +447,48 @@ void StFlow::evalResidual(double* x, double* rsd, int* diag, 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); - + 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); - + 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); From 0cfb26f2af1e02112c29550bedce2838026aa6f6 Mon Sep 17 00:00:00 2001 From: lavrenyukiv Date: Tue, 2 Feb 2021 03:15:25 +0500 Subject: [PATCH 5/7] Update StFlow.cpp - Added coefficients for CO2, H2O, CO, CH4, OH, N2O, NH3, H2O2, H2CO, NO2, NO, HCN, C2H2, C2H4, C2H6, C2N2, C4H2, HCOOH, HNO3, O3 - Temprature range expanded to 200 - 3500 K --- src/oneD/StFlow.cpp | 414 +++++++++++++++++++++++++++----------------- 1 file changed, 256 insertions(+), 158 deletions(-) diff --git a/src/oneD/StFlow.cpp b/src/oneD/StFlow.cpp index 0bc0ab3ab0..a263b70447 100644 --- a/src/oneD/StFlow.cpp +++ b/src/oneD/StFlow.cpp @@ -94,12 +94,14 @@ StFlow::StFlow(ThermoPhase* ph, size_t nsp, size_t points) : setupGrid(m_points, gr.data()); // Find indices for radiating species - m_kRadiating.resize(5, npos); - 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"); + AbsorptionSpeciesList = {"CO2","H2O","CO","CH4","OH","N2O","NH3", + "H2O2","H2CO","NO2","NO","HCN","C2H2","C2H4","C2H6","C2N2", + "C4H2","HCOOH","HNO3","O3"}; + //AbsorptionSpeciesList = {"H2O"}; + + for(std::string spIter : AbsorptionSpeciesList) { + AbsorptionSpeciesMap.insert({spIter, m_thermo->speciesIndex(spIter)}); + } } void StFlow::resize(size_t ncomponents, size_t points) @@ -316,179 +318,275 @@ void StFlow::evalResidual(double* x, double* rsd, int* diag, if (m_do_radiation) { // variable definitions for the Planck absorption coefficient and the // radiation calculation: - double k_P_ref = 1.0*OneAtm; - + doublereal 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}; + std::vector TemperatureOPL = {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}; + + const int OPL_table_size = TemperatureOPL.size(); - // 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]); - } + // Planck optical path length, m + std::vector PlanckOPL_CO2 = {0.028819231725166233, 0.03851088876725997, + 0.03899126980906671, 0.032980228203593136, 0.030267348946179023, 0.030743052568319492, + 0.03349552664432657, 0.03815025367862556, 0.044695166534542856, 0.0533451125340091, + 0.06448815851162609, 0.07867018394357354, 0.09659983965382253, 0.11916430710589461, + 0.14745252150882412, 0.18278437065235484, 0.22674394556072777, 0.2812154755687562, + 0.34842470656989205, 0.43098095501152517, 0.5319261720083194, 0.6547800464563425, + 0.8035960044763693, 0.983014582040391, 1.1983144211620795, 1.4554774859486679, + 1.7612497326430492, 2.123189910047587, 2.549740894200347, 3.050301729813575, + 3.6352339058247307, 4.316015725482628, 5.1051815487411565, 6.016472443851096}; + std::vector PlanckOPL_H2O = {0.007180563211295323, 0.019365261372768307, + 0.034867966945125076, 0.05230634537613443, 0.07276548153611018, 0.09699925744256371, + 0.12530329568293222, 0.15801661948972606, 0.19580389577973367, 0.23972350297492273, + 0.2912064962199706, 0.35202916571833975, 0.4243127637866092, 0.510551491576356, + 0.6136629518326753, 0.7370547168162068, 0.8847050123750548, 1.061254758533107, + 1.2721110105236084, 1.5235621242734583, 1.822906018849437, 2.178591373741365, + 2.6003745205523052, 3.099489796041002, 3.688832242804873, 4.383157080523857, + 5.1993025902907215, 6.156408948413962, 7.276171359454505, 8.583098608019881, + 10.104757796299838, 11.87209392358693, 13.919670784850705, 16.28600423253945}; + std::vector PlanckOPL_CO = {30.32780456023482, 1.5025772744944488, + 0.48568677330878235, 0.31675988274374905, 0.2809201570439709, 0.28946895791957333, + 0.32239055951412515, 0.37417537060525935, 0.4437961211963581, 0.5320760473126184, + 0.6408278564965408, 0.7725511001710472, 0.9303212767894896, 1.1177046866345075, + 1.3387967943387578, 1.598138988644575, 1.9007086226166436, 2.252003287440127, + 2.6578585623777644, 3.1245488139145126, 3.658747887146305, 4.2674749977544435, + 4.958101084203171, 5.7383843347079315, 6.616293182716349, 7.600244912979516, + 8.698845819200764, 9.920965846773354, 11.276052909388223, 12.77335136618167, + 14.42279395291501, 16.234355340729095, 18.218430714508237, 20.3856192883815}; + std::vector PlanckOPL_CH4 = {0.6939169829570981, 0.2276100632203682, + 0.18937636253724896, 0.19815582280865973, 0.2200025240118599, 0.2534843471195877, + 0.3046841472491959, 0.3823986600184083, 0.4987623999795097, 0.6713411345124516, + 0.9259289806598495, 1.3002850424277999, 1.849189501228091, 2.6564473100322914, + 3.8304102615717746, 5.5344722927178935, 7.994950093215335, 11.526582452228748, + 16.563118973976213, 23.69700078531631, 33.73017627335858, 47.73920940736953, + 67.1570036554854, 93.87559477159049, 130.37319123976783, 179.87209418262805, + 246.53023232205777, 335.6763859293513, 454.0906008652992, 610.3469544483086, + 815.2123010446278, 1082.1303102619831, 1427.7830218942922, 1872.756285795192}; + std::vector PlanckOPL_OH = {0.009691578817706287, 0.030049601179558346, + 0.06765464326386915, 0.12708598866615728, 0.21115912062446549, 0.3200686268532878, + 0.4522763199038495, 0.6063045700105023, 0.7818952037702664, 0.980162180031524, + 1.2032232614189733, 1.453779407069445, 1.7348075847988176, 2.0493997120339555, + 2.4006745319830265, 2.7917401665539523, 3.2256754466853086, 3.7054691023958872, + 4.23395822768527, 4.813594141021041, 5.446035279732838, 6.131334910130047, + 6.866656642167146, 7.644231944896166, 8.44868549777347, 9.253875110774946, + 10.020328484366413, 10.694938827176772, 11.215426325696, 11.520789003674768, + 11.56668421121504, 11.33994851879643, 10.864608936023622, 10.19585170664347}; + std::vector PlanckOPL_N2O = {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}; + std::vector PlanckOPL_NH3 = {0.015853113133010365, 0.01926676385432445, + 0.026415824676239527, 0.03897633796851235, 0.05888157907336428, 0.08915105723552613, + 0.13355327539124195, 0.1960730549285638, 0.28048642863977286, 0.3904687215819763, + 0.5303141578042857, 0.7058788022739142, 0.925329453596419, 1.199577110754361, + 1.5425833020806663, 1.9716075385766434, 2.507579090345429, 3.175515361108895, + 4.005017660292852, 5.030732259316466, 6.292943371621803, 7.838191292724922, + 9.719463627047077, 11.997068455181829, 14.738974683444567, 18.02055752855426, + 21.925832857763382, 26.548007401632894, 31.987628299146124, 38.35527359161105, + 45.770040671919816, 54.35989812846817, 64.26225033961265, 75.62383376369114}; + std::vector PlanckOPL_H2O2 = {0.0061240643487373726, 0.014358654955073104, + 0.02372207214832294, 0.03509853861756539, 0.05041287435495022, 0.07158406668799877, + 0.10079925820290134, 0.14072600135394597, 0.1946390235350878, 0.26650037055313247, + 0.3610028583328736, 0.48358138449031224, 0.6403937124726201, 0.8382784165969821, + 1.0846948121183866, 1.3876522056917315, 1.7556322223573986, 2.1975119501537628, + 2.722491030278221, 3.340011592913925, 4.059702755928071, 4.891311219712732, + 5.844656956185561, 6.929574259055882, 8.155886163115925, 9.533355173169655, + 11.071661643874863, 12.780377074549552, 14.668934187405775, 16.746625569074368, + 19.022581230831573, 21.50576621761172, 24.204937981243575, 27.128677559801766}; + std::vector PlanckOPL_H2CO = {1.3851378337328886, 0.29966707757879113, + 0.14782184113680305, 0.11321165888878547, 0.10667714260745761, 0.11262176302472546, + 0.12740720489674168, 0.15069851692677394, 0.18357419050287804, 0.22805478242524305, + 0.2870105594836183, 0.3642229104940912, 0.464518116525414, 0.5939600794802649, + 0.7600869230903206, 0.9722007453399172, 1.2417019725072342, 1.5824881073722192, + 2.011408496739065, 2.5487901402831663, 3.2190405187317457, 4.051345642825396, + 5.080436220225772, 6.347502018651657, 7.901169358594316, 9.798657874500417, + 12.107045287048898, 14.904679352835634, 18.28279597617536, 22.347272505969286, + 27.220616304368882, 33.044168742639215, 39.9804690379451, 48.216066577742545}; + std::vector PlanckOPL_NO2 = {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}; + std::vector PlanckOPL_NO = {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}; + std::vector PlanckOPL_HCN = {0.028635211450903347, 0.03757526481800404, + 0.061354128074583666, 0.099041934820521, 0.14986680872525904, 0.2112860122528754, + 0.2821606517643699, 0.3648298945235665, 0.4645287980339424, 0.5924143164459074, + 0.7542345476291972, 0.9604408381860212, 1.2237054355039148, 1.5594639696461987, + 1.986552912551171, 2.527919503080509, 3.211397414785762, 4.070563895597286, + 5.145712796528703, 6.484886092047459, 8.144985631727018, 10.193065855766827, + 12.707593159716264, 15.779896881489458, 19.51568964842916, 24.03665133695894, + 29.482042907652733, 36.01052705210453, 43.802207627836076, 53.06005801030469, + 64.01228918330375, 76.91425579492252, 92.05099992084168, 109.73826951057255}; + std::vector PlanckOPL_C2H2 = {0.010929525364339873, 0.01392974120512336, + 0.02405443205882381, 0.04386508025067529, 0.07855266240755755, 0.1345672180286686, + 0.21996737264216254, 0.3461990874422214, 0.5305814978967278, 0.798986833061833, + 1.18893005089752, 1.7534946898618935, 2.5663306054757027, 3.727852933362788, + 5.372529800880729, 7.677251759511393, 10.870517832651087, 15.242285828331596, + 21.154028397137477, 29.048837769380295, 39.46135154206141, 53.02686745889668, + 70.48958296253277, 92.71027377239228, 120.6719646353374, 155.48532979799032, + 198.39202277505726, 250.76660329085342, 314.1168538786097, 390.087262270045, + 480.45331608083166, 587.1190933644476, 712.117974644975, 857.6066079444056}; + std::vector PlanckOPL_C2H4 = {0.01654680957667988, 0.02492716585471849, + 0.0460983767688889, 0.09024057897411016, 0.17647152655342777, 0.33699482855742513, + 0.6253345063344479, 1.130036554810777, 1.9963077428931275, 3.4594563059971244, + 5.895932142906171, 9.90032557003133, 16.40006735361465, 26.82443744854234, + 43.35046571234703, 69.25728060994958, 109.43037257639607, 171.07357189711118, + 264.7028534520616, 405.52279353849264, 615.3126797518005, 924.9981677988787, + 1378.1207657380173, 2035.4777754072459, 2981.2965059567937, 4331.43337255344, + 6244.0321979790415, 8933.516689165475, 12688.676814602117, 17895.963748387712, + 25069.504472771623, 34888.98940246357, 48247.80480294268, 66314.77809464224}; + std::vector PlanckOPL_C2H6 = {0.38479552555827473, 0.330335973067736, + 0.4412642688261878, 0.7234904395389986, 1.3047333754840624, 2.4582817867678797, + 4.728954047788753, 9.185141033669707, 17.899543021230603, 34.850795560924084, + 67.58893234939407, 130.26844542962587, 249.09122926085138, 471.9318246172849, + 885.130031885148, 1642.3336783887491, 3013.480869817436, 5466.6396914389925, + 9803.371388627636, 17379.457798083575, 30460.7202083096, 52790.33601327037, + 90480.336797945, 153407.43112853833, 257355.273339689, 427303.61580306507, + 702366.9359563864, 1143267.5642510268, 1843354.9491235327, 2944846.812649131, + 4662659.170473707, 7318745.5983123155, 11391692.329199707, 17587602.579628713}; + std::vector PlanckOPL_C2N2 = {0.005068498559075421, 0.00794066986454325, + 0.011326690006467697, 0.0162428641418049, 0.023999391745476977, 0.03635127840260863, + 0.05584914808886189, 0.08625527503938223, 0.1330651839347766, 0.20418883262957285, + 0.3108436682130422, 0.4687112261455473, 0.699423764955656, 1.0324479058510616, + 1.507453218412811, 2.1772618202142606, 3.1114775402859705, 4.400946829376023, + 6.163150049802422, 8.54873395424021, 11.749306098379188, 16.006768024926618, + 21.624282462568466, 28.979328225842597, 38.538809411076464, 50.8768807730018, + 66.6955409114437, 86.84833705611524, 112.36801734231409, 144.49770280025078, + 184.7272696900555, 234.83366678173417, 296.92842179122425, 373.5097598037753}; + std::vector PlanckOPL_C4H2 = {0.008654385959387364, 0.013221048546382841, + 0.023444997915175448, 0.04492768707054329, 0.09001749431190405, 0.1838640708259934, + 0.3765391487893722, 0.7654115015689522, 1.5353750432524051, 3.0294321701285036, + 5.869712810223213, 11.160509846374135, 20.82201029603419, 38.1284600107405, + 68.56051866929708, 121.13277007136199, 210.4288561712099, 359.67571007269487, + 605.3226740532338, 1003.7744131679052, 1641.1702181707305, 2647.42124377185, + 4216.140762151766, 6632.661018999894, 10313.01303955769, 15857.676154221808, + 24125.02967067159, 36330.98667841451, 54182.78820638309, 80057.91631401845, + 117240.68658474933, 170234.22688919335, 245167.68464599433, 350326.5016591533}; + std::vector PlanckOPL_HCOOH = {0.008500553910125858, 0.009363497122291741, + 0.01283858517065757, 0.01895258140461621, 0.02911801074174551, 0.045793592077955574, + 0.07290068695530225, 0.11656063070671369, 0.18618867990418764, 0.2960655295828281, + 0.46756168130170694, 0.7322445505013475, 1.1361637959006956, 1.745696069866044, + 2.655423918456038, 3.9986624722081285, 5.9613689634328955, 8.800366773348818, + 12.86702911179098, 18.637779834813458, 26.753121652310554, 38.06717802336468, + 53.710058801743436, 75.16625089593722, 104.37193238145707, 143.8355434073111, + 196.78635835589836, 267.3562198557535, 360.8011265516046, 483.76938983787875, + 644.6270384585974, 853.8465812437338, 1124.47299892604, 1472.682051768596}; + std::vector PlanckOPL_HNO3 = {0.006384561215332105, 0.007945818865473717, + 0.011369606874221595, 0.018185332096860546, 0.030917439657766183, 0.05402995712473861, + 0.0952035718822803, 0.16731021826266373, 0.29144824024627225, 0.501507552422749, + 0.8509050408969643, 1.4223515155012039, 2.341789738404136, 3.797990929638042, + 6.069741193008762, 9.563079645542468, 14.861694346418044, 22.794395726651224, + 34.52449765491417, 51.66720541762567, 76.44215125950534, 111.87016843316668, + 162.0251932450192, 232.3537440607607, 330.07832708074693, 464.7023518143656, + 648.6382771969264, 897.9856543786034, 1233.4873962465554, 1681.7015375997942, + 2276.4250569677397, 3060.4203470895163, 4087.4985253227364, 5425.027223873235}; + std::vector PlanckOPL_O3 = {0.06090151715462474, 0.04094792168395479, + 0.0493930255271001, 0.06968082884526133, 0.10301858705107333, 0.1546948114590627, + 0.23334115457989626, 0.35145488255968, 0.5262795613218485, 0.7808920244469918, + 1.1454338063269254, 1.658468204302404, 2.368506909036788, 3.335619739243394, + 4.633232474776676, 6.350081219500141, 8.592252521288646, 11.485497395583318, + 15.17760366231008, 19.84106558507588, 25.67573909058048, 32.91200786023972, + 41.81367130298549, 52.68162087288378, 65.85701739068608, 81.72543958358092, + 100.7206904850768, 123.3288052772987, 150.09273327833398, 181.61685811709106, + 218.5722288990807, 261.7005632383693, 311.82119282475406, 369.8355675948166}; + + std::map> PlanckOPLMap; + + PlanckOPLMap.insert({"CO2", PlanckOPL_CO2}); + PlanckOPLMap.insert({"H2O", PlanckOPL_H2O}); + PlanckOPLMap.insert({"CO", PlanckOPL_CO}); + PlanckOPLMap.insert({"CH4", PlanckOPL_CH4}); + PlanckOPLMap.insert({"OH", PlanckOPL_OH}); + PlanckOPLMap.insert({"N2O", PlanckOPL_N2O}); + PlanckOPLMap.insert({"NH3", PlanckOPL_NH3}); + PlanckOPLMap.insert({"H2O2", PlanckOPL_H2O2}); + PlanckOPLMap.insert({"H2CO", PlanckOPL_H2CO}); + PlanckOPLMap.insert({"NO2", PlanckOPL_NO2}); + PlanckOPLMap.insert({"NO", PlanckOPL_NO}); + PlanckOPLMap.insert({"HCN", PlanckOPL_HCN}); + + PlanckOPLMap.insert({"C2H2", PlanckOPL_C2H2}); + PlanckOPLMap.insert({"C2H4", PlanckOPL_C2H4}); + PlanckOPLMap.insert({"C2H6", PlanckOPL_C2H6}); + PlanckOPLMap.insert({"C2N2", PlanckOPL_C2N2}); + PlanckOPLMap.insert({"C4H2", PlanckOPL_C4H2}); + PlanckOPLMap.insert({"HCOOH", PlanckOPL_HCOOH}); + PlanckOPLMap.insert({"HNO3", PlanckOPL_HNO3}); + PlanckOPLMap.insert({"O3", PlanckOPL_O3}); + + // natural logarithms of reversed optical path lengths + // for linear interpolation in logarithm scale // 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; + + double coef = 0.0; + // loop over all grid points for (size_t j = jmin; j < jmax; j++) { // helping variable for the calculation double radiative_heat_loss = 0; - - // 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 (int k = 0; k < OPL_table_size; k++) { - if (T(x, j) < TemperatureOPL[k]) { - 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; + // temperature table interval search + int T_index = 0; + + for (int k = 0; k < OPL_table_size; k++) { + if (T(x, j) < TemperatureOPL[k]) { + if (T(x, j) < TemperatureOPL[0]) { + T_index = 0;//lower table limit } else { - coef=OPL_log_H2O[OPL_table_size-1]; + T_index = k; } + break; } - 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 (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]; - } + else { + T_index=OPL_table_size-1;//upper table limit } - 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; + // calculation of the mean Planck absorption coefficient + double k_P = 0.0; + + for(std::string sp_name : AbsorptionSpeciesList) { + // absorption coefficient for specie + if (AbsorptionSpeciesMap[sp_name] != npos) { + double k_P_specie = 0.0; + if ((T_index == 0) || (T_index == OPL_table_size-1)) { + coef=log(1.0/PlanckOPLMap[sp_name][T_index]); } else { - coef=OPL_log_OH[OPL_table_size-1]; + coef=log(1.0/PlanckOPLMap[sp_name][T_index-1])+ + (log(1.0/PlanckOPLMap[sp_name][T_index])-log(1.0/PlanckOPLMap[sp_name][T_index-1]))* + (T(x, j)-TemperatureOPL[T_index-1])/(TemperatureOPL[T_index]-TemperatureOPL[T_index-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; + k_P_specie = exp(coef); + + k_P_specie /= k_P_ref; + k_P += m_press * X(x, AbsorptionSpeciesMap[sp_name], j) * k_P_specie; + } } + // 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); From ea5d074f5deb5fcb24b3d89d970fca085da5842d Mon Sep 17 00:00:00 2001 From: lavrenyukiv Date: Tue, 2 Feb 2021 03:18:59 +0500 Subject: [PATCH 6/7] Update StFlow.h - Definitions for radiation model added --- include/cantera/oneD/StFlow.h | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/include/cantera/oneD/StFlow.h b/include/cantera/oneD/StFlow.h index bfae2394f3..c960485529 100644 --- a/include/cantera/oneD/StFlow.h +++ b/include/cantera/oneD/StFlow.h @@ -450,7 +450,9 @@ class StFlow : public Domain1D //! Indices within the ThermoPhase of the radiating species. First index is //! for CO2, second is for H2O. - std::vector m_kRadiating; + std::map AbsorptionSpeciesMap; + std::vector AbsorptionSpeciesList; + // flags std::vector m_do_energy; From 89cc361369a5710abe8b83032073a4d0e58dbd1e Mon Sep 17 00:00:00 2001 From: lavrenyukiv Date: Tue, 2 Feb 2021 04:05:38 +0500 Subject: [PATCH 7/7] Update StFlow.cpp --- src/oneD/StFlow.cpp | 17 ++++++++++------- 1 file changed, 10 insertions(+), 7 deletions(-) diff --git a/src/oneD/StFlow.cpp b/src/oneD/StFlow.cpp index a263b70447..9ee4742af1 100644 --- a/src/oneD/StFlow.cpp +++ b/src/oneD/StFlow.cpp @@ -308,13 +308,16 @@ void StFlow::evalResidual(double* x, double* rsd, int* diag, // 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. - // The model considers only CO2 and H2O as radiating species. Polynomial - // lines calculate the species Planck coefficients for H2O and CO2. The data - // for the lines is taken from the RADCAL program [Grosshandler, W. L., - // RADCAL: A Narrow-Band Model for Radiation Calculations in a Combustion - // Environment, NIST technical note 1402, 1993]. The coefficients for the - // polynomials are taken from [http://www.sandia.gov/TNF/radiation.html]. - + // + // 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 evaluated at temperature range 200 - 3500 K + // and are tabulated. + if (m_do_radiation) { // variable definitions for the Planck absorption coefficient and the // radiation calculation: