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 3 commits into
base: main
Choose a base branch
from
Open

Revive 1D radiation improvements of PR965 #1799

wants to merge 3 commits into from
+613 −78

Conversation

wandadars
Copy link
Contributor

@wandadars wandadars commented Oct 22, 2024 •
edited
Loading

I was looking through the open pull requests and saw #965 . It seemed like it was something that could be pushed across the finish line given the recent updates to Cantera that allow for easier handling of data in the YAML file.

Potentially allows for the discussions in Cantera/enhancements#72 to be realized.

I added some of the data from the original pull request into the gri30.yaml file (attached), and ran the 1D diffusion flame example

gri30.txt

From the discussions in the #965 it seemed like there was more than one possibility for handling that absorption coefficient calculation, so I added a keyword of fit-type which can be table or polynomial. The polynomial type is what was originally used and the table is the log linear interpolation data that was apparently better. PMAC is Plank Mean Absorption Coefficient. I'm pretty clueless on radiation modeling, so if anyone else knows more about this, please chime in.

As an example this would be the YAML entry for a species C2H6 that used the log-linear interpolation tabulated data.

- name: C2H6
  composition: {C: 2, H: 6}
  thermo:
    model: NASA7
    temperature-ranges: [200.0, 1000.0, 3500.0]
    data:
    - [4.29142492, -5.5015427e-03, 5.99438288e-05, -7.08466285e-08, 2.68685771e-11,
      -1.15222055e+04, 2.66682316]
    - [1.0718815, 0.0216852677, -1.00256067e-05, 2.21412001e-09, -1.9000289e-13,
      -1.14263932e+04, 15.1156107]
    note: L8/88
  transport:
    model: gas
    geometry: nonlinear
    well-depth: 252.3
    diameter: 4.302
    rotational-relaxation: 1.5
  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.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]

And for the case where the older polynomial fits (I'm not sure how prevalent/standardized this form is) are used, we specify using:

  radiation:
    model: PMAC
    fit-type: polynomial
    data: [-0.23093, -1.12390, 9.41530, -2.99880, 0.51382, -1.86840e-5]

Tagging the original creator of the pull request and the user who opened the enhancement. @lavrenyukiv @BYUignite

@BYUignite
Copy link

BYUignite commented Oct 30, 2024 via email

@ischoegl ischoegl mentioned this pull request Nov 5, 2024
Closed
@ischoegl
Copy link
Member

Now that Cantera 3.1 is released, this should be ready for further discussion and review.

Linking Cantera/enhancements#218 and tagging @marina8888 as there was interest.

@ischoegl ischoegl mentioned this pull request Dec 28, 2024
Closed
@marina8888
Copy link

marina8888 commented Jan 3, 2025 •
edited
Loading

Hi All,

I've been using the polynomials with high and low temperature ranges in my work by putting them into computeRadiation (excuse the ugly code below, but actually it seems quite similar to #965). Integrating RadLib into Cantera is way above my skill grade, but if no one else is working on this, perhaps, for now, I could rewrite this function for parsing yaml data to Flow1D computeRadiation - in a similar format to the SpeciesThermoFactory.cpp file? Assuming StFlow will be depreciated...

Some coefficients taken from the literature result in solver errors, so I think it will need some tests on the input polynomials, but I don't have a good grasp of the issue at the moment.

Flow1D.cpp:

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;

    // Polynomial coefficients:
 const double c_H2O[6] = {-0.23093, -1.12390, 9.41530, -2.99880, 0.51382, -1.86840e-5};
 const double low_temp_NO[7] = {3.09778e+01, −3.76125e-01, 1.78622e-03, −4.19616e-06, 5.31065e-09, −3.49672e-12, 9.44469e-16};
 const double high_temp_NO[7] = {3.32933e+00, 4.96445e-03, −1.94120e-05, 2.00931e-08, −9.88762e-12, 2.40521e-15, −2.33153e-19};
  const double low_temp_N2O[7] = {9.47313e+01, -1.05713e+00, 4.90429e-03, -1.04877e-05, 1.16578e-08, -6.67397e-12, 1.57665e-15};
  const double high_temp_N2O[7] = {2.72619e+01, 1.26580e-01, −3.75408e-04, 3.78748e-07, -1.86324e-10, 4.55295e-14, -4.43653e-18};
  const double low_temp_NH3[7] = {1.15691e+01, -3.56117e+04, 3.33438e+07, -7.08740e+09, 2.26363e+11, 6.26830e+13};
  const double high_temp_NH3[7] = {−5.10770e-02, −8.61034e+02, 9.23915e+06, −2.94040e+10, 3.65220e+13, −1.29468e+16};
  const double c_H2O[6]  = {0,0,0,0,0,0};
  const double c_H2O[6] = {-0.23093, -1.12390, 9.41530, -2.99880, 0.51382, -0.0000186840};
  const double low_temp_NO[7] = {30.9778, -0.376125, 0.00178622, -0.00000419616, 0.00000000531065, -0.00000000000349672, 0.000000000000944469};
    const double high_temp_NO[7] = {3.32933, 0.00496445, -0.0000194120, 0.0000000200931, -0.00000000000988762, 0.00000000000240521, -0.0000000000000000233153};
    const double low_temp_N2O[7] = {94.7313, -1.05713, 0.00490429, -0.0000104877, 0.0000000116578, -0.00000000000667397, 0.00000000000157665};
    const double high_temp_N2O[7] = {27.2619, 0.126580, -0.000375408, 0.000000378748, -0.000000000186324, 0.0000000000000455295, -0.00000000000000000443653};

    // 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);

    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[0] != 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[0], j) * k_P_H2O;
        }
       // Absorption coefficient for NO (using temperature-dependent ranges)
        if (m_kRadiating[1] != npos) {
            double k_P_NO = 0;
            if (T(x, j) <= 785) {
                // Use low temperature polynomial for NO
                for (size_t n = 0; n <= 6; n++) {
                    k_P_NO += low_temp_NO[n] * pow(1000 / T(x, j), (double)n);
                }
            } else {
                // Use high temperature polynomial for NO
                for (size_t n = 0; n <= 6; n++) {
                    k_P_NO += high_temp_NO[n] * pow(1000 / T(x, j), (double)n);
                }
            }
            k_P_NO /= k_P_ref;
            k_P += m_press * X(x, m_kRadiating[1], j) * k_P_NO;
        }

        // Absorption coefficient for N2O (using temperature-dependent ranges)
        if (m_kRadiating[2] != npos) {
            double k_P_N2O = 0;
            if (T(x, j) <= 820) {
                // Use low temperature polynomial for N2O
                for (size_t n = 0; n <= 6; n++) {
                    k_P_N2O += low_temp_N2O[n] * pow(1000 / T(x, j), (double)n);
                }
            } else {
                // Use high temperature polynomial for N2O
                for (size_t n = 0; n <= 6; n++) {
                    k_P_N2O += high_temp_N2O[n] * pow(1000 / T(x, j), (double)n);
                }
            }
            k_P_N2O /= k_P_ref;
            k_P += m_press * X(x, m_kRadiating[2], j) * k_P_N2O;
        }

        // Absorption coefficient for NH3 (using temperature-dependent ranges)
        if (m_kRadiating[3] != npos) {
            double k_P_NH3 = 0;
            if (T(x, j) <= 1100) {
                // Use low temperature polynomial for NH3
                for (size_t n = 0; n <= 6; n++) {
                    k_P_NH3 += low_temp_NH3[n] * pow(1000 / T(x, j), (double)n);
                }
            } else {
                // Use high temperature polynomial for NH3
                for (size_t n = 0; n <= 6; n++) {
                    k_P_NH3 += high_temp_NH3[n] * pow(1000 / T(x, j), (double)n);
                }
            }
            k_P_NH3 /= k_P_ref;
            k_P += m_press * X(x, m_kRadiating[3], j) * k_P_NH3;
        }

        // 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);

@wandadars
Copy link
Contributor Author

I've been a bit busy the last few months. I recall from my investigation that there was generally multiple ways of obtaining the values that are needed to compute the radiative heat loss. At a high level, the interface between say the Cantera solver and the radiation coding would be something like the species, their mass/mole fraction in a mixture, and the temperature. I think getting a handle on the interface is the highest priority. From the Radlib website it looks like there's 3 major models that are usable.

1.) The Planck Mean (PM) model.
2.) The Weighted Sum of Gray Gases (WSGG) model,
3.) The Rank Correlated Spectral Line Weighted Sum of Gray Gases (RCSLW) model

This is where the folks who work with this library can help to clarify some things. So there's a single global equation that all radiation modeling boils down to trying to get the best estimate of the parameters for? In the current Cantera coding, the final thing that is needed by Cantera is this k_P parameter, which is then used to compute the volumetric heat loss:

double radiative_heat_loss = 2 * k_P *(2 * StefanBoltz * pow(T(x, j), 4) - boundary_Rad_left - boundary_Rad_right);

If that's the case, then somewhere we would have a control for which of the 3 radiation models to use, which then would tell Cantera to setup whatever is needed by the particular model and potentially instantiate a radiation object. From this example it looks like the interface needs to be something like:

//----------------- create radiation objects
 
    rad *planckmean = new rad_planck_mean();
    rad *wsgg       = new rad_wsgg();
    rad *rcslw      = new rad_rcslw(4, T, P, fvsoot, xH2O, xCO2, xCO);
 
    //----------------- compute absorption coefficients and weights
 
    planckmean->get_k_a(kabs_pm, awts_pm,       T, P, fvsoot, xH2O, xCO2, xCO, xCH4);
    wsgg->get_k_a(      kabs_wsgg, awts_wsgg,   T, P, fvsoot, xH2O, xCO2);
    rcslw->get_k_a(     kabs_rcslw, awts_rcslw, T, P, fvsoot, xH2O, xCO2, xCO);

How this k_a quantity relates to the needed k_P is not yet clear to me.

I think something like a factory during initalization that picks the type of radiation model that a user wants to use and creates the object and stores it in something like a generic radiation_object variable name that is then called during execution by the computeRadiation() method in Cantera, with the varying argument lists being handled by an if-statement that checks the model type and provides whatever values are needed to the get_k_a() method and then proceeds through an operation to use the results from this call to generate the k_P value.

@BYUignite
Copy link

In the example noted, k_P likely refers to a Planck mean absorption coefficient that is applied in an optically-thin model to compute the radiative_heat_loss. The equivalent in RadLib would be the
rad *planckmean = new rad_planck_mean();
model selection. More complex and accurate models are the family of weighted sum of grey gases in which the radiation spectrum is written in terms of multiple "gases", which can be thought of as discrete bands of radiation (loosely, and to be confused, unfortunately, with actual gases). So, Planck Mean assumes a single "gas". There are then two parts to a radiation problem (1) the radiation properties for the gray gases and (2) the radiation solver. RadLib provides radiation properties for the gases. Those would then be coupled to a user code for the radiation solver (like optically thin, or two-flux, or discrete ordinates, etc.). There are three models currently: Planck Mean, WSGG (Bordbar), and RCSLW. The first is the simplest, with just one gas. The WSGG hard-codes the number of gases (built-in to the formulation/correlation; and the RCSLW allows a variable number. All three models take input as a local gas state, and output the absorption coefficients (the k's) for each "gas" and the weighting factors (the a's). For the Planck Mean model, the single a value is just a=1 and is not normally used or considered. But both the k's and the a's are used for the WSGG and the RCSLW and show up in the radiative transport equations for spectral (gray-gas style) formulations. I'm happy to discuss some more, including video, etc. as desired. -- David Lignell

@speth
Copy link
Member

speth commented Jan 4, 2025

Glad to see the continued discussion of this. I'd suggest it would be better to share these comments in an "Issue" (Cantera/enhancements#72) rather than as PR comments, as the latter are unfortunately likely to get automatically hidden by GitHub as more commits and code comments are added.

ischoegl
ischoegl reviewed Jan 5, 2025
View reviewed changes
Copy link
Member

@ischoegl ischoegl left a comment •
edited
Loading

Choose a reason for hiding this comment

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

Thank you everyone for the input. While I think that Cantera/enhancements#72 may be a little too specific and potentially beyond the scope of the conversation here, I believe the path forward is relatively clear:

  • Long-term, we would like to use a newRadiationModel factory method that supports RadLib variants as envisioned by @BYUignite.
  • Short-term, we'd like to include additional parameterizations along the lines of @lavrenyukiv's original PR Radiation model update in StFlow.cpp #965 (which appears to have parallels with @marina8888's intent) while moving away from hard-coded radiation parameterizations to YAML as was added in this PR by @wandadars.

I'd suggest moving forward with the original approach of this PR and working out the best possible compromise that allows for future enhancements. Ideally, this would involve a new Radiation1D class that can be instantiated based on YAML directly. It should provide an interface for computeRadiation that replaces the current Flow1D::computeRadiation method.

For the YAML parameterization itself, I am a little hesitant to add this as a new radiation field to species specifications, as it is quite specific to a use case, hard to maintain, and will become unwieldy for large mechanisms. Instead, I'd suggest adding a stand-alone parameterization file radiation_parameters.yaml, perhaps with sections for different models, i.e.:

# radiation_parameters.yaml
PMAC:
- H2O:
  fit-type: polynomial
  data: [-0.23093, -1.12390, 9.41530, -2.99880, 0.51382, -1.86840e-5]
- 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]
[...]

Similar to other entries, this could be linked in the main YAML file, e.g.

# air_radiation.yaml
phases:
- name: gas
  thermo: ideal-gas
  elements: [O, H, N, Ar]
  species:
  - gri30.yaml/species: all
  skip-undeclared-elements: true
  skip-undeclared-third-bodies: true
  kinetics: bulk
  reactions:
  - gri30.yaml/reactions: declared-species  
  transport: mixture-averaged
  radiation:  # <--------------------------------- new entry
  - radiation_parameters.yaml

or, alternatively, passed as a parameter to the RadiationModel constructor? (mismatched species names could be addressed with the Phase::addSpeciesAlias route).

If this is ok with everyone, I'd suggest rebasing the PR and proceeding with a standard review.

@wandadars wandadars force-pushed the 1d_radiation_update branch from 336a341 to 7cccced Compare January 10, 2025 02:52
@wandadars
Copy link
Contributor Author

wandadars commented Jan 10, 2025 •
edited
Loading

Thanks for all the comments everyone. I think @ischoegl 's suggestion is a good solution towards enhancing the way Cantera handles radiative effects.

I'm not very well versed in how the YAML parsing works beyond how I see the usage that happens in the Peng Robinson EOS class and the way we have it in the radiation coding in the PR. I can't quite just call species->input if I have a structure of the radiation_properties.yaml file like the one you have shown where the PMAC is sort of the entry-point for the YAML.

@ischoegl
Copy link
Member

I'm not very well versed in how the YAML parsing works beyond how I see the usage that happens in the Peng Robinson EOS class and the way we have it in the radiation coding in the PR. I can't quite just call species->input if I have a structure of the radiation_properties.yaml file like the one you have shown where the PMAC is sort of the entry-point for the YAML.

Re YAML parsing: have a look at some of the main loaders, e.g.

cantera/src/base/Solution.cpp

Lines 203 to 207 in 2e0d8ac

shared_ptr<Solution> newSolution(const string& infile, const string& name,
const string& transport, const vector<string>& adjacent)
{
auto rootNode = AnyMap::fromYamlFile(infile);
const AnyMap& phaseNode = rootNode.at("phases").getMapWhere("name", name);

and associated calls, e.g.

cantera/src/thermo/ThermoFactory.cpp

Line 183 in 2e0d8ac

void addSpecies(ThermoPhase& thermo, const AnyValue& names, const AnyValue& species)

Reading YAML input is taken care of by AnyMap::fromYamlFile, and the rest is relatively straightforward. Ray has included very good documentation for AnyMap ... see C++ documentation

@wandadars wandadars force-pushed the 1d_radiation_update branch from 67555c5 to f4aabb8 Compare January 13, 2025 19:05
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Reviewers

@ischoegl ischoegl ischoegl left review comments

At least 1 approving review is required to merge this pull request.

Assignees
No one assigned
Labels
1D Flames enhancement
Projects
None yet
Milestone
No milestone
Development

Successfully merging this pull request may close these issues.
5 participants
@wandadars @BYUignite @ischoegl @marina8888 @speth