Commit 8d69ba1d authored by Roberto Porcu's avatar Roberto Porcu Committed by Jordan Musser
Browse files

Improve reactions section (!76)

parent 612a79b2
Loading
Loading
Loading
Loading
+1 −0
Original line number Diff line number Diff line
@@ -58,6 +58,7 @@ To learn more about the implementation, follow the following reference sections:
   references/ManagingGridHierarchy_Chapter
   references/Fluids_Chapter
   references/Particles_Chapter
   references/Chemistry_Chapter
   references/EB_Chapter
   references/hpc
   references/mpmd
+239 −0
Original line number Diff line number Diff line
.. _Chap:Chemistry:

Chemical Reactions
==================

MFIX-Exa categorizes chemical reactions into two classes:

* **Eulerian reactions** are homogeneous fluid reactions, reactions that only
  reference fluid species. These reactions are calculated by looping over grid cells,
  and are defined in units of moles per second per cubic meter
  :math:`[mol \cdot s^{-1} \cdot m^{-3}]`.

* **Lagrangian reactions** include homogeneous particle reactions and
  heterogeneous fluid-particle reactions. These reactions are calculated
  by looping over each particle, and are defined in units of moles per second
  :math:`[mol \cdot s^{-1}]`.

Users implement reaction rates by modifying the header file,
``src/usr/mfix_usr_reactions_rates_K.H``.

The ``operator()`` method for the ``EulerianReactionRates`` class gets two
arguments, ``reactions`` and ``fluid``, which are two classes from which the user
can extract the data needed to implement the corresponding rates.

.. code-block:: cpp

   class EulerianReactionRates
   {
     public:
       /**
        *  Reactions class gives access to the following reactions quantities:
        *   - reaction_idx("reaction_name"), with reaction_name a string matching one of the reactions defined in inputs
        *   - rate(q), with q an integer from 0 to number of Eulerian reactions - 1
        *   - rate("reaction_name"), with reaction_name a string matching one of the reactions defined in inputs
        *   - time(), the current simulation time
        *   - dt(), the current simulation time step
        *   - cell_volume(), this cell geometric volume
        *   - vfrac(), this cell geometric volume fraction
        *
        *  Fluid class gives access to the following fluid quantities:
        *   - shear_viscosity()
        *   - volume_fraction()
        *   - density()
        *   - velocity() and velocity(dir), with dir an integer from 0 to 2
        *   - thermodynamic_pressure()
        *   - temperature()
        *   - specific_enthalpy()
        *   - species_idx("species_name"), with species_name a string matching one of the fluid species defined in inputs
        *   - mass_fraction(n), with n an integer from 0 to number of fluid species - 1
        *   - mass_fraction("species_name"), with species_name a string matching one of the fluid species defined in inputs
        *   - molar_concentration(n), with n an integer from 0 to number of fluid species - 1
        *   - molar_concentration("species_name"), with species_name a string matching one of the fluid species defined in inputs
        *   - molar_mass(n), with n an integer from 0 to number of fluid species - 1
        *   - molar_mass("species_name"), with species_name a string matching one of the fluid species defined in inputs
        *   - average_molar_mass()
        *   - mole_fraction(n), with n an integer from 0 to number of fluid species - 1
        *   - mole_fraction("species_name"), with species_name a string matching one of the fluid species defined in inputs
        */
       template <class Reactions, class Fluid>
       AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
       void operator() (Reactions& reactions,
                        const Fluid& fluid) const
       {
         // User defined Eulerian reactions rates go here
       }
   };

The ``operator()`` of the ``LagrangianReactionRates`` class contains a third argument,
``solids``, which is a class containing data for the particle involved in the reaction.

.. code-block:: cpp

   class LagrangianReactionRates
   {
     public:
       /**
        *  Reactions class gives access to the following reactions quantities:
        *   - reaction_idx("reaction_name"), with reaction_name a string matching one of the reactions defined in inputs
        *   - rate(q), with q an integer from 0 to max number of Lagrangian reactions - 1
        *   - rate("reaction_name"), with reaction_name a string matching one of the reactions defined in inputs
        *   - time(), the current simulation time
        *   - dt(), the current simulation time step
        *   - cell_volume(), cell geometric volume of cell containing this particle
        *   - vfrac(), geometric volume fraction of cell containing this particle
        *
        *  Fluid class gives access to the following fluid quantities:
        *   - shear_viscosity()
        *   - volume_fraction()
        *   - density()
        *   - velocity() and velocity(dir), with dir an integer from 0 to 2
        *   - thermodynamic_pressure()
        *   - temperature()
        *   - specific_enthalpy()
        *   - species_idx("species_name"), with species_name a string matching one of the fluid species defined in inputs
        *   - mass_fraction(n), with n an integer from 0 to max number of fluid species
        *   - mass_fraction("species_name"), with species_name a string matching one of the fluid species defined in inputs
        *   - molar_concentration(n), with n an integer from 0 to max number of fluid species - 1
        *   - molar_concentration("species_name"), with species_name a string matching one of the fluid species defined in inputs
        *   - molar_mass(n), with n an integer from 0 to max number of fluid species - 1
        *   - molar_mass("species_name"), with species_name a string matching one of the fluid species defined in inputs
        *   - average_molar_mass()
        *   - mole_fraction(n), with n an integer from 0 to max number of fluid species - 1
        *   - mole_fraction("species_name"), with species_name a string matching one of the fluid species defined in inputs
        *
        *  Solids class gives access to the following solids quantities:
        *   - id()
        *   - cpu()
        *   - position() and position(dir), with dir an integer from 0 to 2
        *   - radius()
        *   - volume()
        *   - density()
        *   - mass()
        *   - oneOverI(), inverse of the particle momentum of inertia
        *   - statistical_weight()
        *   - velocity() and velocity(dir), with dir an integer from 0 to 2
        *   - temperature()
        *   - specific_heat_capacity()
        *   - species_idx("species_name"), with species_name a string matching one of the solids species defined in inputs
        *   - mass_fraction(n), with n an integer from 0 to max number of solids species - 1
        *   - mass_fraction("species_name"), with species_name a string matching one of the solids species defined in inputs
        *   - molar_concentration(n), with n an integer from 0 to max number of solids species - 1
        *   - molar_concentration("species_name"), with species_name a string matching one of the solids species defined in inputs
        *   - molar_mass(n), with n an integer from 0 to max number of solids species - 1
        *   - molar_mass("species_name"), with species_name a string matching one of the solids species defined in inputs
        *   - average_molar_mass()
        *   - mole_fraction(n), with n an integer from 0 to max number of solids species - 1
        *   - mole_fraction("species_name"), with species_name a string matching one of the solids species defined in inputs
        */
       template <class Reactions, class Fluid, class Solids>
       AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
       void operator() (Reactions& reactions,
                        const Fluid& fluid,
                        const Solids& solids) const
       {
         // User defined Lagrangian reactions rates go here
       }
   };

To simply defining chemical reactions, fluid and solids species-related quantities
can be accessed by calling the variable function with the species name specified
in the inputs file. For example, if the fluid contains a species called ``CO``,
then the species mass fraction for ``CO`` can be accessed as follows:

.. code-block:: cpp

   // Mass fraction of CO species
   amrex::Real const X_CO = fluid.mass_fraction("CO");

This feature allows users to write readable reaction rates, especially when
the number of chemical species is large. However, at runtime this approach
requires the solver to search for the species index associated with the name.
The computational cost of the search can be non-trivial because it is performed
for each iteration of the chemistry solver every time step, for all grid cells
and all particles.

A python-script is provided that substitutes the string references with their
corresponding indices to remove the runtime search. The tool is run by passing
the ``inputs`` file and the reaction header, ``mfix_usr_reactions_rates_K.H``.
The tool creates a backup copy of the original reaction header, appending a time stamp
for when the tool was run, before applying the optimizations.

Example usage:

.. code-block:: bash

  python tools/Chemistry/optimize_usr_rates.py inputs mfix/src/usr/mfix_usr_reactions_rates_K.H


Example rates UDFs
==================

The following code shows a simple reaction rate for methane combustion.
The molar concentrations for oxygen (``O2``) and methane (``CH4``) are accessed
using the names provided in the inputs file. Similarly, the reaction rate
is stored using the name given in the inputs file.

.. code-block:: bash
   :caption: Snippet of inputs defining fluid species and adiabatic flame reaction

   # Fluid model settings
   # -----------------------------------------------------------------------
   fluid.solve = fluid

   fluid.species = CH4 O2 CO2 H2O N2

   fluid.thermodynamic_pressure = 101325.
   fluid.reference_temperature = 298.15

   fluid.species = Air Vapor

   fluid.viscosity.molecular = constant
   fluid.viscosity.molecular.constant = 1.86e-5

   fluid.specific_heat = mixture

   fluid.thermal_conductivity          = constant
   fluid.thermal_conductivity.constant = 0.02662


   # Chemical reactions:
   # -----------------------------------------------------------------------
   chemistry.solve = AdiabaticFlame

   chemistry.AdiabaticFlame.reaction = CH4(g) + 2O2(g) --> CO2(g) + 2H2O(g)


.. code-block:: cpp
   :caption: Snippet of reaction rate header defining reaction rate for adiabatic flame

   class EulerianReactionRates
   {
     public:
       AMREX_GPU_HOST_DEVICE
       EulerianReactionRates () = default;

       template <class Reactions, class Fluid>
       AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
       void operator() (Reactions& reactions,
                        const Fluid& fluid) const
       {

         // Adiabatic flame:
         // CH4 + 2O2 --> CO2 + 2H2O  (mol/m^3.s)
         //---------------------------------------------------------------------//
         // Note: The CH4 combustion rate is artificial and used for the
         // adiabatic flame test case.

         amrex::Real const c_O2  = fluid.molar_concentration("O2");   // [mol/m^3]
         amrex::Real const c_CH4 = fluid.molar_concentration("CH4");  // [mol/m^3]

         amrex::Real k0(2.0e-1); // reaction rate coefficient

         reactions.rate("AdiabaticFlame") = k0 * c_O2 * c_CH4;

       }
   };

.. toctree::
   :maxdepth: 0
+212 −21
Original line number Diff line number Diff line
Chemical Reactions
==================

Enabling the Chemical Reactions solver and specifying model options.

+----------------------------------+-------------------------------------------------------------+----------+-----------+
|                                  | Description                                                 |   Type   | Default   |
+==================================+=============================================================+==========+===========+
| chemistry.solve                  | Specified name(s) of the chemical reactions types or None   | String   |  None     |
|                                  | to disable the reactions solver.                            |          |           |
+----------------------------------+-------------------------------------------------------------+----------+-----------+

The following inputs are defined using the ``chemistry`` prefix.

+--------------------------+---------------------------------------------------------+----------+----------------+
|                          | Description                                             |   Type   | Default        |
+==========================+=========================================================+==========+================+
| solve                    | Specify the name(s) of the chemical reactions or set    |  String  |  None          |
|                          | to None to disable the chemistry solver. The name(s)    |          |                |
|                          | assigned to the chemistry solver are used to specify    |          |                |
|                          | the chemical reactions equations.                       |          |                |
+--------------------------+---------------------------------------------------------+----------+----------------+
| [reaction0].reaction     | Chemical formula for the given reaction. The string     |  String  |  None          |
|                          | given as input must not contain white spaces and        |          |                |
|                          | the reaction direction has to be specified as '-->'     |          |                |
|                          | or '<--'. Chemical species phases must be defined as    |          |                |
|                          | the reaction direction has to be specified as '\-\->'   |          |                |
|                          | or '<\-\-'. Chemical species phases must be defined as  |          |                |
|                          | '(g)' for the fluid phase or '(s)' for the solid phase. |          |                |
+--------------------------+---------------------------------------------------------+----------+----------------+
| mass_balance_tolerance   | Tolerance used to test chemical formula conserves mass. | Real     | 1e-12          |
@@ -81,8 +77,9 @@ is selected. The following inputs can be specified using the
+-----------------------+---------------------------------------------------------------------+--------+-----------+
| ode_max_steps         | The maximum number of substeps for the ODE integration              | int    | 150000    |
+-----------------------+---------------------------------------------------------------------+--------+-----------+
| jacobian_type         | Select how to compute the jacobian for the ODE iterative            | string | Numerical |
|                       | solver. Available choices are:                                      |        |           |
| jacobian_type         | Select how to compute the jacobian for the ODE iterative solver.    | string | Numerical |
|                       | This option is available for ``StiffSolver::BackwardEuler`` and     |        |           |
|                       | ``StiffSolver::VODE``. Available choices are:                       |        |           |
|                       |                                                                     |        |           |
|                       | * ``Numerical`` for a first-order numerical approximation           |        |           |
|                       | * ``Broyden`` for Broyden's method (only for BackwardEuler)         |        |           |
@@ -127,15 +124,209 @@ is selected. The following inputs can be specified using the

Below is an example for specifying chemical reactions for MFIX-Exa.

.. code-block:: none
Example inputs
==============

The following setup shows how to define chemical reactions in the inputs file.
Only ``Ideal Gas`` constraint types can be used when solving chemistry, as the
``IncompressibleFluid`` constraint would be inconsistent with the nature of
a problem where the fluid density is not constant because of mass transfer due
to chemical reactions.

.. code-block:: bash

  mfix.constraint_type  = IdealGasOpenSystem

  mfix.advect_density   = 1
  mfix.advect_enthalpy  = 1
  mfix.solve_species    = 1

  # Species model settings
  #----------------------------------------------------------------------
  species.solve  =  CH4  H2  CO  H2O  CO2  O2  N2  Fe2O3  FeO  inert

  species.diffusivity = constant
  species.diffusivity.constant  =  2.31000e-5  # m^2 / sec

  species.specific_heat = NASA7-poly

  # Methane
  species.CH4.molecular_weight  = 16.04276e-3
  species.CH4.specific_heat.NASA7.a0 =  5.14911468E+00    1.65326226E+00
  species.CH4.specific_heat.NASA7.a1 = -1.36622009E-02    1.00263099E-02
  species.CH4.specific_heat.NASA7.a2 =  4.91453921E-05   -3.31661238E-06
  species.CH4.specific_heat.NASA7.a3 = -4.84246767E-08    5.36483138E-10
  species.CH4.specific_heat.NASA7.a4 =  1.66603441E-11   -3.14696758E-14
  species.CH4.specific_heat.NASA7.a5 = -1.02465983E+04   -1.00095936E+04

  # Hydrogen
  species.H2.molecular_weight = 2.01588e-3
  species.H2.specific_heat.NASA7.a0 =  2.34433112E+00    2.93286575E+00
  species.H2.specific_heat.NASA7.a1 =  7.98052075E-03    8.26608026E-04
  species.H2.specific_heat.NASA7.a2 = -1.94781510E-05   -1.46402364E-07
  species.H2.specific_heat.NASA7.a3 =  2.01572094E-08    1.54100414E-11
  species.H2.specific_heat.NASA7.a4 = -7.37611761E-12   -6.88804800E-16
  species.H2.specific_heat.NASA7.a5 = -9.17935173E+02   -8.13065581E+02

  # Carbon monoxide
  species.CO.molecular_weight = 28.01040e-3
  species.CO.specific_heat.NASA7.a0 =  3.57953350E+00    3.04848590E+00
  species.CO.specific_heat.NASA7.a1 = -6.10353690E-04    1.35172810E-03
  species.CO.specific_heat.NASA7.a2 =  1.01681430E-06   -4.85794050E-07
  species.CO.specific_heat.NASA7.a3 =  9.07005860E-10    7.88536440E-11
  species.CO.specific_heat.NASA7.a4 = -9.04424490E-13   -4.69807460E-15
  species.CO.specific_heat.NASA7.a5 = -1.43440860E+04   -1.42661170E+04

  # Water vapor
  species.H2O.molecular_weight = 18.01528e-3
  species.H2O.specific_heat.NASA7.a0 =  0.41986352E+01    0.26770389E+01
  species.H2O.specific_heat.NASA7.a1 = -0.20364017E-02    0.29731816E-02
  species.H2O.specific_heat.NASA7.a2 =  0.65203416E-05   -0.77376889E-06
  species.H2O.specific_heat.NASA7.a3 = -0.54879269E-08    0.94433514E-10
  species.H2O.specific_heat.NASA7.a4 =  0.17719680E-11   -0.42689991E-14
  species.H2O.specific_heat.NASA7.a5 = -0.30293726E+05   -0.29885894E+05

  # Carbon dioxide
  species.CO2.molecular_weight = 44.00980e-3
  species.CO2.specific_heat.NASA7.a0 =  2.35681300E+00    4.63651110E+00
  species.CO2.specific_heat.NASA7.a1 =  8.98412990E-03    2.74145690E-03
  species.CO2.specific_heat.NASA7.a2 = -7.12206320E-06   -9.95897590E-07
  species.CO2.specific_heat.NASA7.a3 =  2.45730080E-09    1.60386660E-10
  species.CO2.specific_heat.NASA7.a4 = -1.42885480E-13   -9.16198570E-15
  species.CO2.specific_heat.NASA7.a5 = -4.83719710E+04   -4.90249040E+04

  # Oxygen
  species.O2.molecular_weight = 31.99880e-3
  species.O2.specific_heat.NASA7.a0 =  3.78245636E+00    3.66096065E+00
  species.O2.specific_heat.NASA7.a1 = -2.99673416E-03    6.56365811E-04
  species.O2.specific_heat.NASA7.a2 =  9.84730201E-06   -1.41149627E-07
  species.O2.specific_heat.NASA7.a3 = -9.68129509E-09    2.05797935E-11
  species.O2.specific_heat.NASA7.a4 =  3.24372837E-12   -1.29913436E-15
  species.O2.specific_heat.NASA7.a5 = -1.06394356E+03   -1.21597718E+03

  # Nitrogen
  species.N2.molecular_weight = 28.01340e-3
  species.N2.specific_heat.NASA7.a0 =  3.53100528E+00    2.95257637E+00
  species.N2.specific_heat.NASA7.a1 = -1.23660988E-04    1.39690040E-03
  species.N2.specific_heat.NASA7.a2 = -5.02999433E-07   -4.92631603E-07
  species.N2.specific_heat.NASA7.a3 =  2.43530612E-09    7.86010195E-11
  species.N2.specific_heat.NASA7.a4 = -1.40881235E-12   -4.60755204E-15
  species.N2.specific_heat.NASA7.a5 = -1.04697628E+03   -9.23948688E+02

  # Hematite
  species.Fe2O3.molecular_weight = 159.68820e-3
  species.Fe2O3.specific_heat.NASA7.a0 =  1.52218166E-01    2.09445369E+01
  species.Fe2O3.specific_heat.NASA7.a1 =  6.70757040E-02    0.00000000E+00
  species.Fe2O3.specific_heat.NASA7.a2 = -1.12860954E-04    0.00000000E+00
  species.Fe2O3.specific_heat.NASA7.a3 =  9.93356662E-08    0.00000000E+00
  species.Fe2O3.specific_heat.NASA7.a4 = -3.27580975E-11    0.00000000E+00
  species.Fe2O3.specific_heat.NASA7.a5 = -1.01344092E+05   -1.07936580E+05

  # Wustite
  species.FeO.molecular_weight = 71.84440e-3
  species.FeO.specific_heat.NASA7.a0 =  3.68765953E+00    1.81588527E+00
  species.FeO.specific_heat.NASA7.a1 =  1.09133433E-02    1.70742829E-02
  species.FeO.specific_heat.NASA7.a2 = -1.61179493E-05   -2.39919190E-05
  species.FeO.specific_heat.NASA7.a3 =  1.06449256E-08    1.53690046E-08
  species.FeO.specific_heat.NASA7.a4 = -2.39514915E-12   -3.53442390E-12
  species.FeO.specific_heat.NASA7.a5 = -3.34867527E+04   -3.30239565E+04

  # inert
  species.inert.molecular_weight = 92.82399e-3
  species.inert.specific_heat.NASA7.a0 =  1.52218166E-01    2.09445369E+01
  species.inert.specific_heat.NASA7.a1 =  6.70757040E-02    0.00000000E+00
  species.inert.specific_heat.NASA7.a2 = -1.12860954E-04    0.00000000E+00
  species.inert.specific_heat.NASA7.a3 =  9.93356662E-08    0.00000000E+00
  species.inert.specific_heat.NASA7.a4 = -3.27580975E-11    0.00000000E+00
  species.inert.specific_heat.NASA7.a5 = -1.01344092E+05   -1.07936580E+05


  # Fluid model settings
  #----------------------------------------------------------------------
  fluid.solve = fluid

  fluid.thermodynamic_pressure = 101325.
  fluid.reference_temperature = 298.15

  fluid.species =  CH4  H2  CO  H2O  CO2  O2  N2

  fluid.viscosity          = constant
  fluid.viscosity.constant = 4.73e-5  #air @ 1200K

  fluid.specific_heat = mixture

  fluid.thermal_conductivity          = constant
  fluid.thermal_conductivity.constant = 0.024


  # Solids model settings
  #----------------------------------------------------------------------
  solids.types = solid0
  solids.reference_temperature = 298.15

  solids.species  =  Fe2O3  FeO  inert

  solids.specific_heat = mixture


  # Chemistry model settings
  #----------------------------------------------------------------------
  chemistry.solve  =  Hem_CH4  Hem_H2  Hem_CO  Wus_O2

  chemistry.Hem_CH4.reaction = Fe2O3(s)+0.25CH4(g)-->2.FeO(s)+0.25CO2(g)+0.5H2O(g)
  chemistry.Hem_H2.reaction  = Fe2O3(s)+H2(g)-->2.FeO(s)+H2O(g)
  chemistry.Hem_CO.reaction  = Fe2O3(s)+CO(g)-->2.FeO(s)+CO2(g)
  chemistry.Wus_O2.reaction  = FeO(s)+0.25O2(g)-->0.5Fe2O3(s)


  # Initial Conditions
  #----------------------------------------------------------------------
  ic.regions = full-domain

  ic.full-domain.fluid.volfrac       = 0.725
  ic.full-domain.fluid.velocity      = 0.015  0.00  0.00
  ic.full-domain.fluid.temperature   = 1273.
  ic.full-domain.fluid.species.CH4   = 0.1
  ic.full-domain.fluid.species.H2    = 0.2
  ic.full-domain.fluid.species.CO    = 0.1
  ic.full-domain.fluid.species.H2O   = 0.2
  ic.full-domain.fluid.species.CO2   = 0.1
  ic.full-domain.fluid.species.O2    = 0.1
  ic.full-domain.fluid.species.N2    = 0.2

  ic.full-domain.solids  = solid0
  ic.full-domain.packing = pseudo_random

  ic.full-domain.solid0.volfrac  =  0.275
  ic.full-domain.solid0.velocity =  0.00  0.00  0.00

  ic.full-domain.solid0.diameter = constant
  ic.full-domain.solid0.density  = constant

  ic.full-domain.solid0.diameter.constant =  100.0e-6
  ic.full-domain.solid0.density.constant  = 2000.0

  ic.full-domain.solid0.temperature = 1273.
  ic.full-domain.solid0.species.Fe2O3 = 0.3
  ic.full-domain.solid0.species.FeO   = 0.3
  ic.full-domain.solid0.species.inert = 0.4

   chemistry.solve = my_reaction0 my_reaction1

   chemistry.my_reaction0.reaction = Fe2O3(s)+CO(g)-->2.FeO(s)+CO2(g)
   chemistry.my_reaction1.reaction = FeO(s)+0.25O2(g)-->0.5Fe2O3(s)
  # Boundary Conditions
  #----------------------------------------------------------------------
  bc.regions = inflow  outflow

   chemistry.mass_balance_tolerance = 1.e-5
  bc.inflow = mi
  bc.inflow.fluid.volfrac       = 1.0
  bc.inflow.fluid.velocity      = 0.015  0.0  0.0
  bc.inflow.fluid.temperature   = 1273.
  bc.inflow.fluid.species.CH4   = 0.1
  bc.inflow.fluid.species.H2    = 0.2
  bc.inflow.fluid.species.CO    = 0.1
  bc.inflow.fluid.species.H2O   = 0.2
  bc.inflow.fluid.species.CO2   = 0.1
  bc.inflow.fluid.species.O2    = 0.1
  bc.inflow.fluid.species.N2    = 0.2

   chemistry.integrator = StiffSolver::BackwardEuler
   chemistry.integrator.atol = 1.e-8
   chemistry.integrator.rtol = 1.e-9
  bc.outflow = po
  bc.outflow.fluid.pressure =  101325.