Stability of chemical reaction fronts in the vicinity of a blocking state

ВЕСТНИК ПНИПУ. МЕХАНИКА № 3, 2019PNRPU MECHANICS BULLETIN

Various experimental and theoretical observations of the propagation of chemical reaction fronts in deformable solids show the influence of mechanical stresses on the reaction front kinetics (see, e.g., papers on silicon oxidation [1-3], on silicon lithiation [4-9] or on intermetallic compound growth [10-12]). Chemical reaction retardation and even blocking by mechanical stresses were experimentally observed in e.g. [13-18]. A special attention was paid to the influence of mechanical stresses on the stability of the growing interface [19-23]. These problems are actual in the solid electrolyte interfaces modelling for the Li-ion batteries [24, 25].

In general, two major processes control the propagation of a sharp chemical reaction front: (i) difusion of a reactant in the body undergoing a chemical reaction, (ii) consumption of the diffusive reactant by a chemical reaction at the reaction front. Thus, stresses may affect the reaction front propagation via the influence on diffusion or/and the influence on the reaction rate.

In order to couple chemical reaction rate with the mechanical stresses, some models include additional stress dependent cross-effect terms in the diffusion flux which appear basing on generalized expressions of stressdependent scalar chemical potentials [5-7,26-30].

But the velocity of the reaction front may be controlled rather by the reaction rate than by the diffusion (see, e.g., [31, 32] and references therein). In this case the influence of stresses on the reaction rate becomes important. In classical physical chemistry the reaction rate is determined by a scalar chemical affinity equal to a combination of scalar chemical potentials of the reaction constituents [33]. Scalar chemical potentials were derived for the case of phase transformations in gases and liquids where stresses were reduced to a scalar pressure [34].

In the last decades of XXth century it was recognized that in the case of solid phases a chemical potential was a tensor. Tensorial nature followed from the fact that the equilibrium was considered not just in a point of the interface but at oriented surface element [35–38]. Later a tensorial nature of the chemical potential was discussed in [39, 40]. These results could be considered as a prelude to the tenso-rial chemical affinity, and in [41–43] it was shown that in the case of a propagating reaction front the driving force is equal to the normal component of the chemical affinity tensor that, in turn, equals to the combination of the chemical potential tensors of solid constituents and the chemical potential of the diffusive constituent. By this approach, mechanical stresses affect the reaction rate via the chemical affinity tensor. A kinetic equation for the propagating chemical reaction front was formulated in the form of the dependence of the normal component of the reaction front velocity on the normal component of the chemical affinity tensor.

The approach was approved by the solution of a number of boundary value problems for elastic and inelastic solid undergoing chemical reactions and, in particular, demonstrated that stresses can accelerate, retard and even block the reaction (see [44,45] and reference therein). The latter gave reasons to construct so called forbidden zones in strain space formed by strains at which the direct reaction cannot go [46, 47]. The boundary of the forbidden zone is formed by strains at which the reaction front corresponds to chemical equilibrium, and in the present paper we focus on the reaction front behavior when it approaches to the equilibrium. We start with examinations of the stability of the equilibrium reaction front. It should be noted that the stability of the interface between two phases during the phase transformation was previously studied by various approaches, see e.g., [48–57].

In the present paper we follow the procedure of the kinetic stability analysis developed for the case of phase transformations in [56, 57]. The procedure is based on the examinations of the behavior of the front perturbations. If the perturbations grow due to the kinetic equation, then the unperturbed state is unstable. In the case of the equilibrium chemical reaction front this is reduced to the linear stability analysis that includes perturbed jump conditions, diffusion equation and reaction front evolution equation. Far from equilibrium we use numerical simulations.

The paper is organized as follows. In Section 2 we briefly describe the approach and formulate the kinetic equation for the chemical reaction front propagation based on the chemical affinity tensor. In Section 3 results of the numerical simulations of the reaction front propagation are presented. We analyze the interface kinetics for the cases of stable and unstable configurations for the axial symmetry problem of a solid cylinder undergoing a chemical reaction. We explore the transition of the interface motion to the unstable mode. Finally, in Section 4, we study the stresses redistributions due to stability loss which may lead to intensive plastic deformations or fracture.

1.    Chemical affinity tensor and kinetic equation

In general, a chemical reaction between the solid, В – , and a diffusive, В * , constituents can be written as follows:

n - B _ + n * B * ^ n + B ₊ , (1)

where В ₊ is the product of the reaction and the n ’s refer to the corresponding stochiometric coefficients. As shown in Fig. 1, the В * constituent diffuses through the material В + and reacts with the material В – at the sharp interface.

Fig. 1. Schematic representation of the chemical reaction process in deformable solids

Our further considerations of the kinetics of a chemical reaction front are based on the concept of a chemical affinity tensor A, [41–43]. According to the concept, the reaction rate at the oriented surface element with the normal N can be expressed through the normal component ANN of the chemical affinity tensor as го = k* c

1 - exp

where k _* is the kinetic coefficient, c is the diffusive constituent molar concentration, T is the temperature, R is the universal gas constant. Note that formula (2) just generalizes the kinetic equation formulated for the case of a scalar chemical potential [58]. Then the normal component of the reaction front velocity, W , can be expressed from the mass balance via the reaction rate as:

W =

nM

---------го ,

P -

where ρ_– and M _– are the mass density and molar mass of the constituent B _–. Chemical equilibrium, i.e., front propagation blocking, takes place if stresses at the reaction front are such that

^ aw = °-                      ⁽⁴⁾

Further we consider the case of linear elastic solid constituents, so the stresses σ and strains in materials B – and B + are related by Hooke’s law:

C - = C - : e _ , a ₊= C ₊ : ( e + - E ^ch ) ,           (5)

where C – and C + are the stiffness tensors of solid constituents, E ^ch is the chemical transformation strain. We take the chemical potential of the diffusive constituent in the form

c

M * p * = n * ( T ) + RT In —,            (6)

c *

where M* and p* are the molar mass and specific chemical potential of the constituent B„ , c* and n * ( T ) are reference values of the concentration and chemical energy. Then, if to neglect the pressure produced by the diffusive constituent at the reaction front compared with stresses produced in solid constituents by the transformation strain and external loading, in a quasistatic case, the expression of A takes the form (see details in [41–43]):

nMc

A w =---- X + n * RT n —,

P-

where

X = Y +1C- : E- - 1C+ : (E+ - Ech) + C± : (E+ - E-), and у is the temperature dependent combination of chemical energies of the species (free energies of solid constitu- ents in stress-free states and reference chemical potential of the diffusive constituent).

Stresses and strains in Eq. (8) as well as concentration of the diffusive component in the second summand of Eq. (7) have to be found at the chemical interface by solving corresponding equilibrium and diffusion equations. To find stresses and strains one has to solve the following boundary value problem

V- о = 0, u |Г1 = uo, ° - N 1г2 = to,                  (9)

uk = 0, о - N|r = 0, ch                    ch where conditions across chemical reaction front Γ represent displacements and traction continuity, Γ and Γ are the surfaces where Dirichlet and Neumann boundary conditions are applied respectively.

It is assumed that the diffusion process is faster than the chemical reaction, so the stationary diffusion equation is considered. Boundary value problem for the diffusion is then

A c = 0,

D N -V c + n to= 0 at the chemical interface,   (10)

DN - Vc + a(c-c*) = 0 at the reactant supply surface, where D is the diffusion coefficient, α is the surface mass transfer coefficient, c is the diffusive constituent solubility, and stress-dependent reaction rate ω appears in the condition on the moving boundary.

If stresses and concentration at the reaction front are found, then the front velocity can be calculated from kinetic equation (2), (3).

2.    Numerical simulations of the reaction front propagation

A hollow cylinder with internal and outer radii, a and b , undergoing the chemical reaction localized on a cylindrical reaction front of the radius ρ, is considered. The problem is stated in a plain strain formulation. The reaction starts at the outer surface of the cylinder and propagates towards its center (Fig 2).

Fig. 2. A quarter of the hollow cylinder’s cross-section

For choosing the values of material parameters, we follow the stability predictions obtained from the considerations of stress induced phase transformations [56, 57]. It has been shown that in the case of a two-phase solid sphere with a spherical interface the inequality ц_ >  ц₊ , where ц_ and ц₊ are the shear modules of the inner and outer phases, respectively, is to be fulfilled as a necessary condition of the interface stability. The similar result for a two-phase solid cylinder also follows from the stability analysis presented in [59].

In [57] it was also mentioned that a hole in the center of a sphere could act as a destabilizing factor. In a present work, in order to avoid the influence of the hole on the stability of the reaction front, the inner radius of the cylinder was chosen much smaller than the external radius and the predicted radius of the reaction front at the blocking state, reaction blocking radius, so that a = p <  b . The elastic modules taken for numerical simulations for hypothetically stable and unstable cases are listed in Table 1.

Shear and bulk modules of solid constituents used in numerical simulations

Table 1

	X+ , GPa	X_ , GPa	Ц+ , GPa	Ц _ , GPa
Stable set	25	45	10	25
Unstable set	45	25	25	10

The external load U is adjusted such that the reaction blocking occurs at the same radius for stable and unstable configuration. The initial position of the reaction front is taken far from the blocking position.

Numerical analysis was performed using finite elements. Propagation of the chemical reaction front was realized by remeshing the cylinder. Results of the numerical simulations are shown in Fig. 3, where radial coordinates of twenty equally distributed points along the chemical reaction front are plotted as a function of time.

Fig. 3. Radial coordinates of twenty equally distributed points at the reaction front during the chemical reaction. Light and dark gray curves correspond to stable and unstable configurations espectively

For the stable configuration, light grey curves, all the points of the interface have equal radial coordinates, so the interface remains circular and smoothly converges to the analytically predicted reaction front blocking position. In the case of the unstable configuration, the initially circular interface keeps its shape until it approaches the equilibrium position when the loss of stability occurs. This can be concluded from the deviation of the interface points radiuses: dark grey curves in Fig 3.

The shape of the unstable reaction front in the case of initially circular interface is governed by the accuracy of the numerical calculations. An example of such a shape is shown in Fig. 4.

Fig. 4. Shape of the initially circular interface after stability loss in the vicinity of the reaction blocking

In order to demonstrate the instability kinetics, one can introduce an initial predefined perturbation of the reaction front with given shape. Fig. 5 shows the kinetics of the reaction front with introduced sinusoidal-type perturbations with small amplitude and different frequencies for the unstable configuration. The initial position and the developed shape of the reaction front are shown in solid light grey and black curves, respectively. Dashed curves denote intermediate states and the dotted curve states, for the reaction blocking, the position of the interface.

Fig. 5. Developing of the initially predefined sin-type perturbations in the case of unstable configuration. Initial amplitudes of the perturbations in both cases are equal but the frequency in the left configuration is smaller than in the right one

In the case of a stable configuration, the initially introduced perturbations vanish and the reaction front stops at the blocking state. However, in the case of an unstable configuration, the amplitude of the perturbation amplifies while approaching the blocking state. One should note that during propagation the predefined mode of the perturbations remains the same. The growth of the amplitude forms stress concentration areas, which are analyzed in detail in the next section.

3.    Stresses caused by the loss of the interface stability

We consider a static configuration, where the chemical interface has a predefined sinusoidal-type shape near the blocking state, as shown in Fig 6. The unstable set of parameters is used for the simulations.

Various amplitudes and frequencies for the interface shape were analyzed. Stress distributions along the radiuses

Fig. 7. Distribution of the von Mises stresses along corresponding directions

Fig. 8. Distribution of the hoop stresses along corresponding directions

Fig. 6. Unstable configuration of a solid cylinder with defined shape of the interface. The two directions along which the stresses are analyzed are shown with arrows

Fig. 9. Distribution of the radial stresses along corresponding directions

The maximum von Mises stress in OB direction increases with the amplitude and frequency of the interface shape. As shown in the previous section, the amplitude grows with the propagation of the front in the unstable configuration. This may lead to intensive plastic deformations.

One should note that negative hoop stresses (Fig. 8) in both OA and OB directions do not change their sign in the vicinity of the stress concentration due to amplitude growth of the reaction front perturbation. However, an increasing magnitude of the radial stress in the OB direction (Fig. 9) may lead to delamination. In this case not only the stress distribution but also the conditions for the chemical reaction will change.

Conclusions

An approach to study the influence of stresses on the chemical reaction front propagation was developed based on the chemical affinity concept and implemented in a numerical procedure. Stable and unstable chemical reaction front propagations were modeled numerically. Stress state caused by the growth of the instability amplitude was analyzed. It was shown that negative hoop stresses do not change their sign with a growing amplitude of the reaction front instability. Hence, they cannot cause fracture of the cylinder. However, the von Mises stress increases with the perturbation amplitude growth, which may lead to the plasticity and failure.

Acknowledgement

The authors acknowledge the support of the Russian Science Foundation (grant № 19-19-00552). Authors express special gratitude to Dr. Maxim Frolov for his valuable advices regarding numerical stability and convergence.