Reactive Infiltration

Table of Contents

• Framework
  • Representative Volume Element
  • Governing Equations
  • Implementation Details
• Model Building Blocks
  • Product and solid's inner radius
  • Product reaction rate
  • Hermite Smooth Step Function
• Complete Input File

The reactive infiltration physics module is a collection of objects serving as building blocks for composing model's describing reactive infiltration kinetics as a solid $(S)$ material is infiltrated by a liquid $(L)$ to create a product $(P)$. The framework and usage of the model is explained below, with both mathematical formulations and example input files.

Framework

Representative Volume Element

The infiltration process for a material considers a cylindrical Representative Volume Element (RVE) with radius $R$ and span $H$, depicting a solid capillary of radius $r_o$, corresponding to the initial porosity of ${\phi }_0$. The mathematical description of the model uses cylindrical coordinates and assumes axisymmetry around $r=0$.

As the liquid enters the cylinder, it reacts with a solid to form a product with thickness ${\delta }_P$, as shown schematically in the figure below.

The (a) top and (b) cross-section view of the Representative Volume Element (RVE) of a given state of the reactive infiltration process, depicting the outer solid walls, the infiltrated liquid and the product from the chemical reaction

Throughout this process, the solid radius, $r_o$, decreases while the product thickness, ${\delta }_P$, increases. Let $r_i$ denote the inner radius of the product. In addition, define ${\alpha }_i$ (in units of mole per volume) as the amount of substance in the RVE, and ${\mathrm{\Omega }}_i={\displaystyle \frac{M_i}{{\rho }_i}}$ as the molar volume of a material with molar mass $M_i$ (in units of amu, or gram per mole) and density ${\rho }_i$ (mass per volume), with subscripts taking $L$, $S$, and $P$, respectively.

The volume fraction, ${\phi }_i$ of each material is then

$${\phi }_i={\alpha }_i{\mathrm{\Omega }}_i$$

and the RVE porosity (void) is

$${\phi }_v=1-{\phi }_L-{\phi }_P-{\phi }_S$$

Key assumptions made throughout the derivation are:

• Liquid remains liquid over the entire course of reaction, i.e., no phase change.
• The reaction is irreversible.
• Formation of the initial product layer is immediate once liquid comes into contact with solid.
• Once product is formed, reaction rate is primarily limited by the diffusion of liquid through product to the product-solid interface, with a diffusion coefficient $D$.
• The product wall thickness remains uniform during the infiltration.
• The only reaction is between the liquid and the solid.

The following nondimensionalization is applied to the constitutive model for the reactive infiltration process.

$$\begin{array}{rl}& {\overline{\delta }}_P=\frac{{\delta }_P}{R},{\overline{r}}_o={\displaystyle \frac{r_o}{R}}=\sqrt{1-{\phi }_S},{\overline{r}}_i={\displaystyle \frac{r_i}{R}}=\sqrt{1-{\phi }_S-{\phi }_P}.\end{array}$$

The complete state of the RVE is denoted by the tuple $({\phi }_L,{\phi }_S,{\phi }_P)$, with ${\alpha }_L$ as the prescribed force.

Mathematically, it is possible that ${\phi }_L+{\phi }_P+{\phi }_S\ge 1$. Physically, this implies "overflow", aka the prescribed ${\alpha }_L$ is larger than the available voids. Care must be taken at the macroscopic model to avoid or resolve this issue.

Governing Equations

The initial-value problem (IVP) corresponding to the constitutive model is

$$\mathbf{r}=\left\{\begin{array}{c}{r}_{{\phi }_P}\\ r_{{\phi }_S}\end{array}\right\}=\mathbf{0},$$

where $r_{{\phi }_i}$ is the residual corresponding to the volume fraction of the liquid, solid and product. The residuals are defined below. Note an implicit backward-Euler time intergration scheme is used.

• Production rate of the product, dictates by the diffusion of the liquid through the product to react with the solid at the product-solid interface

  $$r_{{\phi }_P}={\displaystyle \frac{{\alpha }_{P,n+1}-{\alpha }_{P,n}}{t_{n+1}-t_n}}-{\displaystyle \frac{2}{{\mathrm{\Omega }}_L}}{\displaystyle \frac{D}{l_c^2}}{\displaystyle \frac{{\overline{r}}_o+{\overline{r}}_i}{{\overline{r}}_o-{\overline{r}}_i}}{R}_L{R}_S$$

  Here $R_L({\phi }_L),R_S({\phi }_S)\in [0,1]$ are the reactivity of liquid and solid, represented by a step function of the void fraction, ${\phi }_i$.
• Reaction rate of the solid, with a (solid, product) reaction coefficient, $(k_S,k_P)$ respectively

  $$r_{{\phi }_S}={\displaystyle \frac{{\alpha }_{S,n+1}-{\alpha }_{S,n}}{t_{n+1}-t_n}}+{\displaystyle \frac{k_P}{k_S}}{\displaystyle \frac{{\alpha }_{P,n+1}-{\alpha }_{P,n}}{t_{n+1}-t_n}}$$

Implementation Details

The above governing equations are fairly general for a wide range of reactive infiltration systems within this framework given the chemical reaction kinetics with appropriate reaction coefficient ( $k_i$) and the corresponding molar volume ${\mathrm{\Omega }}_i$ of the materials.

Example for determining $k_i$ and ${\mathrm{\Omega }}_i$: Consider a reactive infiltration process of liquid (chemical formula $L_x$) and solid ( $S_y$), with the molar mass $M_L$ and $M_S$. The chemical reaction creates a product with chemical formula $L_m{S}_n$, with reaction coefficients $k_S$ and $k_P$.

$$L_x+k_S{S}_y\to k_P{L}_m{S}_n$$

For stoichiometric balance, the reaction coefficients are:

$$\begin{array}{r}{k}_S={\displaystyle \frac{xn}{ym}},k_P={\displaystyle \frac{x}{m}}\end{array}$$

The densities of the liquid, solid, and product are respectively ${\rho }_L,{\rho }_S,{\rho }_P$.

Then, the molar volume is:

$${\mathrm{\Omega }}_L={\displaystyle \frac{x{M}_L}{{\rho }_L}}.{\mathrm{\Omega }}_S={\displaystyle \frac{y{M}_S}{{\rho }_S}},{\mathrm{\Omega }}_P={\displaystyle \frac{m{M}_L+n{M}_S}{{\rho }_P}}.$$

The instantaneous balance also implies that

$${\dot{\alpha }}_S={\displaystyle \frac{k_P}{k_S}}{\dot{\alpha }}_P,$$

Since this constitutive model considers ${\alpha }_L$ as the forcing function for the IVP problem, ${\dot{\alpha }}_L\ne k_P{\dot{\alpha }}_P$.

The following tables summarize the relationship between the mathematical expressions and NEML2 models.

Expression	Model building block	Syntax
${\overline{r}}_i=\sqrt{1-{\phi }_S-{\phi }_P};{\overline{r}}_o=\sqrt{1-{\phi }_S}$	Product and solid's inner radius	ProductGeometry
${\displaystyle \frac{2}{{\mathrm{\Omega }}_L}}{\displaystyle \frac{D}{l_c^2}}{\displaystyle \frac{{\overline{r}}_o+{\overline{r}}_i}{{\overline{r}}_o-{\overline{r}}_i}}{R}_L{R}_S$	Product reaction rate	DiffusionLimitedReaction
$R_L;R_S$	Step Function	HermiteSmoothStep
${\phi }_i={\alpha }_i{\mathrm{\Omega }}_i$	—	Linear Combination

Residual Expression	Syntax
$r_{{\phi }_P}={\displaystyle \frac{{\alpha }_{P,n+1}-{\alpha }_{P,n}}{t_{n+1}-t_n}}-{\displaystyle \frac{2}{{\mathrm{\Omega }}_L}}{\displaystyle \frac{D}{l_c^2}}{\displaystyle \frac{{\overline{r}}_o+{\overline{r}}_i}{{\overline{r}}_o-{\overline{r}}_i}}{R}_L{R}_S$	Linear Combination, Variable Rate
$r_{{\phi }_S}={\displaystyle \frac{{\alpha }_{S,n+1}-{\alpha }_{S,n}}{t_{n+1}-t_n}}+{\displaystyle \frac{k_P}{k_S}}{\displaystyle \frac{{\alpha }_{P,n+1}-{\alpha }_{P,n}}{t_{n+1}-t_n}}$	Linear Combination, Variable Rate

Model Building Blocks

The RVE keeps track of the volume fraction of the liquid, solid and product, ${\phi }_L,{\phi }_S,{\phi }_P$. ${\phi }_L$ directly corresponds to the prescribed liquid substance ${\alpha }_L$ where ${\phi }_{S}and{\phi }_P$ are obtained from solving the IVP problem with ImplicitUpdate.

Product and solid's inner radius

Model object ProductGeometry

This modelcalculates the product and solid's inner radius, $r_i,r_o$. Note that, when there are presence of the product, the solid' inner radius is the outer radius of the product.

$${\overline{r}}_o=\sqrt{1-{\phi }_S}{\overline{r}}_i=\sqrt{1-{\phi }_S-{\phi }_P}$$

Example input file that defines the product and solid's inner radius

[Models]
  [model]
    type = ProductGeometry
    solid_fraction = 'state/phi_s'
    product_fraction = 'state/phi_p'
    inner_radius = 'state/ri'
    outer_radius = 'state/ro'
  []
[]

Product reaction rate

Model object DiffusionLimitedReaction

The reaction rates dictates the creation of the product. Within our assumptions, the growth of the product is controlled by the diffusion of the liquid through the product to react with the solid at the product-solid interface. Considered the same RVE as in figure above, with $c_P$ (units of moles) as the concentration of the product in the RVE. Diffusion equation yields

$${\dot{c}}_P=\mathrm{\nabla }\cdot D\mathrm{\nabla }{c}_P$$

$D$ is the diffusion coeficient (untis of area per time) of liquid through the product. Integrating through the product's volume and applying Green's theorem,

$$\dot{C}=D{\displaystyle \frac{C_{liquid-product}-C_{product-solid}}{{\delta }_P}}{A}_{enclosed}$$

Here,

$$\begin{array}{rl}& \dot{C}={\dot{\alpha }}_P{V}_{RVE}{C}_{product-solid}=0C_{liquid-product}={\displaystyle \frac{{\dot{\alpha }}_L{V}_{RVE}}{{\phi }_L{V}_{RVE}}}={\displaystyle \frac{1}{{\mathrm{\Omega }}_L}}\\ & A_{enclosed}=2\pi (r_o+r_i)H{\delta }_P=r_o-r_i\end{array}$$

The step function $R_L({\alpha }_L)$, $R_S({\alpha }_S)$ as reactivity is introduced to supressed the production rate when either the solid or the liquid is fully consumed. Finally, in normalized variables, and let the $R=l_c$ as the characterisitc length scales of the pores features. The product reaction rate is then,

$${\dot{\alpha }}_P={\displaystyle \frac{2}{{\mathrm{\Omega }}_L}}{\displaystyle \frac{D}{l_c^2}}{\displaystyle \frac{{\overline{r}}_o+{\overline{r}}_i}{{\overline{r}}_o-{\overline{r}}_i}}{R}_L{R}_S$$

Example input file that defines the product reaction rate

[Models]
  [model]
    type = DiffusionLimitedReaction
    diffusion_coefficient = 1e-3
    molar_volume = 1
    product_inner_radius = 'state/ri'
    solid_inner_radius = 'state/ro'
    liquid_reactivity = 'state/R_l'
    solid_reactivity = 'state/R_s'
    reaction_rate = 'state/alpha_rate'
  []
[]

Hermite Smooth Step Function

Model object HermiteSmoothStep

A cubic Hermite interpolation after clamping is used as a smooth step function:

$$R(x_r)=\{\begin{array}{ll}0,& x_r<0\\ 3x_r^2-2x_r^3,& 0\le x_r\le 1\\ 1,& x_r>1\end{array}$$

where

$$x_r={\displaystyle \frac{x-x_l}{x_u-x_l}}$$

Here, $x_l,x_u$ denote the lower and upper bound of the smoothen transition regime.

The Hermite Function is used in the reactive infiltration for the liquid and solid reactivity $R_i$

Example input file that use the smooth step function

[Models]
  [model]
    type = HermiteSmoothStep
    argument = 'state/foo'
    value = 'state/bar'
    lower_bound = '0'
    upper_bound = '0.05'
  []
[]

Complete Input File

A complete example input file for the reactive infiltration is shown below with the appropriate model composition.

ntime = 200
D = 1e-6
omega_Si = 12.0
omega_C = 5.3
omega_SiC = 12.5
chem_ratio = 1.0
 
oCm1 = 0.18867924528 # 1/Omega_C
oSiCm1 = 0.08 # 1/Omega_SiC
 
[Tensors]
  [times]
    type = LinspaceScalar
    start = 0
    end = 1e4
    nstep = ${ntime}
  []
  [alpha]
    type = FullScalar
    batch_shape = '(${ntime})'
    value = 0.01
  []
[]
 
[Drivers]
  [driver]
    type = InfiltrationDriver
    model = 'model'
    prescribed_time = 'times'
    time = 'forces/t'
    prescribed_liquid_species_concentration = 'alpha'
    liquid_species_concentration = 'forces/alpha'
    ic_Scalar_names = 'state/phi_P state/phi_S state/alpha_P state/alpha_S'
    ic_Scalar_values = '0 0.3 0 0.05660377358'
    save_as = 'result.pt'
    verbose = true
    show_input_axis = true
    show_output_axis = true
  []
  [regression]
    type = TransientRegression
    driver = 'driver'
    reference = 'gold/result.pt'
  []
[]
 
[Solvers]
  [newton]
    type = Newton
    verbose = true
  []
[]
 
[Models]
  [liquid_volume_fraction]
    type = ScalarLinearCombination
    from_var = 'forces/alpha'
    to_var = 'state/phi_L'
    coefficients = '${omega_Si}'
  []
  [outer_radius]
    type = ProductGeometry
    solid_fraction = 'state/phi_S'
    product_fraction = 'state/phi_P'
    inner_radius = 'state/ri'
    outer_radius = 'state/ro'
  []
  [liquid_reactivity]
    type = HermiteSmoothStep
    argument = 'state/phi_L'
    value = 'state/R_L'
    lower_bound = 0
    upper_bound = 0.1
  []
  [solid_reactivity]
    type = HermiteSmoothStep
    argument = 'state/phi_S'
    value = 'state/R_S'
    lower_bound = 0
    upper_bound = 0.1
  []
  [reaction_rate]
    type = DiffusionLimitedReaction
    diffusion_coefficient = '${D}'
    molar_volume = '${omega_Si}'
    product_inner_radius = 'state/ri'
    solid_inner_radius = 'state/ro'
    liquid_reactivity = 'state/R_L'
    solid_reactivity = 'state/R_S'
    reaction_rate = 'state/react'
  []
  [substance_product]
    type = ScalarLinearCombination
    from_var = 'state/phi_P'
    to_var = 'state/alpha_P'
    coefficients = '${oSiCm1}'
  []
  [product_rate]
    type = ScalarVariableRate
    variable = 'state/alpha_P'
    rate = 'state/adot_P'
    time = 'forces/t'
  []
  [substance_solid]
    type = ScalarLinearCombination
    from_var = 'state/phi_S'
    to_var = 'state/alpha_S'
    coefficients = '${oCm1}'
  []
  [solid_rate]
    type = ScalarVariableRate
    variable = 'state/alpha_S'
    rate = 'state/adot_S'
    time = 'forces/t'
  []
  ##############################################
  ### IVP
  ##############################################
  [residual_phiP]
    type = ScalarLinearCombination
    from_var = 'state/adot_P state/react'
    to_var = 'residual/phi_P'
    coefficients = '1.0 -1.0'
  []
  [residual_phiL]
    type = ScalarLinearCombination
    from_var = 'state/adot_P state/adot_S'
    to_var = 'residual/phi_S'
    coefficients = '1.0 ${chem_ratio}'
  []
  [model_residual]
    type = ComposedModel
    models = "residual_phiP residual_phiL
              liquid_volume_fraction outer_radius liquid_reactivity solid_reactivity
              reaction_rate substance_product product_rate substance_solid  solid_rate"
  []
  [model_update]
    type = ImplicitUpdate
    implicit_model = 'model_residual'
    solver = 'newton'
  []
  [substance_product_new]
    type = ScalarLinearCombination
    from_var = 'state/phi_P'
    to_var = 'state/alpha_P'
    coefficients = '${oSiCm1}'
  []
  [substance_solid_new]
    type = ScalarLinearCombination
    from_var = 'state/phi_S'
    to_var = 'state/alpha_S'
    coefficients = '${oCm1}'
  []
  [model]
    type = ComposedModel
    models = 'model_update liquid_volume_fraction substance_solid_new substance_product_new'
    additional_outputs = 'state/phi_P state/phi_S'
  []
[]