**1. Introduction**

Currently, lithium-ion batteries (LIB) are one of the most common devices for energy storage [1]. Many works are devoted to understanding the lithium intercalation/deintercalation mechanisms into electrodes and to optimization of LIB components for better electrochemical characteristics (see reviews [2,3]). The most widely accepted electrochemical model of LIB is the so-called pseudo-two-dimensional (P2D) model [4,5], where the intercalation and deintercalation in electrode particles, ion transport in electrolyte and separator are described in terms of the normal di ffusion equations. The generalization to the anomalous di ffusion case is proposed in [6]. However, in several electrode materials, the mutual solubility of the lithiated and delithiated phases is very low, and two phases coexist in a wide range of the state of charge (SOC). Such behavior is typical for charging and discharging in cathode based on lithium iron phosphate and in anode materials based on lithium titanate [7]. For this reason, lithiation and delithiation of the cathode particles can be associated with reversible phase transition in quasi-binary system, where the intercalant concentration can vary in the range from 0 to 1 during charging or discharging. The corresponding theoretical approach is based on the phase-field theory [8–16] that is commonly employed in studies on phase transitions [17–19].

The fact that the electrodes consist of particles of various sizes is often neglected. For example, the classic single particle model (SPM) and P2D model assume the electrodes consisting of spherical particles of identical size. Nonetheless, electrodes are typically made of porous materials with particles of di fferent sizes and shapes. As a result, the medium in LIB is a highly heterogeneous system. Neglecting the particle size distribution (PSD) largely underestimates the capacity in the case of elevated C-rates [20], and the SPM approaches poorly describe the battery performance at higher current densities. PSD is an important factor in the degradation of batteries [21,22], and it may be altered during battery operation due to cracking or agglomeration of particles [21]. In [23], the multiple-particle model indicates that the PSD broadening may lead to higher values of volumetric capacity and energy density. On the other hand, this broadening can cause amplification of electrode polarization [23].

There are several papers (see e.g., [20,23–26]) where the particle size dispersion is considered within the di ffusion approach. In most cases, the actual PSD is described by finite number of particle groups [20,25,26]. The multiple-particle approach meets computational di fficulties related to the need to solve di ffusion equations for each particle bin [26].

In the present paper, we extend the phase-field model based on the Cahn-Hilliard equation to simulate lithium intercalation dynamics in a cathode with particles of distributed size. Due to PSD, during charging/discharging, number of active electrode particles and their interfacial flux density depend on time. The boundary condition with time-dependent flux density for the Cahn-Hilliard equation is formulated from the universal laws established for speed of front propagation and evolution of single-phase concentration in individual particles. Calculation of the time-dependentflux density is performed by the developed general approach that can be employed for analysis of intercalation process in particle ensemble with arbitrary PSD.

#### **2. Phase-Field Intercalation in a Single Spherical Particle**

We start with a simplified phase-field model for a single submicron particle under galvanostatic condition. Expecting to establish some general rules of phase-field intercalation in a spherical particle, we analyze the dynamics for wide ranges of particle size and C-rate. The process of intercalation is simulated in terms of the Cahn-Hilliard (CH) equation [14,16,18,19]:

$$\frac{\partial \mathcal{C}}{\partial t} = M \nabla^2 \left[ \frac{\partial f}{\partial \mathcal{C}} - \kappa \nabla^2 \mathcal{C} \right] \tag{1}$$

where *c* ≡ *<sup>c</sup>*(*<sup>r</sup>*, *t*) is concentration field of intercalating atoms that defines the local state of charge of a cathode particle, *M* is the mobility, *f*(*c*) is the free energy of mixing in a quasi-binary system, κ is the gradient energy coe fficient. In this study, the free energy of mixing *f*(*c*) of quasi-binary solid solution is considered in the regular solution approximation [17–19]:

$$f(c) = \Omega c(1 - c) + k\_B T [c \ln c + (1 - c) \ln(1 - c)]\_+$$

Ω is the interaction parameter, *T* is the temperature, and *kB* is the Boltzmann constant. The influence of elastic strain is neglected. Assuming spherical symmetry of cathode particles, we use the Laplace operator of the form: ∇<sup>2</sup> = ∂2 ∂*r*<sup>2</sup> + 2 *r*∂ ∂*r*.

The insertion and extraction of intercalating atoms are simulated at the temperature of 300 K. The mobility of intercalated atom *M* in a cathode particle is related to the di ffusion coe fficient through the Einstein relation, *M* = *<sup>D</sup>*/(*kBT*). The di ffusion coe fficient of intercalating atom (e.g., lithium) for di fferent cathode materials is *D* = 10−13–10−<sup>15</sup> m<sup>2</sup>/s at the considered temperature [4,8,13]. The interaction parameter Ω for di fferent cathode materials lies in the range of 0.059–0.193 eV/atom [8–16]. The gradient energy coe fficient can be related to the width of equilibrium profile [17,18] that is usually of several nanometers [8,10]. In this study, we employ the di ffusion coe fficient *D* = 10−<sup>14</sup> m<sup>2</sup>/s, interaction parameter Ω = 0.115 eV/atom. The gradient energy coe fficient is usually written as κ = ηΩ*r*<sup>2</sup> 0 [8,14,17,18], where *r*0 is the intermolecular distance and η is the dimensionless coe fficient depending on form of interaction potential [18]. In this study, the value of gradient energy coe fficient is κ = 0.228 eVnm2.

The interaction parameter allows us to calculate phase diagram of the system (Figure 1). The equilibrium values of intercalating atom concentration are *cb*1 = 0.013 and *cb*2 = 0.987 at the temperature of 300 K. The unstable region lies in the interval between *cs*1 = 0.129 and *cs*2 = 0.871.

**Figure 1.** Phase diagram for regular binary solution with interaction parameter of Ω = 0.115 eV/atom.

The CH Equation (1) is solved under natural boundary conditions corresponding to the galvanostatic mode:

$$\left. \frac{\partial \mathbf{c}}{\partial r} \right|\_{0, \mathbb{R}\_0} = 0, \left. \frac{\partial \mu}{\partial r} \right|\_0 = 0, \left. \frac{\partial \mu}{\partial r} \right|\_{\mathbb{R}\_0} = -\frac{j\_r}{M}. \tag{2}$$

Here, μ = ∂ *f*/∂*c* − κ∇2*c* is the chemical potential defining the flux density of intercalating atoms. Initial distribution of the atoms is uniform and is characterized by equilibrium values of concentration, i.e., *c*0 = *cb*1 for insertion and *c*0 = *cb*2 for extraction process. The CH Equation (1) with boundary conditions (2) is solved numerically by using the explicit Euler time integration scheme [27].

Figure 2 demonstrates the evolution of concentration profile of intercalating atoms computed with the introduced values of the model parameter. The insertion and extraction processes are simulated for the cathode particle with radius of *R*0 = 0.1 μm under C-rate of 1C and 10C. Also, we simulated the insertion and extraction dynamics for cathode particle with radius *R*0 = 0.2 μm under C-rate equal to *C*/2.

Analyzing the results of insertion and extraction of intercalating atoms in cathode particles, we can identify some important features of these processes. At the very beginning of the processes, the concentration of intercalant decreases uniformly over the particle volume (see lines 1 and 2 in Figure 2a–f). At this stage, the particle corresponds to the single-phase pattern. The uniformity is explained by fast equilibration of atom distribution due to high diffusion coefficient and small size of the particle. The stage continues until the nearest value of metastability limit (*cs*1 or *cs*2) is achieved. At this stage, the average intercalant concentration *c* changes linearly with time (Figure 3). The approximate dependence can be easily obtained from the CH Equation (1):

$$
\overline{x} = c\_0 - \frac{3j\_r t}{R\_0} \,\text{\,\, \,}\tag{3}
$$

that agrees well with the results of direct solution of CH Equation (Figure 3). The relation (3) is valid at the time interval 0 ≤ *t* ≤ *ts*. From Equation (3) we can calculate the duration of the single-phase stage as *ts* = (*cb*2 − *cs*2)*R*0/3*jr* for extraction and *ts* = (*cs*1 − *cb*1)*R*0/3*jr* for insertion process.

**Figure 2.** Concentration profiles of intercalating atoms for different time points at *T* = 300 K. Corresponding radius of particle and *C*-rate are: (**a**) 0.1 μm, 1C insertion, (**b**) 0.1 μm, 1C extraction, (**c**) 0.1 μm, 10C insertion, (**d**) 0.1 μm, 10C extraction, (**e**) 0.2 μm, 0.5C insertion, (**f**) 0.2 μm, 0.5C extraction. The lines 3, 4, 9 and 10 correspond to passing from single-phase stage to double-phase stage and vice versa.

**Figure 3.** Kinetics of average concentration at the single-phase stage of insertion (1–3) and extraction (4–6) processes. Points are the solution of the CH equation, lines correspond to the linear dependence (3).

After that the concentration profile dramatically changes. Interfacial region depletes in very short time interval and the particle passes to double-phase pattern (see lines 3 and 4 in Figure 2a–f). At this stage, the intercalant redistributes over the particle volume and passes to the nearest equilibrium composition (see lines 4 and 5 in Figure 2a–f). Transition from single-phase to double-phase stage takes place for the short time interval estimated as ~0.2 s and ~0.5 s for particles with size of 0.1 μm (1C) and 0.2 μm, respectively. The longer time interval for larger particle is explained by the presence of additional concentration wave moving toward the particle center (see line 4 in Figure 2a–f) during equilibration of concentration profile. Overpotential of electrode depends on interfacial concentration of intercalating atoms [4,5,8]. Therefore, it is expected that passing from single-phase to double-phase regime can cause the abrupt change of overpotential. The effect could be used for determination of metastability limit.

The next stage is characterized by motion of the concentration wave from the particle interface to the center (lines 5–9 in Figure 2a–f). The position of concentration wave can be associated with coordinate *Rh*(*t*) corresponding to concentration *ch* = (*cb*1 + *cb*2)/2. The position of the wave front *Rh*(*t*) satisfies the linear dependence

$$\frac{R\_0^3 - R\_h^3(t)}{R\_0^2} = \frac{3|j\_r|t}{(c\_{b2} - c\_{b1})'} \tag{4}$$

obtained from the solution of the CH Equation (1). The Equation (4) is valid for *t* > *ts* (after the single-phase stage). The results of simulation by the CH Equation (Figure 2a–f) confirm Equation (4) (see Figure 4). The end of the double-phase stage corresponds to the time point *td* = *<sup>R</sup>*0(*cb*2 − *cb*1)/3*jr* for both insertion and extraction processes. The relation can be employed for calculation of the flux density *jr* corresponding to 1C-rate, if duration of charging or discharging process equals to *td* = 3600 s.

At the end of insertion or extraction process, the system passes over unstable region to another single-phase state that corresponds to fully charged or discharged state of a particle (lines 9 and 10 in Figure 2a–f). The process of insertion or extraction is stopped, when the corresponding equilibrium composition is achieved.

**Figure 4.** Kinetics of concentration wave propagation obtained by the CH Equation (points). Lines represent the linearized dependence (4). Size of particles and applied C-rate are indicated in Figure 2.

Contribution of the single-phase and double-phase stages can be estimated by ratio τ*s* = *ts*/*td*. This ratio can be easily found from Equations (3) and (4):

$$
\tau\_{s1} = \frac{c\_{s1} - c\_{b1}}{c\_{b2} - c\_{b1}}, \ \tau\_{s2} = \frac{c\_{b2} - c\_{s2}}{c\_{b2} - c\_{b1}}.
$$

for charging and discharging, respectively. The value of τ*s* is defined by the phase diagram and depends on temperature only. In case of symmetric phase diagram, this parameters are equal to each other τ*s*1 = τ*s*2. Calculation of the ratio in the range 250–400 K gives contribution of single-phase stage about 10–15%. Simulation of insertion and extraction for lower diffusion coefficient reveals deviation from concentration profiles at the single-phase stage only if the diffusion coefficient is less than ~10−<sup>17</sup> m<sup>2</sup>/s.

Remarkably, variation of insertion and extraction rates in the range 1–10 C at the temperature of 300 K does not modify the two-stage mechanism and relations (3) and (4) remain valid (see Figures 2–4). Thus, Equations (3) and (4) can be used for determination of the composition profile at any time point of insertion or extraction process.

#### **3. Intercalation in Particles of Distributed Size**

Real cathodes consist of particles of various sizes and can be characterized by certain PSD. Under the assumption that the current is equally distributed over the active surface, smaller particles are charged or discharged for shorter time interval and withdraw earlier from the intercalation process. Larger particles have to balance the total current change in the galvanostatic mode. Therefore, flux density at the particle surface turns to depend on time *jr* = *jr*(*t*).

Let us modify the obtained relations with respect to size distribution of cathode particles. If the flux density at the particle interface depends on time after integration of CH equation we obtain for single-phase and double-phase regimes:

$$\begin{aligned} \overline{\boldsymbol{c}}(t) &= \mathbf{c}\_0 - \frac{3}{\mathbb{R}\_0} \int\_0^t \boldsymbol{j}\_r(t)dt, \ 0 \le t \le t\_s, \\\overline{\boldsymbol{R}}\_0^3 - \frac{\mathbf{c}\_k^3(t)}{\mathbf{c}\_0^2} &= \frac{3}{(\mathbf{c}\_{b2} - \mathbf{c}\_{b1})} \int\_0^t \left| \boldsymbol{j}\_r \right| dt, \ t\_s < t < t\_d. \end{aligned} \tag{5}$$

Whereas the total amount of intercalating atoms is conserved, the relation for total time of insertion or extraction process can be written in the form

$$\frac{1}{3}R\_0(c\_{b2} - c\_{b1}) = \int\_0^{t\_d} |j\_\prime(t)|dt.\tag{6}$$

Equation (6) can be used for definition of threshold value *R*min defining the minimal size of particle that can participate in intercalation process

$$\frac{d\mathcal{R}\_{\text{min}}}{dt} = \frac{\mathfrak{d}|j\_{\mathcal{I}}(t)|}{c\_{\mathfrak{k}2} - c\_{\mathfrak{k}1}}.\tag{7}$$

Then we assume that particles with *R* ≤ *R*min are excluded from the intercalation process (they are fully charged/discharged). The total current redistributes over the interface of residual particles with *R* > *R*min. Therefore, the galvanostatic mode can be determined by the following equation:

$$j\_r(t) \int\_{R\_{\min}}^{\infty} w(R) R^2 dR = j\_0 \{ R^2 \},\tag{8}$$

where *j*0 is the flux density at the beginning of intercalation process and *w*(*R*) is the size distribution function of cathode particles, *R*<sup>2</sup> is the squared radius averaged over the whole ensemble of cathode particles.

Combining Equations (7) and (8) leads to the equation for the threshold value *R*min

$$\frac{d\mathcal{R}\_{\text{min}}}{dt} = \frac{3j\_0 \{\mathcal{R}^2\}}{\left(c\_{b2} - c\_{b1}\right) \int\_{\mathcal{R}\_{\text{min}}}^{\infty} w(R)R^2 d\mathcal{R}}.\tag{9}$$

Let us assume that ensemble of cathode particles is described by the gamma distribution

$$w(R) = \frac{R^{m-1}}{a^m \Gamma(m)} \exp\left(-\frac{R}{a}\right) \tag{10}$$

that is usually employed as PSD for different cathode materials [28,29]. Here *a* and *m* are the constant parameters and <sup>Γ</sup>(*m*) is the gamma function. The expectation value and dispersion of the distribution are *R* = *ma* and σ2 = *ma*2, respectively.

Integration of Equation (9) with respect to PSD (10) gives the implicit time-dependence of *R*min

$$\frac{a^3}{\Gamma(m)} \left[ \Gamma(m+3) + \frac{R\_{\text{min}}}{a} \Gamma \left( m + 2, \frac{R\_{\text{min}}}{a} \right) - \Gamma \left( m + 3, \frac{R\_{\text{min}}}{a} \right) \right] = \frac{3j\_0 \left\langle \mathbf{R}^2 \right\rangle t}{\left( \mathbf{c}\_{b2} - \mathbf{c}\_{b1} \right)}. \tag{11}$$

.

The threshold radius *R*min alters from zero to infinity, and the charging/discharging time is determined by

$$t\_{\max} = \frac{a^3 (c\_{b2} - c\_{b1}) \Gamma(m+3)}{3 j\_0 \langle \mathbb{R}^2 \rangle \Gamma(m)}$$

Figure 5 shows the time-dependence of threshold particle radius *R*min in dimensionless coordinates for different values of parameter *m*. The results in these coordinates do not depend on parameter *a* (see Equation (11)).

**Figure 5.** Dynamics of minimal radius *R*min of active particle in process of intercalation. Lines are the results of analytical solution (11) for different values of parameter *m*: (1) *m* = 2, (2) *m* = 4, (3) *m* = 8.

Analyzing Figure 5, we conclude that *R*min depends linearly on time at the beginning of insertion or extraction process that can be described by dependence

$$\frac{R\_{\text{min}}}{\Gamma(\text{R})} = \frac{\Gamma(m+3)}{m\Gamma(m+2)} \cdot \frac{t}{t\_{\text{max}}}.\tag{12}$$

The time interval corresponding to linear dependence is characterized by a constant flux density *jr*. After that the flux density *jr* of intercalating atoms grows very quickly and goes to infinity when time approaches to *t*max.

Taking Equation (7) into account, one can transform Equations (5) into

$$\overline{c}(t) = c\_0 - (c\_{b2} - c\_{b1})\frac{R\_{\text{min}}(t)}{R\_0}, \ 0 \le t < t\_{s\prime} \tag{13a}$$

$$\frac{R\_0^3 - R\_h^3(t)}{R\_0^2} = R\_{\min}(t), \quad t > t\_s. \tag{13b}$$

These equations describe the dynamics of intercalant distribution over the electrode particle at the single-phase and double-phase stage with respect to PSD.

Let us consider the impact of the PSD parameters on the extraction process from the particle ensemble characterized by PSD (10). Here, *R* = 0.1 μm and *m* = 4. We solve the CH Equation (1) for particle with size of *R*0 = 0.2 μm and time-dependent interfacial flux density *jr*(*t*) that defines the boundary condition (2).

The time-dependence of flux density *jr*(*t*) can be calculated by Equations (7) and (11). To represent *<sup>R</sup>*min(*t*) given by Equation (11) in explicit form, we choose the following approximation

$$\frac{R\_{\rm min}}{\langle R \rangle} = \sum\_{m=0}^{n} a\_{m} \left( \frac{t}{t\_{\rm max}} \right)^{2m+1} + b \, \text{arctanh} \left( \frac{t}{t\_{\rm max}} \right) \,. \tag{14}$$

where *n* is the order of approximation, *am* and *b* are fitting parameters. Using this approximation, we obtain the explicit dependence of the interfacial flux density for *n* = 0 in the form

$$j\_r(t) = j\_0 \left( 0.512826 + 0.487174 / (1 - (t/t\_{\text{max}})^2) \right).$$

The initial flux density *j*0 corresponds to C-rate of 1C for cathode particle with average size *R* .

The result of simulation of the extraction process from an individual particle of the ensemble is shown in Figure 6. The radius of this particle is chosen equal to *R*0 = 0.2 μm. The concentration profiles at different times are similar to Figure 2f. Analyzing Figure 6, one can conclude that solution of the CH Equation (1) with time-dependent interfacial flux density agrees well with the general Equation (13b) at the double-phase stage. The size distribution function of the cathode particles causes decrease of the total time of extraction that corresponds to time interval Δ*t* in Figure 6.

**Figure 6.** The time-dependence of concentration wave position *Rh* in dimensionless units for cathode particle with the size of *R*0 = 0.2 μm at the double-phase stage. The particle is the part of ensemble described by PSD given by Equation (10) with *R* = 0.1 μm and *m* = 4. Solid line 1 corresponds to solution of Equation (11). Dashed line 2 is the linear dependence corresponding to constancy of the interfacial flux density. Points 3 and 4 correspond to moments of start and finish of double-phase stage for the particle. Triangle points correspond to solution of the CH equation.

The time-dependent flux density is approximately constant at the single-phase stage, therefore the dependence of average concentration *c* at this stage is almost identical to that in Figure 2f. Generally, the time-dependent flux density can influence the average concentration *c* dynamics at the single-phase stage for large particles. The large particles are expected to deviate from the linear dependence that is observed under constant flux density (see Figure 3).
