Chemical reactions¶

Overview¶

The chemical-reactions module is a collection of objects that serve as building blocks for composing models of macroscale chemical reactions — pyrolysis, reactive infiltration, solidification, and similar control-mass processes. The module provides empirical, homogenized expressions for the reaction rate \(\dot{\alpha}\) as a function of the conversion degree \(\alpha\), plus geometric/volumetric helpers that bridge the reaction’s control mass into a control-volume picture other physics modules (solid mechanics, porous flow, phase-field fracture) can consume.

The catalog carves the domain into two layers:

Reaction mechanisms — typed maps \(\alpha \mapsto \dot{\alpha}\) for a single reacting phase. AvramiErofeevNucleation, ContractingGeometry, and DiffusionLimitedReaction are the three concrete leaves; they share the canonical conversion_degree / reaction_rate variable names, so they can be swapped interchangeably in a composed model.
Geometry and volume bridges — kinematic primitives that turn the reaction state into the volumetric quantities other modules need. CylindricalChannelGeometry gives the dimensionless inner and outer radii of a cylindrical shrinking-core product; EffectiveVolume sums per-component mass-fraction / density contributions into the total composite volume so the reaction’s control mass can be tied into a control-volume framework.

Time integration, predictors, Arrhenius temperature dependence, linear combinations, and the implicit Newton solve are provided by NEML2’s shared primitives. The example below shows how the chemical-reactions types slot into that machinery.

Math¶

A reaction mechanism is a constitutive relation between the conversion degree \(\alpha \in [0, 1]\) and the reaction rate \(\dot{\alpha} = f(\alpha)\). The three mechanisms in this module are

(4)¶\[\begin{align} \dot{\alpha} &= k \, (1 - \alpha)^n &&\text{(contracting geometry)} \\ \dot{\alpha} &= k \, (1 - \alpha) \, \bigl(-\ln(1 - \alpha)\bigr)^n &&\text{(Avrami--Erofeev nucleation)} \\ \dot{\alpha} &= \frac{2 \, D \, R_l \, R_s}{\omega} \, \frac{r_o}{r_o - r_i + \delta} &&\text{(diffusion-limited)} \end{align}\]

where \(k\) is the reaction coefficient (typically Arrhenius, \(k = k_0 \exp(-Q / R T)\)), \(n\) is the reaction order, \(D\) is the diffusion coefficient of the rate-limiting species through a product layer of dimensionless inner / outer radii \((r_i, r_o)\), \(R_l, R_s \in [0, 1]\) are reactivities of the liquid and solid phases, \(\omega\) is the molar volume of the rate-limiting species, and \(\delta\) is a small “dummy thickness” that keeps the rate finite at the start of the reaction.

For a cylindrical shrinking-core geometry parametrized by the volume fractions \(\phi_s\) (solid) and \(\phi_p\) (product), the dimensionless radii are

\[ r_i = \sqrt{1 - \phi_s - \phi_p}, \qquad r_o = \sqrt{1 - \phi_s}. \]

To couple a control-mass reaction to a control-volume framework, the total composite volume is reconstructed from per-component mass fractions \(\omega_i\) and densities \(\rho_i\),

\[ V = \frac{M}{1 - \phi_o} \sum_i \frac{\omega_i}{\rho_i}, \]

where \(M\) is the reference mass and \(\phi_o \in [0, 1)\) is the open volume fraction accounting for material leaving the control mass (for example, gas escaping a pyrolyzing binder). When the input file omits \(\phi_o\), the prefactor reduces to \(M\).

Example: pyrolysis of a binder–particle composite¶

The pyrolysis regression test composes a contracting-geometry reaction with an Arrhenius reaction coefficient, implicit integration of \(\alpha\), and explicit forward-Euler integration of the four mass / volume fractions \((\omega_b, \omega_c, \omega_g, \varphi_o)\) tracking binder, char, gas, and open pore.

# neml2
nbatch = '(1)'
nstep = 100

# reaction mechanism
Y = 0.5835777126099713
# yield
n = 1.0
# reaction order
k0 = 0.04210147513030456
# reaction rate coefficient
Q = 21191.61425138572
# J/mol
R = 8.31446261815324
# J/K/mol

# initial mass fraction
wc0 = 0
# char (residue)
wb0 = 1
# binder (precursor)

# control volume
mu = 0.2
zeta = 0.05

[Tensors]
  [yield_surface]
    type = Python
    expr = 'Scalar(torch.tensor(${Y}, dtype=torch.float64))'
  []
  [endtime]
    type = Python
    expr = 'Scalar([2700.0])'
  []
  [times]
    type = Python
    expr = 'Scalar(torch.linspace(0.0, 2700.0, ${nstep}, dtype=torch.float64).reshape(${nstep}, 1))'
  []
  [T]
    type = Python
    expr = 'Scalar.linspace(300.0, 1500.0, ${nstep})'
  []
[]

[Drivers]
  [driver]
    type = TransientDriver
    model = 'model'
    prescribed_time = 'times'
    force_Scalar_names = 'T'
    force_Scalar_values = 'T'
    ic_Scalar_names = 'wb wc'
    ic_Scalar_values = '${wb0} ${wc0}'
    save_as = 'result.pt'
  []
  [regression]
    type = TransientRegression
    driver = 'driver'
    reference = 'gold/result.pt'
  []
[]

[Models]
  [reaction_coef]
    type = ArrheniusParameter
    reference_value = '${k0}'
    activation_energy = '${Q}'
    ideal_gas_constant = '${R}'
    temperature = 'T'
  []
  [reaction_rate]
    type = ContractingGeometry
    coef = 'reaction_coef'
    order = '${n}'
    conversion_degree = 'alpha'
    reaction_rate = 'alpha_rate'
  []
  [reaction_ode]
    type = ScalarBackwardEulerTimeIntegration
    variable = 'alpha'
  []
  [reaction]
    type = ComposedModel
    models = 'reaction_coef reaction_rate reaction_ode'
  []
[]

[EquationSystems]
  [eq_sys]
    type = NonlinearSystem
    model = 'reaction'
    unknowns = 'alpha'
  []
[]

[Solvers]
  [newton]
    type = Newton
    linear_solver = 'lu'
  []
  [lu]
    type = DenseLU
  []
[]

[Models]
  [predictor]
    type = ConstantExtrapolationPredictor
    unknowns_Scalar = 'alpha'
  []
  [solve_reaction]
    type = ImplicitUpdate
    equation_system = 'eq_sys'
    solver = 'newton'
    predictor = 'predictor'
  []
  [binder_rate]
    type = ScalarLinearCombination
    from = 'alpha_rate'
    weights = '-1'
    to = 'wb_rate'
  []
  [char_rate]
    type = ScalarLinearCombination
    from = 'alpha_rate'
    weights = '${Y}'
    to = 'wc_rate'
  []
  [gas_rate]
    type = ScalarLinearCombination
    from = 'wb_rate wc_rate'
    weights = '-0.2 -0.2'
    to = 'wg_rate'
  []
  [open_pore_rate]
    type = ScalarLinearCombination
    from = 'alpha_rate'
    weights = '${zeta}'
    to = 'phio_rate'
  []
  [binder]
    type = ScalarForwardEulerTimeIntegration
    variable = 'wb'
  []
  [char]
    type = ScalarForwardEulerTimeIntegration
    variable = 'wc'
  []
  [gas]
    type = ScalarForwardEulerTimeIntegration
    variable = 'wg'
  []
  [open_pore]
    type = ScalarForwardEulerTimeIntegration
    variable = 'phio'
  []
  [model]
    type = ComposedModel
    models = 'solve_reaction reaction_rate
              binder_rate char_rate gas_rate open_pore_rate
              binder char gas open_pore'
    additional_outputs = 'alpha'
  []
[]

Explanation¶

The reaction kinetics block — reaction_coef, reaction_rate, reaction_ode composed as reaction — is the entire chemical-reactions contribution:

ArrheniusParameter (reaction_coef) computes the temperature-dependent rate constant \(k(T) = k_0 \exp(-Q/RT)\) from the T force.
ContractingGeometry (reaction_rate) maps the current \(\alpha\) to \(\dot{\alpha} = k (1 - \alpha)^n\) using that coefficient as a nonlinear parameter and the constant order \(n\).
ScalarBackwardEulerTimeIntegration (reaction_ode) writes the residual \(\alpha_{n+1} - \alpha_n - \Delta t \, \dot{\alpha} = 0\) that the Newton solver consumes via the eq_sys ImplicitUpdate block.

The remaining blocks downstream of solve_reaction are conservation relations and time integration from NEML2’s shared primitives: the four ScalarLinearCombination blocks (binder_rate, char_rate, gas_rate, open_pore_rate) implement conservation relations such as \(\dot{\omega}_c = Y \, \dot{\alpha}\) (with \(Y\) the char yield) and \(\dot{\omega}_g = -\mu (\dot{\omega}_b + \dot{\omega}_c)\) with \(\mu\) the trapped-gas fraction, and four ScalarForwardEulerTimeIntegration blocks advance the mass / volume fractions in time. The ConstantExtrapolationPredictor seeds the Newton iterate at each step.

Variable names flow through the composition by string match: the mechanism produces alpha_rate, which the linear combinations consume and rename into wb_rate, wc_rate, wg_rate, phio_rate; the time integrators then pair those rates with the state variables wb, wc, wg, phio; the driver’s initial conditions initialize wb and wc explicitly, while wg and phio start from zero. The final model ComposedModel glues the implicit solve, the rate computations, and the explicit integrations into one forward operator that the TransientDriver walks over the prescribed times history. additional_outputs = 'alpha' keeps the implicit unknown in the result file even though it is not consumed downstream.

The companion reactive-infiltration regression scenario under tests/regression/chemical_reactions/reactive_infiltration/ follows a similar pattern but swaps in CylindricalChannelGeometry + DiffusionLimitedReaction, uses HermiteSmoothStep to give the liquid and solid reactivities \(R_l(\phi_L)\) and \(R_s(\phi_S)\) a smooth on/off behaviour, and Newton-solves on the volume fractions \((\phi_P, \phi_S)\) instead of \(\alpha\).

Chemical reactions¶

Overview¶

Math¶

Example: pyrolysis of a binder–particle composite¶

Explanation¶

See also¶