Micromechanical modelling of mechanochemical processes in heterogeneous materials

There is a range of practical problems where advanced engineering heterogeneous materials undergo chemical transformations. The primary example of such system is energy storage materials, in particular anodes of Li-ion batteries containing active Si particles. The exploitation of such anodes involves extreme volumetric expansion of the active particles during the chemical reaction. The expansion is causing mechanical stress, which, in turn, influences the kinetics of chemical reactions even up to their arrest. A particular reaction between Si and Li is localised, as well as a number of other reactions, such as oxidation or precipitate formation. The model presented in this paper accounts for the kinetics of the reactions in a collection of particles inside a matrix material. The microstructure is modelled using the multiscale mean-field (MF) framework based on the incremental Mori–Tanaka (IMT) method. This is the first application of a multiscale MF technique to modelling reaction front kinetics in particles and linking the intra-particle kinetics with the response of the matrix. A number of physical effects arising from the influence of the deformation mechanisms of the matrix on the kinetics of the intra-particle reactions is investigated. Furthermore, the applicability of the proposed model and the IMT homogenisation scheme is studied by comparison to the full-field simulations in the cases of small and finite strains.


Introduction
In recent years, there has been an emergence of interest in the influence of mechanical stresses on kinetics of chemical reactions taking place at the microscopic scale, which is relevant, for example, for formation of oxide layers on parts of micro-electro-mechanical systems or other materials [1][2][3], for lithiation of Si in Li-ion batteries [4,5] or lead soap aggregate formation in old paintings [6]. Chemical reactions in these specific examples are localised, i.e. they take place at a surface inside a deformable solid. Furthermore, the chemically transformed material has intrinsic volumetric expansion compared to the untransformed material, which causes mechanical stresses. In the case of Si oxidation, crack initiation might take place in the oxide layer, which leads to fatigue failure [1]. When lithiation of Si takes place in Li-ion batteries, there is reaction arrest due to a build-up of mechanical stresses [7].
There are several models of stress-affected chemical reactions [8][9][10][11][12][13][14][15][16]. In the localised reactions, the chemical reaction front moves due to the consumption of the diffusive reactant, which is supplied to the reaction front by a diffusion process. The chemical reaction rate and, thus, the velocity of the front is affected by the stresses, which, in turn, emerge due to transformation of the material as the front moves. There is a similarity between the mechanochemical transformation and the movement of phase boundaries in the classical stress-induced solid-solid phase transformations, where the configurational force, which is driving the phase boundary kinetics, is equal to the jump of the normal component of the Eshelby stress tensor at the phase boundary [17][18][19]. In [12][13][14], it has been shown that the configurational force for a localised chemical reaction front is the normal component of a specific tensor, which has been called the chemical affinity tensor. This tensor is equal to the combination of the chemical potential tensors of the reaction constituents, which, in turn, are equal to the Eshelby stress tensors divided by the reference mass densities.
Theory built on the chemical affinity tensor has been applied to a number of problems with linear elastic constitutive behaviour of materials: an analysis of a planar reaction front [13], a spherical reaction front [20] and a reaction front propagation in a 2D plate with a groove [21]. Furthermore, the theory has been used to simulate the so-called two-phase lithiation of Si particles utilised in Li-ion batteries [22], where the constituent materials undergoing large elasto-viscoplastic deformations were considered. The work has been focused on single particles with stress-free boundary conditions and on the analysis of the influence of various model parameters on the reaction locking effect. When Si particles are used as anode active material in Li-ion batteries, they are placed inside a matrix material, which usually consists of a graphite and a polymer binder, and analysis of the influence of the matrix material on the kinetics of the two-phase lithiation, has not been performed yet.
There are multiple ways of modelling particle-reinforced heterogeneous materials via micromechanical full-field (FF) and mean-field (MF) approaches. For a recent review of the FF homogenisation techniques, the reader is referred to [23], while a recent overview of the MF techniques can be found in [24][25][26]. FF homogenisation techniques have an advantage over the MF techniques in terms of accuracy; however, they require the FF solution of the underlying PDEs. In relatively standard problems, e.g. first-order computational homogenisation in application to mechanics, which is now well-established, this is not an issue, especially considering available computational tools. However, the microstructure considered in this paper undergoes chemical transformation with moving phase boundaries. Such problems are non-trivial themselves and require special techniques for handling the moving phase boundary, e.g. XFEM [27,28], the enhanced gradient FEM [29,30] or CutFEM [31] approaches. In the case of absence of the chemical reaction front, the FF homogenisation has already been formulated for problems of stress-affected diffusion of Li-ions and applied to some test cases in [32][33][34].
The MF techniques rely on a number of assumptions, which simplify the underlying problem. The one of the established MF homogenisation techniques is the Mori-Tanaka inclusion-matrix coupling scheme [35,36]. This scheme has already been used in the context of coupled stress-affected electrochemical modelling of Li-ion batteries [37][38][39], however, only for volumetric reactions and linear elastic materials. Furthermore, this scheme has also been applied to obtain the effective properties of the composite electrodes [40,41], again, in the elastic regime. The Mori-Tanaka scheme has originally been proposed for linear elasticity [35,36], later extended to non-elastic behaviour using an incremental formulation [42] and, finally, has also been formulated for finite deformations [24,25]. However, up to now, the Mori-Tanaka approach has not been applied to heterogeneous materials, where mechanochemical transformation takes place and involves propagating reaction fronts, in both smallstrain and nonlinear finite-strain settings.
Therefore, the aim of this paper is to model the mechanochemical behaviour of particles with localised stress-affected chemical reactions, embedded into a matrix material, using a micromechanical approach based on the Mori-Tanaka model generalised to a finite-deformation regime. The finite-deformation problem formulation is used, as it also allows studying both the case of small and finite strains in a natural way. Furthermore, a validation procedure, which is based on a limit case of small particle volume fraction, is carried out to estimate the accuracy of the incremental Mori-Tanaka (IMT) scheme in the cases of small and finite strains.
It is shown that, the IMT scheme can be considered as useful approach for micromechanical modelling of heterogeneous materials, undergoing mechanochemical transformations. However, this work shows that extra care must be considered generalising the model to finite strains. This is due to inaccurate estimation of the local fields, which highly affect the reaction kinetics, in the finite-strain regime. On the other hand, the validation showed perfect consistency of the model in the case of small strains. In the results section of this paper, the dependency of the intra-particle reaction front kinetics on the matrix material parameters (stiffness, relaxation time, etc) is studied. transformation, the particles become two-phase structures, consisting of the shell (the transformed material) and the core (the untransformed material), which are separated by the chemical reaction front. As the transformation advances, the reaction front moves from the external surface of the particle to its centre. The transformed material has intrinsic volumetric expansion, which causes inhomogeneous stress and strain spatial distribution within the particle and the volumetric expansion of the particle itself.
The model of the composite material is built upon the IMT homogenisation scheme [42]. This scheme was originally designed for a two-phase composite with homogeneous inclusions incorporated into a homogeneous matrix. Since in the considered case, the inclusions have a two-phase structure, the IMT scheme is used to link the matrix strain and the effective inclusion strain, while the intra-inclusion structure and stress-strain state are resolved by a separate sub-model, as schematically shown in figure 1. Thus, the description of the micromechanical model is split into three main parts: the inclusion-matrix coupling, the intrainclusion sub-model and the constitutive laws of the materials involved in the model.
To achieve computational efficiency, the micromechanical model uses the assumption that the inclusion geometry remains spherical throughout the transformation and the deformation processes. This allows using the spherical symmetry and solving 1D equations for the intra-inclusion quantities. Although the underline equations become 1D, the model is still presented in tensorial form.

Inclusion-matrix coupling and volume-averaging
For completeness of the paper, the ideas behind the IMT scheme are briefly summarised. The IMT scheme is a modification of the classical Mori-Tanaka (CMT) MF homogenisation scheme [35,36], which has been designed for the particle-reinforced composites with linear elastic constituents. The CMT scheme is a way of coupling inclusion and matrix stress-strain states and linking these states to the macroscopic (volume-averaged) quantities. It must be emphasised that the CMT scheme operates specifically with the volume-averaged matrix stress/strain and does not resolve the local details. The CMT scheme is built using three assumptions: • the inclusion strain is assumed to be spatially homogeneous (the matrix strain is spatially homogeneous by definition, as it is the volume average), • the inclusion strain is assumed to be given by the solution of the Eshelby's problem, • in which the strain applied at the infinity (the far-field strain) is assumed to be equal to the matrix strain in the composite. The IMT scheme, in turn, generalises the CMT scheme for the case of arbitrary constitutive behaviour of the phases. It is founded on the assumption that at any moment of time, the linearised with respect to time inclusion and matrix strains are coupled by the CMT scheme.
In the IMT scheme, the tangent stiffness tensors are used in the Eshelby's solution. This implies the following important assumption of the IMT scheme: • the tangent stiffness of the matrix is assumed to be spatially homogeneous.
Although the scheme is called 'incremental', it is often written in a rate form, e.g. [43,44], which is presented below. When the rate form is discretised in time, the incremental form is recovered. Although the IMT scheme in application to finite deformations and nonlinear material behaviour contains a number of relatively restrictive assumptions, it is wellestablished and a satisfactory quantitative match of this scheme with the FF simulations has been observed in a number of cases, see literature reviews in [24,25].
2.1.1. Spherical inclusions-effective quantities. As mentioned above, during mechanochemical transformation, the intra-inclusion stress and strain are spatially inhomogeneous. Determination of these quantities is discussed in section 2.2. The IMT coupling scheme, however, operates with homogeneous quantities. Therefore, the effective inclusion Cauchy stress and deformation gradient 7 are introduced as respectively, where I is the second-order unit tensor, T is the traction at the edge of the particle, r a is the current radius of the particle, R 0 is the initial radius of the particle and J a is the volume change of the particle (details are given in section 2.2.1).
where f and f 0 are the current and the initial (reference) volume fractions of the inclusion, respectively, tensors s m and F m are the average Cauchy stress and the average deformation gradient of the matrix, respectively. Such averaging is common in micromechanics, e.g. [45]. The current and reference volume fractions are related as and D m are the macroscopic, the inclusion and the matrix rate of deformation tensors, respectively. The derivation of the equations is given in appendix A. It must be emphasised that equation (7) is only valid for the case when the inclusion is spherical and the macroscopic deformations are exclusively volumetric.
2.1.4. Inclusion-matrix coupling. As discussed above, in the original IMT scheme (formulated for small strains), the inclusion strain is coupled to the matrix strain using the solution of the Eshelby's problem. Furthermore, in the original IMT scheme, the inclusion and the matrix strains are usually written as functions of either the incremental macroscopic strain or the incremental macroscopic stress. Here, first, the rate formulation of the IMT scheme is used, second, the finite-deformation framework is used. There can be two distinct cases-the macroscopic deformation-controlled case and the macroscopic stresscontrolled case.
When F ⟨ ⟩ is imposed (the deformation-controlled case), the following relations 8 are used: : , , 4 a 4 a 4 m 4 m are fourth-order stress/strain concentration tensors, which are functions of the tangent stiffness tensors of the inclusion and the matrix, current stresses of the inclusion and the matrix, current volume ratios of the inclusion and the matrix, the Eshelby tensor of the spherical inclusion in the matrix with averaged properties. The expressions for these tensors are derived in appendices B and C.
As the particle expands during the transformation of the shell, the effective stress-free rate of deformation of the particle D* is introduced in the scheme. It is defined as the effective inclusion rate of deformation under zero traction rate, The inclusion is modelled using the mechanochemical model published previously [22]. This model resolves the kinetics of the chemical reaction front using the chemical affinity tensor [12][13][14], which acts as the configurational driving force for the front and was derived from the balance laws and the dissipation inequality for materials with arbitrary rheology.
2.2.1. Mechanics. The model equations are given in a tensorial form, although a specific form of the deformation gradient is used in the model, due to the spherical symmetry restriction of the particle. More specifically, the current position of a material point inside the particle is denoted as r R ( ), which is a function of the reference position of the material point R. Thus, the intra-particle deformation gradient is .
The intra-inclusion Cauchy stress, in turn, is denoted as s R i ( ). Furthermore, as the particle consists of the core and the shell, the quantities within the core and the shell are denoted using subscripts '−' and '+', respectively, where the position of the reaction front mapped onto the reference configuration of the untransformed material is denoted as Γ. It is useful to denote the traction at the edge of the particle and the current total radius of the particle as s = T R e e : , 15 The intra-inclusion stress and deformation are found by solving the equilibrium equation with either imposed traction boundary conditions or imposed displacement boundary conditions. The former boundary condition is used to calculate the effective stress-free expansion of the particle (details are given in section 3.3.1), while the latter boundary condition is used to calculate the current stress distribution within the particle coupled to the matrix. Equation (17) is solved with respect to the unknown function r R ( ). Furthermore, as the shell undergoes the transformation strain, the displacement and traction continuity conditions are enforced at = G R . It must also be noted that operator ∇ is defined with respect to the current configuration.
2.2.2. Diffusion. The reaction front moves due to the consumption of the reactant, which diffuses through the transformed material from the outer surface of the particle to the reaction front. The diffusion flux is projected onto the reference configuration of the untransformed material and the quasi-stationary diffusion equation is considered, where n * andn are the stoichiometric coefficients of the reactant and the untransformed material, respectively, in the chemical reaction; ris the density of the untransformed material in the reference (undeformed) state and -M is the molar mass of the untransformed material; V is the velocity of the reaction front and is discussed in section 2.2.3; D is the diffusion coefficient of the reactant through the transformed material; α is the surface mass transfer coefficient; c * is the solubility of the reactant in the transformed material. Operators ∇ 0 and Δ 0 are defined with respect to the reference configuration of the untransformed material.
2.2.3. Chemical reaction. The chemical reaction results in the propagation of the chemical reaction front from the outer surface of the particle towards its centre. The expression for the normal component of the reaction front velocity with respect to the reference configuration of the untransformed material is given by where k * is the kinetic coefficient, R g is the universal gas constant, θ is the temperature and A NN is the normal component of the chemical affinity tensor [12][13][14]. If the pressure produced at the interface by the diffusive constituent is neglected (as it is small compared to stresses produced by the transformation strain), then this component is given by where γ is the temperature-dependent combination of the chemical energies of the constituents, g 3 is the volumetric expansion ratio of the material due to the chemical transformation, -P is the first Piola-Kirchhoff stress of the untransformed material, the double square brackets denote the jump of the quantity at the reaction front. Furthermore, -W and + W are the elastic strain energy densities per unit volume of the untransformed and the transformed materials, respectively, with respect to their reference configurations. It must be emphasised that + W is represented in the reference configuration of the transformed material, since the constitutive laws for the transformed material are defined in its reference configuration.

Mechanical constitutive laws
As mentioned above, the micromechanical model uses an assumption that the inclusion remains spherical throughout the transformation and the deformation processes. Therefore, the effective inclusion stress and deformation gradient, s a and F a , respectively, are exclusively volumetric. Hence, if the macroscopic stress and deformation gradient, s ⟨ ⟩ and F ⟨ ⟩, respectively, are also volumetric, then the matrix stress and deformation gradient, s m and F m , respectively, will also be volumetric, according to (3) and (4). Therefore, to simplify the model, the micromechanical model is limited to volumetric s ⟨ ⟩ and F ⟨ ⟩ and, thus, the nonvolumetric components of s m and F m are excluded. It must be emphasised here that in the Mori-Tanaka scheme, the matrix stress and strain are volume-averaged quantities. Therefore, locally, the matrix can deform deviatorically, while the volume-averaged deviatoric stress is zero. This is implicitly used in the scheme when the Eshelby's solution is adopted for the particle.
where -K is the bulk modulus and -J is the volume change of the material of the core. For completeness of the equations, the expression for -P used in equation (22) is given in appendix D.
2.3.2. Particle shell. The shell of the particle consists of the transformed material, which is modelled as the finite-strain nonlinear viscoelastic material with hardening. The total deformation gradient is decomposed into elastic, viscoplastic and chemical parts, where F g is the transformation strain (analogous to thermal expansion of the material), which is taken to be isotropic. The following elastic strain energy density is used: where + K is the bulk modulus, G h is the hardening modulus and )is the total shear modulus of the material. It must be emphasised that the elastic strain energy density is given in the reference configuration of the transformed material. Here B ē and B ā denote the elastic and the total (with respect to the reference configuration of the transformed material) isochoric Finger tensors (definitions are given in appendix D), respectively. This corresponds to the following Cauchy stress: where superscript 'd' denotes the deviatoric part of a tensor. Finally, the evolution of the plastic deformation is prescribed using the plastic rate of deformation tensor (definition is given in appendix D): ( ) respectively, where λ e and λ p are the elastic and plastic stretches, which correspond to volumetric deformations. Analogous to the intra-particle constitutive laws, the following Cauchy stress is used: where K h is the volumetric hardening modulus. The total bulk modulus of the material is The model is completed with the following law of the evolution of the plastic deformation: where, analogous to the intra-particle constitutive laws, t m and q m are material parameters. Parameter σ 0 is just a normalisation constant and is chosen to be the same as in section 2.3.2. ( ) The evaluation of the derivatives at a constant deformation rate results from the consideration of the strain-rate dependent materials. This is important because the tangent can be understood as the stiffness of the material in response to small-strain perturbations at the current (deformed) state, and, since the stress-strain response of the material changes with the strain rate, the tangent is also strain-rate dependent. Using K a and K m , the following tangent compliance tensors for the inclusion and the matrix, respectively, are obtained: being the basis vectors. Here, G a and G m are the tangent shear moduli of the inclusion and the matrix, respectively. As mentioned above, the inclusion is allowed to deform only volumetrically, therefore, its effective shear modulus is infinite. Thus, when concentration tensors G H Q R , , , 4 a 4 a 4 m 4 m are calculated, limit  ¥ G a is taken. Furthermore, since the matrix material is modelled as volumetrically nonlinear viscoelastic material, the deviatoric state of the matrix material is assumed to be elastic, therefore, G m is a material parameter.

Numerical solution procedure
The micromechanical model relies on the resolution of the intra-particle stress fields and the reaction front kinetics. This part is described in detail in section 3.6 of [22], where the intraparticle kinetics of an isolated stress-free particle has been considered. Here, only the general idea of the intra-particle solution scheme is outlined. Afterwards, the global computational scheme (i.e. the micro-macro coupling) is discussed. The scheme is non-trivial and the motivation for such choice is explained first, with the discretised set of equations provided afterwards. In this section, only macroscopic strain-controlled case is considered, i.e. equations (8) and (9), since the macroscopic stress-controlled case is identical from the implementational point of view.

Intra-particle quantities
Equation (18) is solved analytically, since the problem is spherically symmetric, and the closed form of concentration c as a function of velocity V is obtained. Next, equation (17) is rewritten in the reference configuration of the untransformed material and is discretised using the finite-element method. Since the problem is spherically symmetric, 1D mesh with the equidistant nodal spacing is used; the finite-element mesh is fixed in the time domain. Finally, the chemical reaction front can be located only at the nodal positions and the reaction front movement is resolved by the reaction front moving from one nodal position to a neighbouring node. Since the extent of the movement is known, which is the element size, the time required for the move, which is the time step, becomes unknown and is obtained by solving the discretised version of equation (21). Thus, within a time step, the unknowns of the problem are the nodal solution and the time step. The discretisation in time of the equation for the velocity and the constitutive relations is implicit. The equations are solved using the Newton-Raphson method.

Necessity for implicit/explicit hybrid scheme
The micro-macro coupling equations (8) and (9) are given in the rate form and there are several ways of obtaining their discrete form. The straightforward approach would be the implicit discretisation in time. However, this would lead to an inefficient scheme due to one specific interdependency, which is explained in this subsection.
The micro-macro coupling equations contain the concentration tensors, which depend on K a and K m , and the stress-free rate of deformation D*. Since it is not possible to calculate K a analytically for non-trivial constitutive laws, K a must be calculated numerically by small perturbation of the boundary conditions of equation (17) and by subsequent solution of the intra-particle model. Similarly, D* cannot be calculated analytically either and must be calculated numerically by solving the intra-particle model with the moving front and zero applied traction rate, equation (10). Therefore, at each moment in time, at least three different intra-particle stress states must be obtained by solving the intra-particle model: the stress state with the actual boundary conditions coming from the micro-macro coupling equations, the stress state with the perturbed boundary conditions to calculate K a and the stress state with the zero traction rate boundary conditions to calculate D*.
Equations (8) and (9) can be discretised implicitly or explicitly, which corresponds to tensors G H Q R , , , 4 a 4 a 4 m 4 m being taken at the current or at the previous increments, respectively. If the implicit formulation is used, the above-mentioned problems of obtaining three different intra-particle stress states become interdependent via the boundary conditions and the time steps. The solution of the problems, in this case, involves operations with large matrices and a significantly more complex code compared to the case when the problems are solved separately. To be able to solve the problems sequentially, equations (8) and (9) can be discretised explicitly.
During the development of the numerical implementation of the micromechanical model, at first, an explicit representation of the micro-macro coupling equations has been implemented. This, however, resulted in a numerically unstable model. This instability of the explicit scheme was similar to the phenomenon encountered in stiff systems, where the time step must be extremely small for the numerical stability. Here, however, as explained in section 3.1, the time step is varying with time when the intra-particle model is resolved, since the time step is connected to the reaction front 'jumping' from one node to another on the mesh with equidistantly-spaced nodes. Thus, when the velocity drops, the time step increases, which initiates the instability in the fully explicit scheme. The problem of instability has been resolved by employing the implicit/explicit hybrid discretisation scheme of the micro-macro coupling equations, which is explained below.

Computational scheme
The scheme starts by assuming that all quantities at the previous time stept j 1 are known. The aim of the scheme is to find the incremental changes of the inclusion and the matrix strains and to find the intra-inclusion stress and strain distributions at the current time step, subject to the transformation strain of the inclusion, i.e. volumetric change, and the imposed macroscopic stress or deformation.
3.3.1. Intra-particle quantities and stress-free expansion. From the previous time increment t j 1 to the current time increment t j , the reaction front inside the particle moves by ΔR, which is the finite-element length in the intra-particle model. The positions of the reaction front at t j 1 and t j can be defined as Gj 1 and Γ j , respectively. The time step Δt that is required for this movement results from the solution of the discretised version of the equation for the velocity (21), as explained in section 3.1. It is important to note that since the velocity is stressdependent in this model, the time step also becomes stress-dependent and, thus, also dependent on the traction applied to the particle.
The intra-particle quantities depend on the position of the reaction front. Furthermore, when traction is applied to the particle, the intra-particle quantities also become dependent on the traction. Therefore, the current radius of the particle can be represented as the function of the front position and the traction, ). On the other hand, when the displacement boundary conditions are applied to the particle, it is convenient to denote the traction at the edge of the particle as the function of the front position and the imposed particle radius, G T r , a ( ). Using the above definitions, the incremental stress-free expansion of the particle within a time step can be written as where -T j 1 denotes the traction at the previous time increment. In equation (40) the reaction front moves from Gj 1 to Γ j , which causes the expansion of the particle, while the traction is kept constant, i.e. the traction time derivative is zero. For convenience, the following short notation is introduced: (40) is the solution at the previous time step, i.e. already known, while the determination of r a * involves solving the intra-particle model with the imposed traction boundary conditions for equation (17). The time step, with which solution r a * is calculated, is denoted as Dt*. The physical meaning of this time step is the time required for the reaction front to move by one element length while the traction applied to the particle stays constant. This is important as, obviously, Thus, the micro-macro coupling equations are discretised using an implicit/explicit hybrid scheme, where bulk moduli K a and K m , are taken at the intermediate time increment, while all other functions are taken at the previous time increment. Such choice avoids solving the intraparticle problem for finding K m as a root of a nonlinear equation. As mentioned above, the strain increments in (41) and (42) are incremental strains taking place fromt j 1 to t j *. However, to solve the intra-inclusion problem, incremental strains from where e D a denotes the inclusion incremental strain taking place fromt j 1 to t j . Here, derivative e ¶ ¶D D t * * is found computationally, as explained in section 3.3.3. Using e D a , the boundary condition for the intra-inclusion problem can be formulated: Time step Δt, which enters the boundary condition, is unknown and is solved for in the Newton-Raphson loop for the intra-inclusion solution. Such choice allows decoupling the calculation of the current intra-inclusion stress-strain state and the stress-free incremental expansion of the inclusion. The computational scheme is completed by formulation of the discrete equations for K a and K m . Equation for K m can be obtained analytically and is provided in appendix E. Due to the intra-particle model given by the set of nonlinear PDEs, K a can only be obtained numerically.
To obtain it, the perturbation of the traction at t j * is defined as where dr a is a small perturbation of the radius and r a * is already known, as it has been used for calculating e D * . The first term in expression for dT is obtained by solving the intra-particle model with the displacement boundary conditions for equation (17) and with time step d D + t t * instead of solving the equation for the velocity, which is explained below. This leads to the discretisation of equation (36) for K a at t j *, The way how the time step is calculated in the intra-particle model is modified for the calculation of the first term in expression for δT and it is required for enforcing the constant deformation rate when K a is evaluated. To enforce = r r const a ȧ according to (36), time δt it takes for the particle to expand from r a * to d + r r a a * should be found from The derivative of r a * by Dt* is calculated numerically, which involves solving the intraparticle model with the perturbed time step. Furthermore, this derivative is related to the derivative of the incremental expansion with respect to the time step in (43) 2. Time t j *. (a) Find the incremental stress-free expansion of the particle e D * by solving the intraparticle model under previous time step traction -T j 1 , equation (40). (b) Find the derivative of e D * by Dt*, equation (47), by solving the intra-particle model with perturbed time step. (c) Find the effective bulk modulus of the particle K a , equation (46), by solving the intraparticle model with perturbed radius. (d) Solve nonlinear equations for the tangent bulk modulus of the matrix K t j m * ( ) and the matrix incremental strain e D m , equation (42). (e) Calculate the incremental expansion of the particle e D a at t j *, equation (41).

Time t j .
(a) Extrapolate e D a at t j * to e D a at t j , equation (43). (b) Apply displacement boundary conditions to the particle based on e D a , equation (44), and find the current intra-particle stresses and strains by solving the intra-particle model. (c) Update volume fraction of the particle according to equation (5).

Results
Micromechanical modelling of heterogeneous materials with particles undergoing chemical transformation allows investigating the influence of the material composition (e.g. particle volume fraction) and the properties of the matrix material on the kinetics of the localised chemical reaction taking place within the particles. This section aims at investigating these effects and is structured as follows. First, the reaction front kinetics is discussed for a free particle (i.e. not inside the matrix), section 4.1. Next, the IMT model is applied and validated for small strains by performing a comparison with the FF FE simulations, section 4.2. Afterwards, effects of the matrix properties on the kinetics of the chemical reaction front are studied for the case of small strains, sections 4.3 and 4.4. This is followed by a finite-strain example of lithiation of Si particles inside anodes of the next-generation lithium-ion batteries, section 4.5. Next, an effect matrix relaxation on kinetics of the reaction is also studied, section 4.6. Finally, a comparison to the FF simulations is performed for the finite-strain setting in section 4.7.
It is well known that any type of MF modelling has a discrepancy with the FF models due to restrictive underlying assumptions. Usually, the comparison is performed at the macroscopic level by comparing the stress-strain dependencies, e.g. [25], and the discrepancy is usually quantitatively acceptable. On the other hand, when the local quantities are compared, e.g. inclusion-averaged stresses, significant discrepancies can be observed. In the case of the proposed micromechanical model, the local quantities are essential, as they determine the intra-particle lithiation kinetics. Therefore, it is useful to compare the proposed MF and the FF models, before using the MF model to predict various physical effects.

Velocity profile
Since the reaction front in the particle is assumed to remain spherical, the velocity is given by the radial component V which is the same for all points of the front. As the reaction front moves from the edge of the particle towards its centre, the velocity changes due to the buildup of mechanical stresses. Here, as in the previous study [22], the materials of the particle core and the particle shell are taken to be nonlinear elastic and elasto-viscoplastic, respectively. Typical profiles of the velocity V in the particle without the surrounding matrix material are presented in figure 2. The exact shape of the curve depends on the transformed/ untransformed material parameters. In general, the profile has three distinct regions: the increase of the velocity, approximately constant velocity and rapid decrease of the velocity to zero corresponding to the reaction arrest. Depending on the parameters of the problem, these regions can either more or less pronounced, as for example can be seen by comparing figures 2(a), (b) and 6(a). The latter figure is discussed in section 4.5.2. Such complex behaviour results from the interplay of the components of the chemical affinity tensor and is discussed in the previous publication by the authors [22] and earlier in [13,20].
For stress-free spherical particle, given that the transformed phase behaves as elastoplastic material, the velocity profile qualitatively follows the behaviour of the radial stress at the reaction front. Initially, the radial stress at the front is tensile and increases, which accelerates the chemical reaction, as in the case of linear elastic spherical particle [20], where it has been shown that the tensile radial stress in the shell of the particle with stress-free boundary conditions leads to the increase of the reaction front velocity. In the case of an elastoplastic shell, at some point, the increase of the stress slows down, since the formation of the new material requires stretching already expanded and plastically-deformed shell, which superposes compressive stress onto the radial stress at the front. This leads to the decrease of the stress to negative values-the compressive regime, which slows down the reaction. Although stresses at the reaction front enter only one of the terms of the normal component of the chemical affinity tensor A NN , equation (22), such increase-decrease pattern is further magnified by the exponent of A NN within the equation for the velocity, equation (21).

Model validation for small strains
As mentioned in the introduction, the FF 3D finite-element models of the considered composite material with particles undergoing a localised chemical reaction are relatively difficult to construct and solve due to the non-trivial problem of the resolution of the reaction front kinetics. To avoid constructing full 3D models, a limited comparison can still be performed in a simplified spherically-symmetric setting. When the volume fraction of the particle in the MF model approaches zero, the model represents an isolated particle in the infinite matrix (which, in the case of linear elasticity, is the classical Eshelby's inclusion problem). The Figure 2. The dependence of the reaction front velocity on the normalised reaction front position measured from the edge of the particle. Comparison of an isolated particle, the particle inside the matrix according to the MF model and the particle inside the matrix according to the FF model. Subfigures (a) and (b) use different parameter sets, which can be found in tables F3 and F4, respectively. corresponding problem in the FF setting is the particle encapsulated into a very thick 'coating' and with the thickness of the 'coating' approaching infinity. Although such comparison is a limiting case of a volume fraction approaching zero, it is still useful and shows whether the MF and FF models match at least for small volume fractions (dilute regime).
The focus of this paper is on local kinetics of the chemical reaction, which is affected by local quantities, such as tractions applied to the particles. Therefore, the comparison which is performed here is non-standard, since usually only the macroscopic quantities are compared. The reaction front velocity has been calculated for three different cases: an isolated particle (i.e. without the surrounding matrix material), the particle inside the matrix using the MF model and the particle inside the matrix using the FF model. The first calculation is necessary to show that the presence of the matrix has an effect on kinetics, while the other two provide the comparison between the models. 4.2.1. Linear elastic matrix. When the matrix material is linear elastic and the particle volume fraction tends to zero, the exact match between the MF and the FF models is expected because both of them become the classical Eshelby's inclusion problem. This is expected irrespectively of the constitutive behaviour of the particle, as the Mori-Tanaka scheme takes the entire transformation strain of the particle for the calculation of the stresses and strains in the matrix. In figure 2(a), the particle velocity as a function of reaction front position is plotted for three above-mentioned cases. The important result here is that there is a perfect match between the MF and the FF models, as expected for the linear elastic matrix. Parameters that were used in the simulations can be found in tables F1, F3 and F5. In the case of the MF model, the volume fraction has been taken to be f 0 =0.001. In the case of the FF model, the relation between the radius of the particle with the 'coating' and the particle has been 10:1.

Elasto-viscoplastic matrix.
When the matrix material is not linear elastic, the exact match between the MF and the FF models cannot be expected anymore even for the infinitely small volume fraction, due to assumptions of the MF model. In figure 2(b), the kinetics of the reaction front in the particle placed inside the elasto-viscoplastic matrix is compared for the proposed MF and the FF model. Again, three cases are considered to show that the matrix indeed has an effect on the kinetics, which is studied below in more detail, and it is seen that the discrepancy between the MF and FF is unnoticeable. It must be emphasised that the velocity depends on local quantities, such as tractions applied to the particle, which do not usually match well when the MF-type and FF models are compared. Therefore, the match observed in figure 2(b) is specifically valuable and shows the validity of the proposed model for small strains. Thus, it can be concluded that at least for small volume fractions and small strains the proposed model has been validated and matches well the FF simulations. Parameters that were used in the simulations can be found in tables F1, F4 and F6. The same initial volume fraction in the MF model and the particle/coating ratio as in figure 2(a) have been used.

Influence of matrix stiffness on kinetics
Having validated the model under original Mori-Tanaka assumptions (i.e. low particle volume fraction, small strains), it is now possible to study various effects arising from the interaction of the deformation mechanisms of the matrix and the reaction front kinetics within the small-strain setting. To perform simulations in this subsection, the bulk modulus K 0 m was assumed to range from 250 MPa to 32 GPa in logarithmic scale. The volumetric hardening modulus K h was taken to be 1/5 of the bulk modulus K 0 m for each simulation. Finally, the relaxation time τ m was fitted for each bulk modulus, such that true strain at the yield point was 0.7% at the strain rate of --10 s 2 1 . The volumetric stress-strain behaviour of the considered matrix materials is shown in figure 3(a). The parameters corresponding to the simulations be found in tables F1, F4 and F7. The initial volume fraction of f 0 =0.2 was taken for the simulations.

Influence of the matrix bulk modulus on the front velocity.
When the bulk modulus of the matrix is increased, the velocity decreases, as shown in figure 3(b). This result is obvious, since higher matrix bulk modulus leads to higher constraints from the matrix on the particle, hence higher stresses in the particle, which slow down kinetics. The somewhat more complicated result to interpret is the nonlinear effect of the matrix bulk modulus on the velocity. Initially, up to 1 GPa, the effect of the bulk modulus is small, however, it becomes significant for the intermediate values and diminishes again at large values.
The decay of the influence of the bulk modulus at high values can be explained by the matrix becoming so stiff after a certain threshold value, that the particle almost cannot expand the matrix volumetrically. However, it should be emphasised that the extremely large bulk modulus does not mean that the matrix becomes essentially a fixed displacement boundary condition for the particle. The reason for this is that only an average far-field behaviour of the matrix is represented in the model and the Mori-Tanaka matrix-particle coupling is used, which relies on the Eshelby's solution for the particle. In the Eshelby's problem with the incompressible matrix, the particle can expand due to the shear mechanism of the matrix, while the shear stress at the infinity is zero. Here, the Mori-Tanaka assumptions are essential, since the scheme assumes that the matrix stresses of the Mori-Tanaka composite are the stresses at the infinity in the Eshelby's solution. Thus, in the Mori-Tanaka scheme, the inclusion is allowed to expand into the almost incompressible matrix due to the shear of the matrix, while the average far-field stresses (which are simply called 'the matrix stresses' in the Mori-Tanaka scheme) remain volumetric and small due to almost near incompressibility. figure 3(b) is carefully examined, it can be seen that there is a kink in the velocity profile in the region of the reaction front position between 0 and 0.2. This kink is the result of the onset of the plastic deformation in the matrix, must be considered. The Mori-Tanaka coupling scheme provides the behaviour that is between the Reuss bound ('sequential' connection of materials) and the Voigt bound ('parallel' connection of materials). To understand the observed behaviour of the Mori-Tanaka composite, an auxiliary 1D case can be examined-two linear elastic materials with different stiffnesses and with the first material undergoing transformation strain (e.g. thermal expansion). The second material does not have any transformation strain. Furthermore, materials have equal volume fractions. The materials are coupled using the Reuss and the Voigt schemes and the total stress is zero. When these materials are coupled sequentially, everything is trivial-all stresses are zero and the total strain is equal to the transformation strain. When these materials are coupled in a parallel way, it is easy to see that This means that when stiffness of the second material E 2 is small, the total strain is close to the transformation strain and the stress in the second material is close to zero. For large E 2 , the total strain becomes small and the stress in the second material becomes large (close to E 1 times the transformation strain). Furthermore, the strain of the first material is equal to the total strain. Since the Mori-Tanaka composite is expected to be between the Reuss and the Voigt bounds, the effects of E 2 should be less profound than in the Voigt scheme, but follow it qualitatively in this example with two materials. Exactly these effects are qualitatively observed in figures 4(d) and (a), for the total strain of the composite and the matrix stress, e c and s m , respectively. With the increase of the matrix bulk modulus, the total strain decreases, while the matrix stress increases approaching a certain limit. In the simulations, due to the elasto-viscoplasticity of the matrix and the selected problem parameters, the elastic matrix strain is close to 0.7% (depending on the strain rate). According to the auxiliary 1D example above, the matrix strain should also decrease with the matrix bulk modulus. This means that the matrix plastic strain should also decrease with the increase of the matrix bulk modulus, as seen in figure 4(b).

Stresses and strains in the matrix. When
The decrease of the matrix plastic strain and the increase of the matrix stresses in figures 4(b) and (a), respectively, also match from the quantitative point of view. For example, when the bulk modulus increases from 250 MPa to 8 GPa (32 times increase), the matrix plastic strain decreases from 3.57% to 0.63% (5.7 times decrease) at front position 0.5. From this it can be concluded that the stresses should definitely increase, as the increase is the bulk modulus in significantly higher than the decrease of the strains.

Stress relaxation.
There is another interesting effect that can be observed in figure 4(a). At some point, the stress relaxation takes place and the matrix stress starts decreasing. The reason for this is that matrix is modelled as a viscoelastic material with highly nonlinear stress-dependent viscosity, which effectively behaves as elastoplastic material. However, when simulation times become large, the viscous effects start dominating and the stress relaxation takes place. Figure 3(b) shows that there is a region between front positions of 0.6 and 0.9 when the velocity of the front rapidly slows down and approaches zero. This means that the time required for the front to move by certain distance increases nonlinearly, which eventually leads to stress relaxation. Furthermore, the asymptotic decrease of the velocity to zero is directly linked to the nonlinear viscoelasticity of the transformed material (particle shell). This effect is further discussed in the previous study [22] and compared to the case of linear elastic particle where the velocity does not approach zero asymptotically [20].
The qualitative changes in the behaviour of the matrix stress can be seen more clearly when the average matrix tangent K m is plotted in figure 4(c). Due to nonlinear constitutive law, this tangent depends on the average deformation of the matrix, which, in turn, depends on the deformation of the inclusion, which is governed by the reaction front position. In the initial section of the plot, tangent decreases slightly due to nonlinearity of the material. After this, there is a rapid elastoplastic transition, which is also seen as a kink in figure 4(a). The transition point is followed by the almost linear decrease of the tangent, which now roughly corresponds to the volumetric hardening modulus of the matrix. Starting from front position of 0.6, which coincides with the velocity rapidly decreasing, figure 3(b), tangent again changes the slope and goes below zero, corresponding to the stress decrease. Finally, there is one more change of slope in the tangent around front position of 0.8, which coincides with the change of sign of the second derivative of the front velocity, i.e. when the velocity is already small (below -1 nm s 1 ) and it starts approaching zero asymptotically. This last kink in the tangent modulus also corresponds to the kinks around 0.8 in figures 4(a) and (b) for the matrix stress and the matrix plastic strain, respectively, which is especially visible for the large bulk modulus. As a final point, the obvious fact can be noted that the starting point of the matrix tangent in figure 4(c) corresponds to the matrix bulk modulus.

Influence of particle volume fraction on kinetics
The reaction front velocity V increases with the increase of the initial volume fraction f 0 , which has been varied from 0.05 to 0.3 linearly, as shown in figure 5(b). The parameters corresponding to the simulations be found in tables F1, F4 and F8. To explain the behaviour of the velocity of the reaction front as a function of the particle volume fraction, an auxiliary 1D example from section 4.3.2 with two linear elastic materials can be considered. However, now it is assumed that the materials have different volume fractions, with the volume fraction of the first material, which has undergone the transformation strain, denoted as f 1 and the volume fraction of the second material denoted asf 1 1 . The Reuss connection is again trivial-all stresses are zero. In the Voigt scheme, the total stress becomes the sum of stresses multiplied by the corresponding volume fractions. It is easy to show that in the Voigt connection of the materials, Given that E 1 >E 2 , this means that with the increase of the volume fraction, the magnitude of the stress of the first material decreases, while the stress remain compressive. The magnitude of the stress of the second material, in turn, increases with the increase of f 1 . Again, since the Mori-Tanaka composite is expected to be between the Reuss and the Voigt bounds, similar trends are expected for it. This was indeed observed in the simulations, where the decrease of the magnitude of the compressive inclusion stress was observed with the increase of the initial volume fraction (results are not shown in figures). This leads to the increase of the velocity in figure 5(b) with the volume fraction. Furthermore, as the magnitude Figure 5. The dependence of the current volume fraction scaled by the initial volume fraction of the particle (a) and the reaction front velocity (b) on the normalised reaction front position measured from the edge of the particle at various values of the initial volume fraction of the particle.
of the matrix stress increases with the increase of the initial inclusion volume fraction, while the matrix bulk modulus remains constant, the elastoplastic transition in the matrix, which causes the kink around front position of 0.1, starts happening earlier, i.e. smaller transformation strain of the inclusion is required to reach the yield stress of the matrix. Furthermore, it must be noted that the sensitivity of the front kinetics to the initial volume fraction is much more prominent for stiffer matrix materials (results are not shown in figures).
Since the composite undergoes finite deformations, the volume fraction of the inclusion evolves during the front movement within the inclusion. This evolution is plotted in figure 5(a) for the current volume fraction scaled by the initial volume fraction, f/f 0 . It is interesting to note that this quantity highly depends on the initial volume fraction. To understand this, the auxiliary 1D example with two linear elastic materials can be invoked again. Now, however, the strains are the quantities of interest. In the case of Reuss coupling, it is obvious that e e e = = , 0 1 t r a n s f 2 and e e = f . As the Mori-Tanaka composite demonstrates the intermediate behaviour between these two schemes, the qualitative trend for the scaled volume fraction is expected to repeat the qualitative trend for the Voigt scheme. This is confirmed in figure 5(a), where f/f 0 decreases with the increase of f 0 . There are of course two distinct kinks on the plot. The first one is related to the elastoplastic transition of the matrix, the second one is due to the onset of the matrix stress relaxation.

Model application to finite strains-lithiation of Si particles
Having studied some of the physical effects arising from the matrix deformation in the smallstrain setting, it is useful to study the model within the finite-strain setting. For this purpose, lithiation of Si particles inside battery anode matrix material is considered. The parameters related to Si particles (mechanical properties of the transformed and the untransformed phases of the particle, the parameters related to the diffusion process and the chemical reaction) have been identified in the previous study by the authors [22]. The parameters of the matrix material are somewhat more challenging to identify and a literature study of the mechanical properties of the battery anode matrix materials is given below. 4.5.1. Model parameters. Typical lithium-ion battery anodes consist of a mixture of graphite, voids, polymer binder and a liquid electrolyte that allows for the transport of Li ions. As discussed in the introduction, the next-generation battery anodes also contain Si particles which act as an active material for storage of Li ions. During lithiation, Si undergoes large volumetric expansion of up to 300% [5], which leads to mechanical stresses affecting the kinetics of the reactions. Therefore, from the mechanical point of view, the most important effect for the analysis is the stress-affected kinetics inside the particle, while the matrix interferes in this process as a 'boundary condition' for the particle. This justifies the consideration of an average matrix, without making a distinction between graphite, binder or voids. This, however, leads to a difficulty of estimating the effective properties of such averaged matrix.
Since, from the physical point of view, the matrix material is a porous or a granular foamlike structure made of graphite and polymer binder, the chosen constitutive law for the matrix, section 2.3.3, has the volumetric plasticity, which corresponds to the irreversible collapse of the walls in the matrix material, and the volumetric hardening. The deviatoric deformation of the matrix is not explicitly modelled here, although the shear modulus of the matrix participates in the formation of the Eshelby tensor used in the Mori-Tanaka coupling scheme, which assumes that the average deviatoric deformation of the matrix remains elastic. Such constitutive behaviour requires five parameters: the bulk and the shear moduli, the bulk hardening modulus, the power law parameter and the relaxation time.
Literature on the volumetric mechanical behaviour of the graphite-binder anode materials is very limited. In [46], density as a function pressure has been measured for anodes consisting of 90%-92% graphite and 8%-10% PVDF binder (in the experimental papers discussed in this subsection, weight percentages were given). It is easy to convert densitypressure data to volume-stress dependence, where the volumetric stress is implied. Thus, from the data of figure 2 of [46] it is possible to extract the bulk modulus of 56 MPa for the mixture. In [47], density versus pressure has been measured for anodes prepared of three different types of coated/uncoated graphite particles and 10% PVDF binder. Using the same procedure, from the data of table 1 of [47] it is possible to extract the bulk moduli of 81.5, 78 and 153 MPa for three different types of source material. Finally, in [48], density versus pressure has been measured for anodes consisting of 50% graphite 35% PVDF binder and 15% acetylene black. Similarly, from the data of figure 3 of [48] it is possible to obtain the bulk modulus of 7.5 MPa. It can be seen that there is a more than an order of magnitude spread in the average bulk modulus of the anode matrix material. Furthermore, none of the studies characterised the volumetric yield and post-yield behaviour.
The uniaxial compressive deformation of the graphite-binder anode materials including the post-yield behaviour has been characterised in [49]. The prepared samples contained 4% PVDF binder. From the data of figure 2 of [49] it is possible to obtain the Young's modulus of 150 MPa, the hardening modulus of 31 MPa and the yield stress of 67 MPa, although the images of the samples clearly showed the signs of fracture during the post-yield behaviour, which means that the post-yield parameters are not very reliable. Furthermore, in [50] the Young's moduli of materials, which were made of 89% of graphite, 8% PVDF and 3% acetylene black using different production process, were measured to be ranging from 258 to 693 MPa.
Due to a large spread of matrix properties found in literature, which depend on the manufacturing conditions as well as on the composition of the matrix, the parametric study was performed in sections 4.3 and 4.4. For the anode material the bulk modulus K 0 m used in the simulations ranged from 12.5 to 3200 MPa in logarithmic scale. Of course, values of 3200 MPa are well beyond those of the considered materials, but the results are provided to demonstrate the asymptotics of the behaviour. The shear modulus G m was assumed to be 400 MPa. These properties correspond to the Young's modulus ranging from 103 to 1152 MPa, which covers the experimentally observed range. Using the assumption that the volumetric behaviour is similar to the tensile behaviour, the volumetric hardening modulus K h was always taken to be 1/5 of the bulk modulus K 0 m . Since no data on strain-rate sensitivity of the matrix material is available, the power law parameter q m was assumed to be 20, which gives almost strain-rate independent material. Finally, the relaxation time τ m was fitted for each bulk modulus, such that true strain at the yield point was 3% at the strain rate of 2. Effects of matrix stiffness and particle volume fraction on kinetics. The velocity profiles resulting from variation of the matrix stiffness and the particle volume fraction are shown in figure 6. It can be seen that in the case of finite deformations, qualitatively exactly the same effects are observed as in the case of small deformations: (A) the increase of the matrix bulk modulus decreases overall velocity and leads to the arrest of the chemical reaction earlier, (B) the elastoplastic transition in the matrix (which corresponds to the kink in the velocity profile) takes place further from the initial position of the reaction front with the increase of the matrix bulk modulus, (C) there are two asymptotic behaviours with respect to the matrix bulk modulus-first is the approach of the velocity profile to that of the stress-free particle for small matrix bulk modulus, second is the approach a certain profile corresponding to a particle in an incompressible matrix for large matrix bulk modulus, (D) the increase of the particle volume fraction increases overall velocity profile and leads to the arrest of the chemical reaction later, (E) there is one asymptotic behaviour with respect to the particle volume fraction-the velocity profile approaches that of the stress-free particle for large volume fractions, (F) the elastoplastic transition in the matrix is shifted to positions closer to the initial position of the reaction front with the increase of the particle volume fraction. In figure 6, it can also be seen that under current model parameters, the velocity is not significantly affected by the matrix, as the matrix is relatively compliant compared to the particle. This results in a conclusion that currently used graphite-voids-binder matrix materials discussed in section 4.5.1 should affect the kinetics of the chemical reaction only slightly, while the major effect of the stresses on the kinetics results from the build-up of stresses in the particle shell. Since the observed physical effects for the case of finite deformations are qualitatively identical to the effects for the case of small deformations, the analysis of the stresses and strains in the matrix material, section 4.3.2, and stress relaxation, section 4.3.3, remains valid for this case.

Influence of matrix relaxation on kinetics
As discussed above, the elastoplastic transition in the matrix material is manifested as a kink in the kinetics of the reaction front. Such effect should be observed when any change of a deformation mechanism of the matrix takes place. In this subsection, a system with the linear viscoelastic matrix material is studied, where the matrix material can undergo stress Figure 6. The dependence of the reaction front velocity on the normalised reaction front position measured from the edge of the particle when the matrix bulk modulus is varied and f 0 =0.1 (a) and when the particle volume fraction is varied (b) for the case of finite strains.
relaxation. The simulation parameters used in this subsection are summarised in tables F1, F2 and F10. The initial volume fraction has been taken f 0 =0.1. In figure 7, the reaction front velocity profile is plotted for different shear moduli and relaxation times of the matrix. It is clearly seen that the profiles become non-monotonous and now, there are two competing characteristic processes in the system with their own characteristic times: (A) the reaction front retardation due to the material transformation and the build-up of stresses and (B) the relaxation of the matrix stresses.
Higher shear modulus of the matrix corresponds to higher matrix stresses and, therefore, higher sensitivity of the velocity to the stress relaxation mechanism. Obviously, the absolute value of the velocity for a given relaxation time is lower for a higher shear modulus.
As summarised in section 2.3.3, the rheological model of the matrix is similar to the standard linear solid model, i.e. two branches connected in parallel: the spring and dashpot branch and the hardening spring branch. In the simulations, the relaxation time was ranging from 0.1 to 100 s. This means that at the lower limit, the dashpot relaxes almost instantaneously and the matrix is effectively represented by only the hardening spring. At the higher limit, the dashpot stays 'locked' for the most part of the simulation, i.e. the matrix is effectively represented by two parallel springs, and 'unlocks' only when the velocity drops to small values, which corresponds to the simulation time increasing rapidly. Such 'unlocking' for high τ m leads to an interesting effect-in figure 7(b), the velocity rapidly decreases almost linearly to a certain threshold value, beyond which the profile becomes almost horizontal and corresponds to the reaction front slowly creeping towards the centre of the particle. This is the extreme case of the 'kink' in the velocity profile due to change of the deformation mechanism of the matrix.
At the intermediate values of the relaxation time, initially, the profile is following the 'locked' dashpot case, and after the stress relaxation process takes over, the matrix stress rapidly decreases, which removes the surplus pressure from the particle and the velocity increases rapidly. With time, the velocity approaches the 'unlocked' dashpot case. Here terms 'locked' and 'unlocked' are used in an approximate way, as the dashpot relaxation is exponential.

Comparison to the FF simulations for finite strains
Finally, it is useful to perform a similar comparison between the proposed MF model and the FF simulations, as done in section 4.2, but now for the case of finite deformations. In figure 8 an example of such comparison is provided for the traction at the edge of the particle and the reaction front velocity. In the case of the FF model, the total radius of the 'coated' particle is denoted as R C . Thus, the decrease of R 0 /R C represents the approach to the case of an isolated particle in the infinite matrix. In the case of the MF model, such approach takes place when f 0 tends to zero. It can be seen that there is certain discrepancy in the solutions with the FF model having the higher absolute values of stresses in the particle, therefore, lower velocity and earlier reaction arrest. The parameters used for these simulations are given in tables F1, F2 and F11.
The differences are quantitatively significant and result from the MF model not being able to resolve the local stresses around the particle. The selected set of the matrix parameters for the comparison is somewhat intermediate within the range of parameters used in section 4.5.2. For a different set of parameters the discrepancy might change, and it can be speculated the it can increase for a stiffer matrix, due to higher absolute values of the stresses in this case. Due to highlighted differences, the Mori-Tanaka coupling scheme in application to the particles undergoing large mechanochemical transformation must be viewed as the qualitative tool, although for in the case of small strains, the MF and the FF models matched both qualitatively and quantitatively even for nonlinear time-dependent constitutive behaviour of the matrix.

Conclusions
There is a number of advanced engineering heterogeneous materials, which undergo complex localised chemo-mechanical transformations. Since the kinetics of the chemical reactions is stress-affected, accurate and computationally efficient way of accounting for the behaviour of material components (e.g. matrix) is required. In this paper, a multiscale approach based on the MF homogenisation, which captures the movement of the reaction front in spherical particles and accounts for the particle-matrix coupling, has been proposed and examined. This has been achieved by combining a single two-phase particle, in which the kinetics of the reaction front is explicitly described and governed by the chemical affinity tensor, with a matrix using the Mori-Tanaka approach that is generalised for finite deformations and timedependent inelastic materials. Figure 8. The dependence of the traction at the edge of the particle (a) and the reaction front velocity (b) on the normalised reaction front position measured from the edge of the particle. The mean-field model (MF) with the various initial volume fraction of the inclusion is compared to the full-field model (FF) with the various proportion of the inclusion radius to the total radius.
The major advantage of the model proposed in this paper is that it allows performing investigations of the effects related to the interaction of the particle undergoing mechanochemical transformation and the matrix with complex rheology. An application of the model to large volumetric deformation of the Si-based composite anodes of Li-ion batteries showed that in the case when the matrix is viscoelastic or elasto-viscoplastic, the change in the deformation mechanism of the matrix due to the build-up of stresses, caused by the particles' expansion, impacts the kinetics of the chemical reaction with a visible kink in the velocity. When matrix is viscoelastic, there are two competing processes in the system-the retardation of the reaction due to build-up of stresses and the stress relaxation in the matrix. Such competition causes non-monotonous kinetics of the chemical reaction, which follows deceleration-acceleration-deceleration pattern.
The proposed micromechanical model is computationally fast compared to the FF schemes. The model has been compared to the FF simulations in the regime when the particle volume fraction tends to zero. In the case when the particle expansion remained in the smallstrain regime, an excellent match between the MF and the FF models has been observed, even for the elasto-viscoplastic constitutive behaviour of the matrix. It must be emphasised, that the compared quantity has been the velocity of the chemical reaction front, which is strongly sensitive to local fields at the interfaces. Usually, when MF-type and FF models are compared, local quantities demonstrate poorer match than global quantities. Therefore, the achievement of the validity of the proposed model in the small-strain regime becomes of high significance. Similar comparison in the finite-strain regime, however, resulted in an observable discrepancy between the models, which could have been expected due to the restrictive MF assumptions. Although the proposed model cannot qualitatively capture the extreme finite-strain regime of 300% volumetric deformation, for which the validation has been performed, it can still be useful for predicting overall trends with computational efficiency.
There are studies, where the comparison between MF and FF models in the finite-strain regime is performed for effective composite properties, e.g. [24,25], and the MF models are shown to provide good estimates. Here, it has been shown that certain caution must be exercised when MF-type models are generalised to finite deformations, as the predictions of the local fields at the interfaces, which can be important to some physical problems, such as mechanochemistry, might somewhat deviate from the FF models. Hence, reliable simulation of the matrix-particle interaction in the finite-strain regime still requires FF simulations, while in the small-strain regime, the proposed model is an excellent alternative to the FF models.

·˙)
When equation (4) is differentiated, the expression for the velocity gradient can be obtained, where L a and L m are the velocity gradients of the inclusion and the matrix, respectively. For purely volumetric deformations, Taking the symmetric part of L ⟨ ⟩ and substituting expression for χ, equation (7) is obtained.

Appendix B. Incremental formulation of volume-averaging
The time derivative of the deformation gradient can be represented as where subscripts n andn 1 denote the current and the previous time increments, respectively. This is applicable to any deformation gradient-the macroscopic, the inclusion and the matrix. Furthermore, such representation introduces two distinct configurations-the current and the previous time increment configurations. Therefore, it is possible to construct a mapping = D -F F F , 4 9 Here, the infinitesimal linear strain e D , which is applied to the previous time increment configuration, is introduced. Using this, it is easy to see that equations (6) and (7)

Appendix C. Tensors in the Mori-Tanaka scheme
The Mori-Tanaka homogenisation scheme relies on the solution of the Eshelby's problem (ellipsoidal inclusion in the infinite matrix). In the general case, the inclusion has different stiffness tensor than the matrix, it has undergone stress-free transformation strain (e.g. thermal expansion) and, furthermore, there is strain applied at the infinity. The solution to this problem is well known and can be found in many textbooks, e.g. [51,52] where e 1 is the inclusion strain, e ¥ is the strain applied at infinity 9 , e * is the stress-free transformation strain of the inclusion, C 4 1 and C 4 2 are stiffness tensors of the inclusion and the matrix, respectively, and S 4 is the Eshelby tensor. The idea of the IMT scheme is to apply the CMT equations to the incremental form of the volume-averaging equations. Since the matrix and the inclusion undergo finite deformations in this paper, the incremental volume-averaging equations are non-trivial and are derived in appendix B. Thus, the CMT scheme it applied to equations (52) and (53).
In the CMT scheme, it is assumed that the Eshelby's solution for the inclusion is valid, equation (54), in which e ¥ is equated to to the matrix strain of the CMT composite. Therefore, using equation (54) where the notation was deliberately changed from subscripts '1' and '2' to superscripts 'e' and 'm', respectively, to distinguish the Eshelby's problem from the IMT scheme. Here, C they become the current tangent stiffness tensors. The Eshelby tensor, S 4 , is a function of C 4 m . Here, e D * denotes the incremental stress-free expansion of the inclusion, i.e. such expansion of the inclusion within the increment that the increment of the inclusion stress is zero.
By the assumptions of the CMT scheme, section 2.1, it is always possible to relate e D m to e D ⟨ ⟩ and e D * using the concentration tensors: When limit Δt → 0 is taken, these equations become equations (11) and (12). The Eshelby tensor for spherical inclusions can be found in many textbooks and is given by n The first Piola-Kirchhoff stress tensor of the untransformed material used in equation (22)