Computational Modeling of Heterogeneity of Stress, Charge, and Cyclic Damage in Composite Electrodes of Li-Ion Batteries

Figure 2 shows the typical microstructural of a commercial NMC cathode and the sketch of the half-cell model. The microstructure of NMC (Fig. 2a) was imaged using a scanning electron microscopy (SEM, T330, JEOL). Figure 2b shows the cross-section of the NMC cathode after mechanical polishing (EcoMet 3000, Buehler). Grit sizes of 3 and 1 μm were used for the initial polishing and then 0.05 μm colloidal silica suspension was used for final polishing until a mirror-like finish was obtained.^45,46 The active particles are embedded in the porous matrix of carbon additives and polymeric binders. Upon charging, Fig. 2c, Li extracts from the NMC lattice and splits into Li-ions and electrons governed by the charge transfer kinetics on the particle surface. Electrons flow through the conductive carbon-binder network to the cathode current collector, and transport to the anode side through the external circuit. To maintain electroneutrality, Li-ions diffuse through the electrolyte in the pore space of electrodes and separator and reach the anode side to neutralize the electrons. The above processes are reversed if the cell is under discharging. The electrochemical processes were well described by the theoretical framework of Doyle et al.⁴⁷ In our model, we use the subscript s and l to denote the solid active cathode particles and the liquid electrolyte, respectively, to distinguish the different domains.

Zoom In Zoom Out Reset image size

Electrolyte

Li-ion transport in the electrolyte follows the mass conservation:

$$\begin{eqnarray}&&\displaystyle \frac{\partial {c}_{l}}{\partial t}+{\rm{\nabla }}\cdot {{\boldsymbol{N}}}_{l}=0,\end{eqnarray} \tag{ 1 }$$

where ${c}_{l}$ is Li-ion concentration in the electrolyte, ${{\boldsymbol{N}}}_{l}$ is the Li-ion flux including the contributions of both diffusion and migration:

$$\begin{eqnarray}&&{{\boldsymbol{N}}}_{l}=-{D}_{l}{\rm{\nabla }}{c}_{l}+\displaystyle \frac{{{\boldsymbol{i}}}_{l}{t}_{+}}{F},\end{eqnarray} \tag{ 2 }$$

where ${D}_{l}$ represents the Li-ion diffusivity in the electrolyte, ${{\boldsymbol{i}}}_{l}$ the electrolyte current density, ${t}_{+}$ the transference number of cations, and $F$ the Faraday's constant. The electrolyte current density was derived by the prior work.⁴⁴

$$\begin{eqnarray}&&{i}_{l}=-{\kappa }_{l}{\rm{\nabla }}{\phi }_{l}+\displaystyle \frac{2{\kappa }_{l}RT}{F}\left(1+\displaystyle \frac{\partial lnf}{\partial ln{c}_{l}}\right)\left(1-{t}_{+}\right){\rm{\nabla }}ln{c}_{l},\end{eqnarray} \tag{ 3 }$$

where ${\kappa }_{l}$ is the electrolyte conductivity, ${\phi }_{l}$ the electric potential in the electrolyte, $R$ the gas constant, $T$ the temperature, and $f$ the mean molar activity coefficient of the electrolyte. The conservation of charge requires

$$\begin{eqnarray}&&{\rm{\nabla }}\cdot {{\boldsymbol{i}}}_{l}=0.\end{eqnarray} \tag{ 4 }$$

Active particles and conductive matrix

Li concentration within the active particle, ${c}_{s},$ is determined by

$$\begin{eqnarray}&&\displaystyle \frac{\partial {c}_{s}}{\partial t}+{\rm{\nabla }}\cdot {{\boldsymbol{N}}}_{s}=0.\end{eqnarray} \tag{ 5 }$$

The flux of Li, ${{\boldsymbol{N}}}_{s},$ is given by

$$\begin{eqnarray}&&{{\boldsymbol{N}}}_{s}=-M{c}_{s}{\rm{\nabla }}\mu ,\end{eqnarray} \tag{ 6 }$$

where $M$ is the mobility of Li, and $\mu$ is the chemical potential of Li in the host. When the stress effect on the chemical potential is taken into consideration, the chemical potential of Li is written as:^19–23

$$\begin{eqnarray}&&\mu ={\mu }_{0}+RTln\left(\displaystyle \frac{{c}_{s}}{{c}_{s\_\max }}\right)-{\rm{\Omega }}{\sigma }_{h},\end{eqnarray} \tag{ 7 }$$

where ${\mu }_{0}$ is a reference value, ${c}_{s\_\max }$ the maximum Li concentration in the active materials, ${\rm{\Omega }}$ the partial molar volume of Li in the host, and ${\sigma }_{h}$ the hydrostatic stress. Substituting Eq. 7 into 6, Eq. 5 can be rewritten as follows with a definition ${D}_{s}=MRT$ as the diffusion coefficient,

$$\begin{eqnarray}&&\displaystyle \frac{\partial {c}_{s}}{\partial t}+{\rm{\nabla }}\left(-{D}_{s}{\rm{\nabla }}{c}_{s}+\displaystyle \frac{{D}_{s}{\rm{\Omega }}{c}_{s}}{RT}{\rm{\nabla }}{\sigma }_{h}\right)=0.\end{eqnarray} \tag{ 8 }$$

The electric current density in particle and matrix follows the Ohm's law:

$$\begin{eqnarray}&&{{\boldsymbol{i}}}_{s}=-{\kappa }_{s}{\rm{\nabla }}{\phi }_{s},\end{eqnarray} \tag{ 9 }$$

where ${\kappa }_{s}$ is the electrode conductivity and ${\phi }_{s}$ is the potential in the composite electrode. Due to the charge balance, the electrode current density follows

$$\begin{eqnarray}&&{\rm{\nabla }}\cdot {{\boldsymbol{i}}}_{s}=0.\end{eqnarray} \tag{ 10 }$$

Electrode-electrolyte interface

The charge transfer kinetics occurs at the interface between the electrolyte and the active materials and generates the local current. The local current density, ${i}_{loc},$ is often described by the Butler-Volmer equation:

$$\begin{eqnarray}&&{i}_{loc}={i}_{0}\left(\exp \left(\displaystyle \frac{{\alpha }_{a}F\eta }{RT}\right)-\exp \left(-\displaystyle \frac{{\alpha }_{c}F\eta }{RT}\right)\right),\end{eqnarray} \tag{ 11 }$$

where ${\alpha }_{a}$ represents the anodic transfer coefficient, ${\alpha }_{c}$ the cathodic transfer coefficient, and ${i}_{0}$ the concentration-dependent exchange current density given by

$$\begin{eqnarray}&&{i}_{0}=F{\left({k}_{c}\right)}^{{\alpha }_{a}}{\left({k}_{a}\right)}^{{\alpha }_{c}}{\left({c}_{s\_{\rm{\max }}}-{c}_{s}\right)}^{{\alpha }_{a}}{\left({c}_{s}\right)}^{{\alpha }_{c}}{\left(\displaystyle \frac{{c}_{l}}{{c}_{l\_{\rm{ref}}}}\right)}^{{\alpha }_{a}},\end{eqnarray} \tag{ 12 }$$

where ${c}_{l\_ref}$ is the reference Li-ion concentration in the electrolyte, ${c}_{s\_surf}$ the Li concentration at the electrode surface, and ${k}_{a}$ and ${k}_{c}$ are the rate constants for the anodic and cathodic reactions, respectively. The charge transfer rate is related with the overpotential $\eta ,$ which is defined as:^26–28

$$\begin{eqnarray}&&\eta ={\phi }_{s}-{\phi }_{l}-{E}_{eq}-\displaystyle \frac{{\rm{\Omega }}{\sigma }_{h}}{F},\end{eqnarray} \tag{ 13 }$$

where ${E}_{eq}$ is the equilibrium potential for Li reaction in the active material. The hydrostatic stress ${\sigma }_{h}$ is expressed as $\left({\sigma }_{11}+{\sigma }_{22}+{\sigma }_{33}\right)/3.$ The last term in Eq. 13 includes the stress effect on the charge transfer reaction which accounts for the shift of the equilibrium potential by stresses.

Mechanical deformation is induced when Li-ions are inserted into or extracted from the active materials. The inhomogeneity of Li distribution and the mechanical interactions between the active particles and the conductive matrix generate stresses and cause subsequent mechanical degradation in the cell. The diffusion induced stresses are calculated using an analogy to thermal stresses.^48,49 The continuum media of the conductive matrix, containing carbon blacks, polymeric binders, and pores, is treated as a linear elastic material with an elastic modulus of 4 GPa.⁴⁶ In the domain of active particles, the total strain ε includes the elastic strain ε_e and Li-concentration induced strain ε_d:

$$\begin{eqnarray}&&{\boldsymbol{\varepsilon }}={{\boldsymbol{\varepsilon }}}_{e}+{{\boldsymbol{\varepsilon }}}_{d}.\end{eqnarray} \tag{ 14 }$$

The diffusion induced strain ${\varepsilon }_{d}$ can be written as:

$$\begin{eqnarray}&&{{\boldsymbol{\varepsilon }}}_{d}=\displaystyle \frac{1}{3}{\rm{\Omega }}\left({c}_{s}-{c}_{{s}_{0}}\right){\boldsymbol{I}},\end{eqnarray} \tag{ 15 }$$

where $\left({c}_{s}-{c}_{{s}_{0}}\right)$ is the concentration change of Li relative to the initial stress-free state. The deformation kinematics assumes an isotropic expansion in each direction upon Li reactions. This assumption is not correct to treat the deformation in a single crystal of layered cathode such as NMC. Nevertheless, the isotropic assumption is not entirely unreasonable to calculate the strain in the granular NMC secondary particles which are composed of agglomerates of primary particles. The constitutive equation describing the stress-strain relationship is given by:

$$\begin{eqnarray}&&{\boldsymbol{\sigma }}={\boldsymbol{C}}:{{\boldsymbol{\varepsilon }}}_{e},\end{eqnarray} \tag{ 16 }$$

where ${\boldsymbol{C}}$ is the stiffness matrix. For linear elastic materials, the stiffness matrix is reduced to two independent elastic constants—the elastic modulus and Poisson's ratio as listed in Table I. The stress field is governed by the equation of momentum equilibrium:

$$\begin{eqnarray}&&{\rm{\nabla }}\cdot {\boldsymbol{\sigma }}=0.\end{eqnarray} \tag{ 17 }$$

Table I.  Parameters used in the computational modeling.

Parameter	Value	Parameter	Value
${\kappa }_{s}$ (S·m−1)	104 50	${E}_{eq}$ (V)	4.43–2.78y+8.17y2 + 9.84 y3−214.51y4 + 777.77y5−1290.62y6 + 1034.38y7−324.15y8, y = ${c}_{s\_surf}$/${c}_{s\_\max }$
${\kappa }_{l}$ (S·m−1)	1.147 51
R (J·mol−1·K−1)	8.314	${D}_{l}$ (m2·s−1)	1 × 10−10 51
T (K)	293	${D}_{s}$ (m2·s−1)	7 × 10−15 52
F (A·s·mol−1)	96485	${{\rm{V}}}_{{\rm{total}}}$ ($\mu {{\rm{m}}}^{3}$)	8297
t+	0.363 51	E (GPa)	140 (NMC particles); 4 (Conductive matrix) 46
∂lnf /∂ln ${c}_{l}$	0.43 50
kc (m·s−1)	2 × 10−11 53	υ	0.3 46
ka (m·s−1)	2 × 10−11 53	Ω (m3·mol−1)	1.84 × 10−6
$\alpha$	0.5 53	${K}_{0}$ (N·m−2)	2 × 1015 44
${{\rm{SoC}}}_{\max }$	0.97 54	${k}_{d}$ (N·m−2)	0.375 × 1015
${{\rm{SoC}}}_{\min }$	0 54	${R}_{ref}$ (Ω· m2)	2 × 10−3 44
${c}_{s\_\max }$ (mol·m−3)	49000 54	${\lambda }_{ref}$ (nm)	10 44
${c}_{l\_ref}$ (mol·m−3)	1	${C}_{dl}$ (F·m2)	0.2 55
${c}_{l\_ini}$ (mol·m−3)	1000	${{\epsilon }}_{l}$ (1)	0.5

We take the approach of 3D representative volume element, Fig. 3, to examine the bias of mechanical stresses to the electrochemical responses in terms of the cell potential, surface overpotential, Li distribution, and cell capacity. The RVE model is based on the 3D X-ray phase contrast tomography data of NMC provided by Yang et al.⁸ as shown in Fig. 3a. The RVE maintains the primary morphological features of the original NMC cathode, which include the particle size polydispersity, particle spatial distribution, and inter-particle contact. We use the software Avizo to label the active NMC particles and then output the coordinate and size information of each particle. 1015 particles are distinguished in the tomography data in Fig. 3a. These data are set as the reference for selecting the representative volume element. Then, a RVE containing 22 particles (Fig. 3b) is selected and built as a Standard Tessellation Language (STL) data which can be imported into COMSOL. The RVE is incorporated into a half-cell configuration. Figure 3c shows the meshed geometry of the particles. We use statistical analysis to compare the RVE and the original NMC composite to ensure that the RVE is truly representative to the microscopic features of the cathode. Figures 3d, 3e, and 3f show the probability distribution functions of the particle diameter, the distance (surface to surface) of the nearest neighbors, and the average and standard deviation of the distance from the 1st, 2nd, and 3rd neighbors for individual particles, respectively.

Zoom In Zoom Out Reset image size

In the modeling, the empty space of the RVE is filled with a homogenized matrix which consists of the electrolyte and the carbon-binder matrix. Considering that the matrix is impermeable for Li-ions, only the electron transport is considered in the carbon-binder matrix. Furthermore, since the liquid electrolyte is insulating for electrons, only Li-ion transport is allowed in the electrolyte. The electrolyte and the conductive matrix can be treated as a superimposed continuum. In this case, the mass balance of Li-ions in the electrolyte component of the homogenized matrix (Eq. 1) should be modified as,

$$\begin{eqnarray}&&{{\epsilon }}_{l}\displaystyle \frac{\partial {c}_{l}}{\partial t}+{\rm{\nabla }}\cdot {{\boldsymbol{N}}}_{l}=0,\end{eqnarray} \tag{ 18 }$$

where ${{\epsilon }}_{l}$ is the volume fraction of electrolyte in the homogenized matrix (i.e. porosity of the matrix). In the NMC electrode in Fig. 3a and the reconstructed RVE in Fig. 3b, the volume fractions of active particles, the liquid electrolyte, and the carbon-binder matrix are about 40%, 30%, and 30%, respectively, resulting in the porosity of the homogenized matrix ${{\epsilon }}_{l}$ being 0.5. The pure phase transport properties including the matrix conductivity ${k}_{s},$ the electrolyte conductivity ${k}_{l},$ and the ionic diffusivity in the electrolyte ${D}_{l}$ should be replaced by the effective transport properties in the homogenized matrix. The Bruggeman relationship⁵⁶ is employed to evaluate the effective transport properties,

$$\begin{eqnarray}&&\begin{array}{l}{\kappa }_{s\_eff}={\left(1-{{\epsilon }}_{l}\right)}^{1.5}{\kappa }_{s},\\ {\kappa }_{l\_eff}={{{\epsilon }}_{l}}^{1.5}{\kappa }_{l},\\ {D}_{l\_{\rm{eff}}}={{{\epsilon }}_{l}}^{1.5}{D}_{l}.\end{array}\end{eqnarray} \tag{ 19 }$$

We incorporate the reconstructed RVE model into a half-cell configuration (Fig. 3c) for modeling the electrochemistry and mechanics fields. Due to the high electric conductivity of Li metal, the thickness of Li metal anode has a negligible impact on the cell voltage and polarization, and thus is set as zero. To simulate the galvanostatic charging and discharging, a constant current is applied at the cathode current collector side and its magnitude is determined by the specific cyclic rate. In the charging process, the following boundary conditions for ${{\boldsymbol{N}}}_{l},$ ${{\boldsymbol{N}}}_{s},$ ${{\boldsymbol{i}}}_{l}$ and ${{\boldsymbol{i}}}_{s}$ are prescribed at the interface between the electrolyte and active particles:

$$\begin{eqnarray}&&\begin{array}{l}{\bf{n}}\cdot {{\boldsymbol{N}}}_{l}=-\displaystyle \frac{{i}_{loc}}{F}\\ {\bf{n}}\cdot {{\boldsymbol{N}}}_{s}=\displaystyle \frac{{i}_{loc}}{F}\\ {\bf{n}}\cdot {{\boldsymbol{i}}}_{l}=-{i}_{loc}\\ {\bf{n}}\cdot {{\boldsymbol{i}}}_{s}={i}_{loc},\end{array}\end{eqnarray} \tag{ 20 }$$

where the unit normal vector ${\bf{n}}$ points from the electrolyte/matrix to the active particles.

The insulation and zero-flux boundary condition is applied at $x=0,\,x=L,\,z=0,$ and $z=H$ in Fig. 3c:

$$\begin{eqnarray}&&\begin{array}{l}{\bf{n}}\cdot {{\boldsymbol{N}}}_{l}=0\\ {\bf{n}}\cdot {{\boldsymbol{N}}}_{s}=0\\ {\bf{n}}\cdot {{\boldsymbol{i}}}_{l}=0\\ {\bf{n}}\cdot {{\boldsymbol{i}}}_{s}=0,\end{array}\end{eqnarray} \tag{ 21 }$$

where the unit normal vector ${\bf{n}}$ points outward of the model.

The current collector contacts with NMC electrode ($y=0$), the boundary conditions are prescribed by:

$$\begin{eqnarray}&&\begin{array}{l}{\bf{n}}\cdot {{\boldsymbol{N}}}_{l}=0\\ {\bf{n}}\cdot {{\boldsymbol{N}}}_{s}=0\\ {\bf{n}}\cdot {{\boldsymbol{i}}}_{l}=0\\ {\bf{n}}\cdot {{\boldsymbol{i}}}_{s}={i}_{app},\end{array}\end{eqnarray} \tag{ 22 }$$

where the unit normal vector ${\bf{n}}$ points from the electrode to the current collector. ${i}_{app}$ represents the applied current density when the cell is under galvanostatic charging/discharging. It corresponds to the macroscopic current density applied on the current collector in experiments. During charging, ${i}_{app}\gt 0.$ The value of ${i}_{app}$ is set by:

$$\begin{eqnarray}&&{i}_{app}=\displaystyle \frac{\left({{\rm{SoC}}}_{{\rm{\max }}}-{{\rm{SoC}}}_{{\rm{\min }}}\right){c}_{{s}_{{\rm{\max }}}}\cdot {{\rm{V}}}_{{\rm{total}}}\cdot {C}_{{\rm{rate}}}}{F},\end{eqnarray} \tag{ 23 }$$

where ${{\rm{SoC}}}_{\max }$ and ${{\rm{SoC}}}_{\min }$ are the maximum and minimum states of charge in the electrode, respectively, ${{\rm{V}}}_{{\rm{total}}}$ the electrode volume, and ${C}_{{\rm{rate}}}$ the C rate which represents that the cell takes $1/{C}_{{\rm{rate}}}$ hour to fully charge or discharge.

For the Li metal surface ($y={W}_{1}+{W}_{2}$), the boundary conditions are given by:

$$\begin{eqnarray}&&\begin{array}{l}{\bf{n}}\cdot {{\boldsymbol{i}}}_{l}={i}_{Li}\\ {\bf{n}}\cdot {{\boldsymbol{N}}}_{l}=\displaystyle \frac{{i}_{Li}}{F},\end{array}\end{eqnarray} \tag{ 24 }$$

where the unit normal vector ${\bf{n}}$ points from the Li metal to the cathode. ${i}_{Li}$ is the current density on the anode side, which is equal to ${i}_{app}$ due to the conservation of electrons on the cathode and anode sides.

When we calculate the mechanical stress, fixed displacement boundary conditions are applied at $x=0,\,x=L,\,\,z=0,z=H,\,\,y=0,$ and $y={w}_{1}+{w}_{2}$ in Fig. 3c.

The half-cell model is implemented into COMSOL Multiphysics (Version 5.3). The built-in time-dependent solver MUMPS (Multifrontal Massively Parallel sparse direct Solver) is used to solve the weak forms of the nonlinear and coupled electrochemistry and mechanics equations as outlined above. We first model a cycle of charging and discharge of the cell to explore the effects of stresses on the voltage profiles, overpotential, and equilibrium potential at the particle surfaces at two different charging rates. A positive current density ${i}_{app}$ is applied to the current collector during charging until the cell potential reaches 4.3 V. Then, the applied current density is reversed in sign upon discharging until the cell potential drops to 3.0 V. To make a quantitative assessment of the stress effect, we compare the simulation results of the fully coupled electro-chemo-mechanics model against the reference output without considering the stress effects on Li diffusion (Eq. 8) and the surface charge transfer (Eq. 13). The parameters used in the simulation are listed in the Table I.

Figure 4a shows the cell potential profiles under the charging rates 1 C and 5 C. The x-axis is plotted by a normalized time $t/\tau ,$ where $t$ is the Li reaction time and $\tau$ is the theoretical time for NMC to reach the fully charged/discharged capacity. Under the charging rates 1 C and 5 C, the values of $\tau$ are 3600 s and 720 s, respectively. The solid and dashed lines represent the results with and without considering the stress effects. It is clear that stress plays a significant role in modifying the cell potential response. In the case where the stress effect is not considered, the cell potential increases and decreases faster, leading to a shorter potential plateau and a lower yield capacity. This effect is magnified at a higher C-rate. To understand this phenomenon, we first look into the bias of mechanical stresses on the surface overpotential $\eta$ as defined in Eq. 13. We select a point at the surface of an NMC particle and plot the evolution of $\eta$ during charging in Fig. 4b. The overpotential remains a positive value, slowly decreases over time and drops to zero at the cutoff voltage. The stress at the particle surface modifies the overpotential, nevertheless, the difference caused by stresses (∼3 mV and ∼5 mV for 1 C and 5 C, respectively) is rather small compared to the magnitude of the electrochemical overpotential (∼40 mV and ∼80 mV for 1 C and 5 C, respectively). To isolate the stress effects on Li diffusion and the surface charge transfer, we turn off the stress regulation on Li diffusion (not considering the stress term in Eq. 8) but only consider the stress bias on the overpotential, and then plot the cell potential profiles in Fig. 4c for the two charging rates. Not surprisingly, the stress on the particle surface only changes the voltage profiles very minorly. Albeit small, it may be interesting to observe an opposite trend of the stress-regulation in the voltage profiles at the two charging rates. At the slow charging (black lines in Fig. 4c), the model including the stress effect reaches the cutoff voltage earlier and gives a relatively smaller capacity than the stress-free model. This trend is reversed at the high charging rate (red lines in Fig. 4c). When the NMC particles are subject to charging (Li extraction), a tensile stress is developed at the surface of the particles because of the depletion of Li. In the writing of the overpotential, $\eta ={\phi }_{s}-{\phi }_{l}-{E}_{eq}-\tfrac{{\rm{\Omega }}{\sigma }_{h}}{F},$ a tensile hydrostatic stress decreases the overpotential for the charge transfer reaction. For a slow charging, the redox reaction at the interface between the electrode and the electrolyte is relatively slow and therefore the Li source at the interface is sufficient. The deliverable capacity of the cell is essentially determined by the interfacial reaction. Because of the fixed applied current density, the overall reaction rates for the models with and without the stress effect are about the same. To maintain this overall reaction rate, the electric potential in the stress-coupled model would increase to compensate the decrease of the overpotential caused by stresses. Therefore, the cell potential in this scenario would reach the cutoff potential earlier. On the other side, when the cell undergoes a fast charging, Li on the particle surface is depleted faster. The retarding effect of stress on the surface charge transfer will in fact increase the Li concentration on the particle surface compared to the stress-free model. As a result, the equilibrium potential ${E}_{eq}$ decreases much slower during charging when stress is included. Figure 4d shows the profiles of ${E}_{eq}$ and $\tfrac{{\rm{\Omega }}{\sigma }_{h}}{F}$ at a selected location during charging. ${E}_{eq}$ is a monotonic function of the surface state of charge (SoC), which is defined as ${c}_{s\_surf}/{c}_{s\_\max }.$ The inset in Fig. 4d compares the surface SoC with and without including the stress effect. The comparison confirms that the surface SoC is higher at a given charging time and therefore ${E}_{eq}$ increases slower in the stress-coupled model. Figure 4d also shows the stress-biased overpotential item $\tfrac{{\rm{\Omega }}{\sigma }_{h}}{F}.$ Although it increases over time, the magnitude is much smaller than the ${E}_{eq}.$ In the fast charging case, we can conclude that the different electrochemical response is primarily due to the shift of the equilibrium potential. Again, by comparing the results in Figs. 4a and 4c, we find that the stress effect on charge transfer reaction is not the major mechanism regulating the cell potential response. The cell capacity and potential can be slightly tuned by the stress state at the particle surfaces depending on the subtle competition between the stress-biased ${E}_{eq}$ and $\tfrac{{\rm{\Omega }}{\sigma }_{h}}{F}$ at different charging rates, but the major contribution of mechanical stresses on the electrochemical response comes from the stress-regulated ion diffusion.

Zoom In Zoom Out Reset image size

To examine in detail the charge heterogeneity regulated by stresses and their mutual effects, we select an arbitrary NMC particle and plot the Li profile and the stress field during charging (Li extraction) under different C-rates in Fig. 5. The x-axis represents the normalized location x/R where x is the distance to the center of the particle and R is the radius of the particle. Here the radial direction is arbitrarily chosen within the particle. The solid and dashed lines of the same color show the results of Li profiles and the stress field with and without considering the stress effect on Li diffusion and surface charge transfer along the same radial direction in the same particle. The data are plotted at different states of charge of 3.5 V, 4.0 V, and 4.3 V. When the time scale of Li diffusion is longer than the charging time, both Li concentration and the stress field show a gradient distribution. In Figs. 5a (1 C rate) and 5c (5 C rate), Li concentration decreases over time and the region near the surface depletes Li faster. When the stress effect is included in the chemical potential of Li diffusion, Li distributes more homogenously. The stress effect is much more significant at the higher charging rate. The concentration gradient causes a mismatch of the lattice and therefore generates the intercalation stress. The equivalent (von-Mises) stress within the particle are shown in Figs. 5b and 5d. At the same state of charge, the magnitude of the stress in the stress-decoupled model (dashed lines) is always much higher than that in the stress-diffusion coupled model (solid lines). The stress and Li (charge) distribution mutually influence each other in the way that charge heterogeneity promotes the mechanical stresses, while the gradient of the stress field regulates the Li redistribution. This mutual effect is more significant when the cell is under fast charging and discharging.

Zoom In Zoom Out Reset image size

Now we look at the charge distribution and the stress field at the scale of the composite electrode. The morphology and the microstructural details, including the particle distribution, inter-particle contact, and the interactions between the active particle and the conductive matrix, determine the Li trajectories and the stresses. The cell level model is an ideal system to examine the charge heterogeneity and stress profiles relevant to the performance of batteries. Figure 6 shows the overall features of Li concentration and the equivalent stress after a cycle of charging and discharging of the cell under 1 C rate. Li distributions within individual particles are roughly spherically symmetric because the Li profile is dominated by the diffusion process. The deviation from the perfect-spherically symmetric solution is because of the stress effect induced by the mechanical interactions between the active particles and the matrix as well as the particle contact. The results derived from the models with and without considering the stress effect are compared in Figs. 6a and 6b. The cross-sectional views of various particles of different sizes are displayed in the dotted boxes—a small particle in the upper panel and a cluster of particles in the lower panel. In the Li distributions in Fig. 6a, we show that residual Li is retained in the active particles after a complete cycle which is a major reason of the decay in the deliverable capacity. Li retention in the large particles is more significant because of the longer diffusion length. In the stress-coupled model, the stress gradient in the particles drives Li diffusion. During the charging process, the surface regime of the particles experiences a tensile stress while the core regime is subject to a compressive stress. Such stress gradient facilitates the outward flow of Li. Likewise, the local stress alters the sign during discharging and therefore the stress gradient promotes the inward Li insertion. In both cases, the stress field facilitates Li transport and reduces the charge heterogeneity as seen in Fig. 6a. When the charge distribution is less heterogeneous, the internal stress is also minimized. Figure 6b plots the field of the equivalent stresses. We see that the stresses in the stress-coupled model are smaller, more homogeneous, and less likely to reach the extreme value. Figure 7 shows the results derived from the 5 C charging condition. When a larger electrochemical driving force (overpotential) is supplied, the charge heterogeneity is much more significant. The Li profile vastly deviates from the conventional core–shell distribution in the way that, first of all, the mechanical stresses generated by the lattice mismatch and the particle-matrix interactions result in different local current densities, and second, the charging protocol influences the Li distribution within individual particles—in a large particle for instance (red box in Fig. 7a), the core regime remains of high Li concentration after the completion of charging and discharging because of the Li retention at fast charging. The residual Li can be re-activated by the stresses. In the lower panel of Fig. 7a, we see that the compressive stress at the core reduces the Li gradient and therefore promotes the Li activity in the cell cycles. The reduction of the Li heterogeneity will in turn decrease the magnitude of the stresses. One point may be worth noting in Figs. 7b and 7d is that the local contacts of the inter-particles and the particle-matrix interface can experience drastic stresses at fast charging. The solid-solid contact domains marked in Figs. 7b and 7d are always subject to larger stresses which will most likely induce mechanical failure in those areas.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

To reduce the computational cost, we reconstruct a reduced model containing irregular-shaped NMC particles and similar microstructural features of the 3D X-ray phase contrast tomography data. Avizo is used to directly convert the original 3D tomography data into Standard Tessellation Language (STL). The open-source tool MeshLab is used to smooth the surface of NMC particles via Laplacian smoothing filter. The reduced model has a size of 25 μm × 25 μm × 36 μm as shown in Fig. 8a. The particle ensemble is then incorporated into the half-cell configuration. Figure 8b shows the meshed geometry of the particles in the finite element modeling.

Zoom In Zoom Out Reset image size

To model the particle detachment from the conductive matrix, we calculate the gap between the active particles and the matrix (denoted as $\lambda$) by prescribing a spring layer at the interface. The value of $\lambda$ increases when the NMC particles shrink in volume during Li extraction. The stiffness of the spring is denoted as ${K}_{int}.$ The exact value of ${K}_{int}$ and also the mechanical properties of the interface in general are unknown in the existing literature. However, a prior study⁵⁹ demonstrated a gradual loss of the interfacial reliability in the electrodes over cycles which infers that the contact stiffness reduces its value due to the repetitive size change of the active particles during electrochemical cycles. To serve the purpose of demonstration, we have taken a simplified phenomenological model of ${K}_{int}$ which is defined as a function of cyclic number:

$$\begin{eqnarray}&&{K}_{int}={K}_{0}-{k}_{d}\left({n}_{cyc}-1\right),\end{eqnarray} \tag{ 25 }$$

where ${K}_{0}$ is the initial stiffness. The coefficient ${k}_{d}$ represents the rate of reduction of the stiffness over cycles, and ${n}_{cyc}$ is the cyclic number. This simplified model can be refined or modified when experimental data becomes available. Here we take the stiffness ${K}_{int}$ in a range of $2\,\times \,{10}^{15}\,({\rm{N}}\,{{\rm{m}}}^{-2})$ in the pristine state of the electrode to $5\,\times \,{10}^{14}\,({\rm{N}}\,{{\rm{m}}}^{-2})$ at the completion of the cyclic modeling. The goal of this modeling is not to replicate the long-term cycling of batteries, but is to demonstrate how the mechanical damage at the interface deteriorates the cell performance. Here we have only simulated 5 completed cycles which is probably sufficient to draw a qualitative correlation between the mechanical aging and capacity degradation of batteries. At the formation of the interfacial gap, the interfacial resistance increases. We take a previously adopted model to describe the interfacial resistance ${R}_{gap}$⁴⁴:

$$\begin{eqnarray}&&{R}_{gap}={R}_{ref}\cdot \left({e}^{\tfrac{\lambda }{{\lambda }_{ref}}}-1\right),\end{eqnarray} \tag{ 26 }$$

where ${R}_{ref}$ and ${\lambda }_{ref}$ are phenomenological parameters of a reference interface resistance and a reference interfacial gap, respectively. Their values are list in the Table I. The interfacial resistance alters the driving force of surface charge transfer reaction. The overpotential drops by ${\rm{\Delta }}\eta ={R}_{gap}\cdot {i}_{loc}$ when the interfacial resistance evolves.

The electrochemical impedance spectroscopy (EIS) is a powerful tool for investigation of aging processes and determining the state of health of Li-ion batteries.^60–62 When a perturbation of potential of varying frequencies is applied on the electrode, the impedance response can be used to characterize the electrochemical properties. The response at the high frequency is related to the capacitance, electrochemical reactions, and local resistances, while the low-frequency response is contributed by the mass transport within the electrolyte and active particles.

To evaluate the regulation of interfacial failure on the cell performance, we conduct a series of EIS simulations to investigate the electrochemical properties after different number of cycles. The simulation is carried out by applying a perturbation of the potential of varying frequencies at the current collector in the half-cell model in Fig. 8. An electrochemical study is performed in COMSOL. The previous time-dependent variable is recast into the frequency domain, as defined as follows:

$$\begin{eqnarray}&&n={n}_{0}+Re\left\{\tilde{n}\cdot {e}^{2\pi f\cdot it}\right\},\end{eqnarray} \tag{ 27 }$$

where $n$ is the applied voltage, ${n}_{0}$ the initial value around which the perturbation takes place, $\tilde{n}$ the complex perturbation, $i$ the imaginary unit, $f$ the frequency in Hz, and t the time. The cell impedance, $Z({\rm{\Omega }}\cdot {m}^{2}),$ is calculated at the current collector of the cathode as follows:

$$\begin{eqnarray}&&Z=\displaystyle \frac{\tilde{{\phi }_{s}}}{{\boldsymbol{n}}\cdot \tilde{{{\boldsymbol{I}}}_{s}}},\end{eqnarray} \tag{ 28 }$$

where $\tilde{{\phi }_{s}}$ is the complex potential, ${\boldsymbol{n}}$ the normal vector at the boundary, and $\tilde{{{\boldsymbol{I}}}_{s}}$ the current density.

The impedance is modeled after the cell reaches the fully charged state 4.3 V. The potential of the current collector is set as a sinusoidal perturbation around 4.3 V with a 10 mV amplitude for the frequencies varying from 10 mHz to 1 kHz.⁵⁵ A double-layer capacitance ${C}_{dl}$ is applied at the interface between the active particles and the matrix.⁶³ The resistance ${R}_{gap}$ defined in Eq. 26 is implemented at the interface. Additional parameters used for the calculation of interfacial debonding and EIS analysis are listed in Table I.

To quantify the impact of the interfacial damage on the cell performance, we first determine the interfacial distance between the NMC particles and the conductive matrix over multiple cycles, and then conduct numerical modeling of EIS experiments at the charged states. Upon charging, NMC particles shrink in the volume and therefore the interfacial gap increases. When Li re-inserts in NMC during discharging, the active particles expand in size and the interfacial gap decreases. The comparison of the interfacial gap, averaged over all the active particle surfaces, are plotted in Fig. 9a under the 1 C and 5 C charging rates. The interfacial gap reaches the peak value at the completion of charging and reduces to a value close to zero at the end of a complete cycle. The magnitude of the interfacial gap depends on the charging rate. At a relatively slow charging 1 C, the interfacial gap reaches a peak value over 60 nm. Overall the interfacial gap at 1 C is larger than that in the fast charging 5 C case. This is because a larger number of Li are inserted or extracted from the NMC particles at the slow charging and therefore the mismatch strain at the interface is more significant. Another important feature is that the change of the interfacial gap, especially under fast charging, is not reversible. This is because of the asymmetry of the charging and discharging capacities in the cell within a complete cycle. The asymmetry is induced by the charge redistribution modulated by the mechanical damage. We will discuss this aspect in more detail in next figure. The interfacial gap evolves over cycles and will eventually lead to the fade of the capacity. Here we use the Nyquist plot, a complex plot of the electrochemical impedance for a cell, to evaluate the transport properties for the short-term or long-term processes. Figure 9b shows the Nyquist plots (solid lines) simulated with and without considering the interfacial debonding in the model. All the Nyquist plots show a typical semicircle at the high frequency domain which is attributed to the charge transfer reactions, and a Warburg tail at low frequency which is due to the Li diffusion in the active particles. The transport properties for different electrochemical processes can be extracted from the Nyquist plots using an equivalent circuit model. The dashed lines represent the fitting curves using a circuit model shown in the inset of Fig. 9b, where R_Ω represents the electric ohmic resistance from the electrolyte, the carbon-binder matrix, and the interface between the active and inactive materials, R_ct the resistance for the surface charge transfer, C_dl the capacitance for the electric double layer, and Z_w the Warburg impedance for Li diffusion in the bulk active material. Table II summarize the fitted impedance parameters for the models with and without considering the interfacial debonding. Both R_Ω and R_ct significantly increase when the interfacial debonding occurs over cycles. The loss of contact of the active materials to the conductive agents impedes electron transport and the charge transfer reaction. Ultimately, this interfacial degradation results in a high polarization and significant capacity loss in batteries.

Zoom In Zoom Out Reset image size

Table II.  Comparison of fitted EIS parameters with and without considering interfacial debonding at the fully charged states (4.3 V) in different cycles.

	${R}_{{\rm{\Omega }}}({\rm{\Omega }}\cdot {{\rm{m}}}^{2})$	${R}_{ct}({\rm{\Omega }}\cdot {{\rm{m}}}^{2})$	${C}_{dl}({\rm{F}}\,{{\rm{m}}}^{-2})$	${Z}_{w}({\rm{\Omega }}\cdot {{\rm{m}}}^{2}\,{{\rm{s}}}^{-1/2})$
No debonding	0.00039	0.01515	1.69	0.00230
Debonding/1st charging	0.02471	0.02172	0.2036	0.00274
2nd charging	0.03985	0.02361	0.1871	0.00292
3rd charging	0.07099	0.02631	0.1692	0.00321
4th charging	0.14360	0.03019	0.1507	0.00370
5th charging	0.34520	0.03554	0.1321	0.00454

Figure 10 shows the evidence of the charge redistribution and Li retention over cycles modulated by the mechanical damage at the interface. The interfacial debonding within the composite is heterogeneous—the degree of damage depends on the local Li concentration (volumetric strain) and the stress state. Figures 10a and 10b demonstrate the Li profile in an arbitrarily selected particle at the end of each cycle. A comparison is made for the models with and without considering interfacial debonding. The particle detachment from the conductive matrix modifies the local resistance of charge transfer and therefore induces a dynamic nature of the charge/Li distribution. When the interfacial damage occurs, Li is retained within the particles at each charged state and further accumulates while the electrochemical cycle proceeds. The residual Li becomes inactive in the electrochemical process which ultimately results in a gradual downgrading of the deliverable capacity of the cell. The impact of the interfacial damage on the cell polarization, voltage response, and capacity retention is detailed in Fig. 11. The cell potential over the first few cycles is plotted in Fig. 11a. The potential responses at two charging rates 1 C and 5 C are compared. At the fast charging rate, the cell experiences more polarization and delivers less Li capacity compared to the slow charging condition. The black and red lines show the results with and without considering the interfacial debonding, respectively. When the interfacial contact between the active particles and the conductive matrix is maintained, the cyclic voltammetry reaches a steady state at the 2nd cycle due to the incomplete Li extraction/insertion in the 1st cycle. When the interfacial failure is included, the interval between each charging and discharging decreases as the cycles proceed because of the Li retention within the active particles as discussed above. The reduction of Li activity caused by the interfacial delamination results in a steady increase of the polarization and minimization of the deliverable capacity of the cell. Such effect is magnified at the high charging rate. Figure 11b plots the fade of cell capacity over cycles. The capacity ratio is calculated by the actual capacity in each charging and discharging divided by the theoretical value. The cell capacity in general drops after the first cycle because of the diffusion limited Li insertion/extraction. Consistent with the voltage profile, the cell capacity maintains a nearly constant in the following cycles when the active particles remain well contacted with the conductive agent. However, when the interfacial stability is considered in the picture, the mechanical failure becomes a determining factor of the capacity fade in battery cycles. Note that in the modeling of the cyclic behavior, we have used the fully coupled theory on stress, diffusion, and charge transfer. The result, when the interfacial debonding is not considered, shows that the stress regulation on Li/charge distribution is reversible in the cycles. It is the mechanical damage at the interface making the charge distribution, capacity, and voltage response irreversible in the cyclic process.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size