4.4.12.2. CHEMEQ Reference

A CHEMEQ reaction mechanism is specified in a nested block inside an Aria Material. The CHEMEQ block should generally define the following things:

Species names involved in the mechanism

Reactions in the mechanism

Species phases (only has an effect when integrating with a pressurization zone)

Concentration/energy release units

Begin Aria Material ReactingFoo

  # ...

  BEGIN PARAMETERS FOR CHEMEQ MODEL Foo
    species names  are Epoxy     CO2  CH4
    species phases are Condensed Gas  Gas
    Energy Release Units are Per Unit Mass

    Begin Reaction R1
      Reaction is Epoxy -> CH4 + CO2
      Rate Function = Arrhenius A = 1e9 Ea = 25000 R = 8.314
      Concentration Function = Standard mu = automatic
      Heat of Reaction = -1e8
    End

    Concentration Units = Mass Fractions
  END
End

The CHEMEQ model creates an ODE at each integration point that is solved to advance species concentrations and used to calculate an energy source term. The CHEMEQ model is not solved by default unless you include it in an energy source term. Using the chemeq_heating source term automatically locates the CHEMEQ model on the specified block and activates it.

EQ Energy for temperature on block_1 using Q1 with Mass Diff Src
Source for Energy on block_1 = chemeq_heating

Given a reaction rate $\omega_r$ and stoichiometric coefficient $\nu_{i,r}$ for species i in reaction r (negative for reactants, positive for products), the time rate of change of the concentration of species i ( $c_i$ ) is

$\frac{d c_i}{d t} = \sum_r \omega_r \nu_{i,r}$

At each nonlinear iteration, CHEMEQ solves this ODE from $t^n$ to $t^{n+1}$ at the current temperature to determine the total heat release rate, $\dot{q}$ , which is calculated using the start and end point of this integration and the per-species formation enthalpies ( $H_{f,i}$ ). Optionally, this may be multiplied by density depending on your reaction units.

$\dot{q} = \sum_i H_{f,i} \frac{c_i^{n+1} - c_i^n}{\Delta t}$

For a complete list of CHEMEQ commands, see the CHEMEQ command reference.

4.4.12.2.1. Defining Reactions

Reactions are defined in one or more Begin Reaction blocks inside the CHEMEQ model block. There is a legacy syntax for specifying reactions by supplying matrices of coefficients (described in Legacy Reaction Specification), but that approach is no longer recommended.

Each reaction is defined by an overall rate that is the product of one or more terms, a stoichiometric relationship, and a heat release. The syntax for specifying these reaction blocks is described in Reaction Blocks.

4.4.12.2.2. CHEMEQ Solver & Initial Conditions

Each CHEMEQ model requires a CHEMEQ Solver block that refers to it by name. For example, for the CHEMEQ model Foo shown earlier you would need to provide the following solver block in the Aria region scope

# Sierra > Aria Procedure > Aria Region
Begin CHEMEQ Solver for Foo
  ODE SOLVER = CVODE ADAMS 12 FUNCTIONAL
  species Epoxy = 1.0
  species CO2   = 0.0
  species CH4   = 0.0
End

The CHEMEQ Solver block specifies the ODE solver to use and the initial conditions for all species in the specified mechanism. For most mechanisms CVODE ADAMS 12 FUNCTIONAL is a good option. If the chemistry model is particularly stiff, using CVODE BDF 5 NEWTON may be more robust, but it is also slightly slower.

The initial concentrations provided must be in the same units assumed in the CHEMEQ model - Aria does not do any unit conversions for these.

4.4.12.2.3. Activation & Deactivation

The CHEMEQ Solver block can also be used to specify activation and deactivation conditions for the CHEMEQ mechanism (to have it turn on or off at specific temperature, temperature change rates or a more generic criteria).

# Sierra > Aria Procedure > Aria Region
Begin CHEMEQ Solver for Foo
  ODE SOLVER = CVODE ADAMS 12 FUNCTIONAL
  species Epoxy = 1.0
  species CO2   = 0.0
  species CH4   = 0.0

  Activation temperature = 350.0
  Deactivation Temperature = 600.0 terminate
End

Once the temperature hits the specified deactivation condition, the CHEMEQ solve stops and the simulation either stops or continues, depending on the option selected. The deactivation mechanism looks for a nodal value that reaches the specified deactivation condition (temperature, temperature rate or a generic condition). An example for a generic condition is to measure the reaction rate (which is correlated with the stiffness of the underlying chemistry ODE). The following steps show how to set this generic condition up:

The Aria material model that contains the targeted CHEMEQ reaction also needs to define a user expression for deactivation. Because we plan to use the reaction rate, a SPECIES_PRODUCTION model needs to be declared to define the time rate of change of the targeted species:

Begin Aria Material EM1
  ...
  #This user expression is needed for the generic deactivation condition
  user expression = scalar_string_function user_tag = "deactivation_fcn" f = "abs(epoxy_reaction_rate) > 1e-3"
  ...
  #This species_production model is only needed to define our reaction rate
  SPECIES_PRODUCTION of Epoxy = From_Time_Rate_Of_Change
  ...
End

In order to make epoxy_reaction_rate a usable expression, we need to add a postprocess command within the equation system (or Aria region)

Begin aria region myRegion
  ...
  Begin equation system Energy
    ...
    Postprocess species_production of Epoxy on all_blocks as epoxy_reaction_rate
    ...
  End
  ...
End

Now that epoxy_reaction_rate is properly defined for the user expression deactivation_fcn, we can use this within our CHEMEQ Solver block to define a generic deactivation condition. Notice that in the following block, we define 3 methods of deactivation. This is valid, and the first criteria to be reached will deactivate the CHEMEQ model:

# Sierra > Aria Procedure > Aria Region
Begin CHEMEQ Solver for Foo
  ODE SOLVER = CVODE ADAMS 12 FUNCTIONAL
  species Epoxy = 1.0
  species CO2   = 0.0
  species CH4   = 0.0

  Activation temperature = 350.0
  Deactivation Temperature = 600.0 terminate
  Deactivation Temperature rate = 20 terminate
  Deactivation expression = "deactivation_fcn" terminate
End

The deactivation option is often used for energetic materials when the reaction rates start to reach thermal runaway conditions. At that point, the cost of solving the chemistry model becomes prohibitive and it often makes sense to switch to a reduced model for the remainder of the burn. The two common options used for this are level set burns and uniform heat release.

To use a uniform heat release after deactivation, you need to specify Post Deactivation Release Time and Post Deactivation Heat Release to control how long the source should be active for and how much heat to release. When using a pressurization zone you should also specify Post Deactivation Gas Release (in moles per volume) to produce an appropriate pressurization.

In the example below, there will be a uniform energy source term of 1e9 W/m³ for 10 seconds after the CHEMEQ model deactivates.

Begin CHEMEQ Solver for Foo
  # ...

  Post Deactivation Release Time = 10.0
  Post Deactivation Heat Release = 1e9
  Post Deactivation Gas Release = 1.0

  Deactivation temperature = 550.0 continue

  # ...
End

4.4.12.2.4. Reaction Units

Some of the burden of managing the units properly for CHEMEQ falls to the user. There are currently two commands required to indicate the units being used in CHEMEQ: Energy Release Units and Concentration Units.

When you use the Energy Release Units are Per Unit Mass command, the total heat release term is multiplied by density (from the material model definition) and is calculated as

$\dot{q} = \rho \sum_i H_{f,i} \frac{c_i^{n+1} - c_i^n}{\Delta t}$

where $H_{f,i}$ is formation enthalpy for species i (either provided or computed from reaction heat releases). When you use Energy Release Units are Per Unit Volume the energy source term does not include density, and is calculated as

$\dot{q} = \sum_i H_{f,i} \frac{c_i^{n+1} - c_i^n}{\Delta t}$

This choice is directly related to your choice of concentration units in your chemical mechanism. If the concentration units are mass fractions, you should use the Energy Release Units are Per Unit Mass option. If the concentration units are mass densities, you should use the Energy Release Units are Per Unit Volume option.

There is also a Concentration Units = X command that is required in order to extract CHEMEQ concentrations out into material properties. The options for X include mass fractions, mole fractions, density, or molar.

4.4.12.2.5. Post-Processing

Sometimes you may want to extract species concentrations from a CHEMEQ model to use in material properties or to post-process. CHEMEQ creates element fields of its species concentrations that are available for output, but are not easily usable in most Aria post-processors.

The best way to use these quantities is by using the FROM_CHEMEQ model to pull them into the appropriate material property in the Aria Material block. You must have specified the concentration units using the Concentration Units command in the CHEMEQ model in order to do this. No unit conversions are done, but using an improper conversion will result in an error. For example, if the concentration units are mass fractions, you can only use FROM_CHEMEQ with the mass fraction property.

For example, to extract the species mass fractions from CHEMEQ and use them to define thermal conductivity you could do

Begin Aria Material ReactingFoo

  # ...

  Species Names = Cs CO2 O2
  Mass Fraction Cs  = From_Chemeq
  Mass Fraction CO2 = From_Chemeq
  Mass Fraction O2  = From_Chemeq

  Thermal Conductivity = Scalar_String_Function f = "4 + 0.1*Cs"

  BEGIN PARAMETERS FOR CHEMEQ MODEL Foo
    species names are Cs CO2 O2

    # ...

    Concentration Units = Mass_Fractions
  END
End

To use them in a post-processor after defining them with FROM_CHEMEQ you can use the standard syntax

Postprocess average of expression mass_fraction of CO2 on block_1 as avg_co2

4.4.12.2.6. Legacy Reaction Specification

Warning

This section is provided for reference to help understand and convert legacy CHEMEQ models. We do not recommend using this syntax for new models.

Instead of using reaction blocks, older CHEMEQ models provided a grid of concentration exponents and stoichiometric coefficients along with arrays of activation energies, pre-exponential terms, and heat releases.

For example, the following legacy input

BEGIN PARAMETERS FOR CHEMEQ MODEL DEMO
  number of reactions is 3
  species names  are A B C D
  Energy Release Units are Per Unit Mass

  # Define Arrhenius terms and heat release
  Steric Coefficients are        0.0    0.0    0.0
  Log Preexponential Factors are 48.7   37.3   28.1
  Activation Energies are        200000 300000 400000
  Energy Releases are           -1e6    2e7    3e5

  # Define the reaction rate concentration dependence
  Concentration Exponents for A  are 1.0 0.0 0.0
  Concentration Exponents for B  are 0.0 1.0 0.0
  Concentration Exponents for C  are 0.0 0.0 2.0
  Concentration Exponents for D  are 0.0 0.0 0.0

  # Define the reactions
  #  A -> B
  #  B -> C
  #  A + C -> 2D
  Stoichiometric coefficients for A are -1.0  0.0 -1.0
  Stoichiometric coefficients for B are  1.0 -1.0  0.0
  Stoichiometric coefficients for C are  0.0  1.0 -1.0
  Stoichiometric coefficients for D are  0.0  0.0  2.0

  Concentration Units are Mass Fractions
END

can be converted to a block style below

BEGIN PARAMETERS FOR CHEMEQ MODEL DEMO
  species names  are A B C D
  Energy Release Units are Per Unit Mass

  Begin Reaction R1
    Reaction is A -> B
    Rate Function = Arrhenius A = {exp(48.7)} Ea = 200000 R = 8.314
    Concentration Function = Automatic
    Heat of Reaction = 1e6
  End

  Begin Reaction R2
    Reaction is B -> C
    Rate Function = Arrhenius A = {exp(37.3)} Ea = 300000 R = 8.314
    Concentration Function = Automatic
    Heat of Reaction = -2e7
  End

  Begin Reaction R3
    Reaction is A + C -> 2D
    Rate Function = Arrhenius A = {exp(28.1)} Ea = 400000 R = 8.314
    Concentration Function = C 2.0
    Heat of Reaction = -3e5
  End

  Concentration Units are Mass Fractions
END

Note

The sign convention for the legacy heat of reaction is the opposite of the block-style reaction convention.

The legacy style required all reactions to use the same basic rate function, while the block-style format provides significantly greater flexibility.