diff --git a/docs/source_docs/index.rst b/docs/source_docs/index.rst index ba50691c2607baec0b64601024c4b3ae75ad62df..fd5c22545597b6741cd81e29740549fcff64dea4 100644 --- a/docs/source_docs/index.rst +++ b/docs/source_docs/index.rst @@ -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 diff --git a/docs/source_docs/references/Chemistry_Chapter.rst b/docs/source_docs/references/Chemistry_Chapter.rst new file mode 100644 index 0000000000000000000000000000000000000000..c167a4092592ecb36d24b83059d2ce30b3d2d7eb --- /dev/null +++ b/docs/source_docs/references/Chemistry_Chapter.rst @@ -0,0 +1,239 @@ +.. _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 + 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 + 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 + 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 diff --git a/docs/source_docs/user_guide/inputs/chemical_reactions.rst b/docs/source_docs/user_guide/inputs/chemical_reactions.rst index 6e67ba877274c0dd07b4f3b4767db4b85f615537..b72dd3d92015d5dd5c836d50341182ce37d07440 100644 --- a/docs/source_docs/user_guide/inputs/chemical_reactions.rst +++ b/docs/source_docs/user_guide/inputs/chemical_reactions.rst @@ -1,24 +1,20 @@ 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.