Effect of the Particle Size Distribution on the Cahn-Hilliard Dynamics in a Cathode of Lithium-Ion Batteries

: The phase-field model based on the Cahn-Hilliard equation is employed to simulate lithium intercalation dynamics in a cathode with particles of distributed size. We start with a simplified phase-field model for a single submicron particle under galvanostatic condition. We observe two stages associated with single-phase and double-phase patterns typical for both charging and discharging processes. The single-phase stage takes approximately 10%–15% of the process and plays an important role in the intercalation dynamics. We establish the laws for speed of front propagation and evolution of single-phase concentration valid for different sizes of electrode particles and a wide range of temperatures and C-rates. The universality of these laws allows us to formulate the boundary condition with time-dependent flux density for the Cahn-Hilliard equation and analyze the phase-field intercalation in a heterogeneous cathode characterized by the particle size distribution. intercalation dynamics in a cathode with particles of distributed size. Considering the simplified phase-field intercalation in a single submicron spherical particle under galvanostatic condition, we established universal behavior observed for different sizes of electrode particles and a wide range of temperatures and C-rates. The universal rules were formulated for speed of front propagation and evolution of single-phase concentration. Two stages associated with single-phase and double-phase patterns are typical for both charging and discharging processes. The single-phase stage takes approximately 10%–15% of the process and plays an important role in the intercalation dynamics.


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-twodimensional (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 diffusion equations. The generalization to the anomalous diffusion 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][9][10][11][12][13][14][15][16] that is commonly employed in studies on phase transitions [17][18][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 different 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][24][25][26]) where the particle size dispersion is considered within the diffusion approach. In most cases, the actual PSD is described by finite number of particle groups [20,25,26]. The multiple-particle approach meets computational difficulties related to the need to solve diffusion 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 timedependentflux density is performed by the developed general approach that can be employed for analysis of intercalation process in particle ensemble with arbitrary PSD.

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]: where   , c c t  r 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 coefficient. In this study, the free energy of mixing   f c of quasi-binary solid solution is considered in the regular solution approximation [17][18][19]:  is the interaction parameter, T is the temperature, and B k 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: 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 diffusion coefficient through the Einstein relation, The diffusion coefficient of intercalating atom (e.g., lithium) for different cathode materials is D = 10 −13 -10 −15 m 2 /s at the considered temperature [4,8,13]. The interaction parameter  for different cathode materials lies in the range of 0.059-0.193eV/atom [8][9][10][11][12][13][14][15][16]. The gradient energy coefficient 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 diffusion coefficient D = 10 -14 m 2 /s, interaction parameter Ω = 0.115 eV/atom. The gradient energy coefficient is usually written as 2 0 r     [8,14,17,18], where 0 r is the intermolecular distance and  is the dimensionless coefficient depending on form of interaction potential [18]. In this study, the value of gradient energy coefficient is 0.228   eVnm 2 . The interaction parameter allows us to calculate phase diagram of the system (Figure 1). The equilibrium values of intercalating atom concentration are  The CH Equation (1) is solved under natural boundary conditions corresponding to the galvanostatic mode: Here,  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 R0 = 0.1 µm under C-rate of 1C and 10C. Also, we simulated the insertion and extraction dynamics for cathode particle with radius R0 = 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 ( that agrees well with the results of direct solution of CH Equation (Figure 3). The relation (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 doublephase 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 obtained from the solution of the CH Equation (1). The Equation (4)  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 2 a-f). The process of insertion or extraction is stopped, when the corresponding equilibrium composition is achieved. This ratio can be easily found from Equations (3) and (4): for charging and discharging, respectively. The value of 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.

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   r r j j 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: 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 Equation (6) can be used for definition of threshold value min R defining the minimal size of particle that can participate in intercalation process Then we assume that particles with min R R  are excluded from the intercalation process (they are fully charged/discharged). The total current redistributes over the interface of residual particles with min R R  . Therefore, the galvanostatic mode can be determined by the following equation: where 0 j is the flux density at the beginning of intercalation process and ) (R w is the size distribution function of cathode particles, 2

R
is the squared radius averaged over the whole ensemble of cathode particles. Combining Equations (7) and (8) (9) Let us assume that ensemble of cathode particles is described by the gamma distribution (10) that is usually employed as PSD for different cathode materials [28,29]. Here a and m are the Integration of Equation (9) with respect to PSD (10) gives the implicit time-dependence of min R The threshold radius min R alters from zero to infinity, and the charging/discharging time is      Figure 5 shows the time-dependence of threshold particle radius min R in dimensionless coordinates for different values of parameter m . The results in these coordinates do not depend on parameter a (see Equation (11)).
The time interval corresponding to linear dependence is characterized by a constant flux density r j . After that the flux density r j of intercalating atoms grows very quickly and goes to infinity when time approaches to max t .
Taking Equation (7) into account, one can transform Equations (5) into 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, µm . We solve the CH Equation (1) for particle with size ofR0 = 0.2 µm and time-dependent interfacial flux density ) (t j r that defines the boundary condition (2). The time-dependence of flux density ) (t j r can be calculated by Equations (7) and (11). To represent ) ( min t R given by Equation (11) in explicit form, we choose the following approximation , arctanh max 0 (14) where n is the order of approximation, m a and b are fitting parameters. Using this approximation, we obtain the explicit dependence of the interfacial flux density for 0  n in the form The initial flux density j0 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 toR0 = 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. 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).

Conclusions
We extended the phase-field model based on the Cahn-Hilliard equation to the lithium intercalation dynamics in a cathode with particles of distributed size. Considering the simplified phase-field intercalation in a single submicron spherical particle under galvanostatic condition, we established universal behavior observed for different sizes of electrode particles and a wide range of temperatures and C-rates. The universal rules were formulated for speed of front propagation and evolution of single-phase concentration. Two stages associated with single-phase and double-phase patterns are typical for both charging and discharging processes. The single-phase stage takes approximately 10%-15% of the process and plays an important role in the intercalation dynamics.
The particle size distribution causes the time-dependence of the interfacial flux density under galvanostatic conditions. The analytical results for a particle ensemble are presented for the case of the gammadistribution of particle sizes, and can be easily generalized to other PSDs. The universality of the established laws allowed us to formulate the boundary condition with time-dependent flux density for the Cahn-Hilliard equation. Numerical solutions of the corresponding boundary-value problem agree well with the obtained universal relations.
Thus, we obtained thegeneral analytical relations describing the intercalation process in ensemble of submicron electrode particle characterized by the arbitrary PSD. Other approaches based on the diffusion equation for intercalating atoms can consider finite number of particle bins only and have higher computational cost.
The Cahn-Hilliard equation and obtained approximate relations on the phase-field intercalation can be used for modification of the well-known SPM and P2D models of LIBs.