3D Heterogeneous Model for Electrodes in Lithium-Ion Batteries to Study Interfacial Detachment of Active Material Particles and Carbon-Binder Domain

Abstract

Mechanics plays a crucial role in the performance and lifespan of lithium-ion battery (LIB) cells. Thus, it is important to address the interplay between electrochemistry and mechanics in LIBs, especially when aiming to enhance the energy density of electrodes. Accordingly, this work introduces a framework for a fully coupled electro-chemo-mechanical heterogeneous 3D model that allows resolving the inhomogeneities accompanied by electrochemical and mechanical responses of LIB electrodes during operation. The model is employed to numerically study the mechanical degradation of a nickel manganese cobalt (NMC) cathode electrode, assembled in a half-cell, upon cycling. As opposed to previous works, a virtual morphology for a high-energy electrode with low porosity is developed in this study, which comprises distinct domains of active material (AM) particles, the carbon-binder domain (CBD), and the pore domain to resemble real commercial electrodes. It is observed that the mechanical strain mismatch between irregularly and randomly positioned AM particles and the CBD might lead to local contact detachment. This interfacial gap, in combination with the diminishing contact strength over cell cycling, continuously deteriorates the electrode performance upon cycling by impedance rise and capacity drop. In agreement with previous experimental reports, the presented simulation results exhibit that the contact loss mostly takes place in the regions closer to the separator. Eventually, the resulting gradual capacity drop and change in impedance spectrum over cycling, as the consequence of interfacial gap formation, are discussed and indicated.

1. Introduction

Lithium-ion batteries have established themselves as the predominant energy storage system due to their superior energy and power densities, as well as their extended operational lifespan. Therefore, they have become the prevailing choice for energy storage across a diverse spectrum of applications. These applications encompass electric vehicles (EVs), consumer electronics, energy storage systems, and various eco-friendly industries [1]. Nonetheless, given the ever-growing demand within industries for heightened energy density, enhanced performance, and cost-effective solutions, ongoing optimization investigations remain a vital necessity in the pursuit of advancing LIB technology [2].

Given the substantial time and resources required for the development and characterization of new lithium-ion cell types in both research laboratories and commercial cell manufacturing facilities, computational modeling and simulation emerge as efficient tools, offering fast and in-depth insights at a microscopic level into the intricate processes occurring within lithium-ion electrodes during charging and discharging. In this context, the continuum model, a widely employed physics-based framework [3,4], simplifies the complex electrode microstructure by treating it as a continuum through the application of approximation methods. This simplification substantially reduces computational costs. However, it is essential to acknowledge that the existing insufficiencies inherent in the continuum model limit its capability in a range of applications [5,6,7].

The microstructural architecture and spatial arrangement of domains constituting Lithium-ion battery electrodes exert profound influence over the complicated interplay between diverse electrochemical and mechanical processes during LIB operation [8]. Consequently, to attain a comprehensive understanding of the effect of electrode structure on LIB performance, it becomes imperative to construct a virtual representation that emulates the actual, heterogeneous electrode structure. Creating a model from tomography imaging is an expensive and time-intensive endeavor, entailing the inevitable destruction of samples, thereby imposing limitations on iterative optimization and model enhancement, especially in the context of prototype electrodes or cells [9,10,11]. Additionally, the data derived from imaging techniques might be insufficient for comprehending battery electrochemistry and mechanics due to various constraints. In contrast, the computational modeling approach offers a conceptual design in the initial phases of research projects, facilitating optimization with respect to desired functionalities. In other words, virtual material testing and experimentation with diverse electrode properties can be conducted efficiently and cost-effectively [10].

Developing a heterogeneous model provides deep insights into localized non-uniformities within the electrode structure, leading to uneven utilization of active materials. Such inhomogeneity may result in degradation phenomena, like lithium plating on anode material particles [12,13], particularly those near the separator interface during fast charging of electric vehicles. Additionally, this model enables investigations into optimized battery operating parameters and contributes to the advancement of microstructural engineering strategies [5,14,15,16] aimed at mitigating aging phenomena within the electrode.

Mechanics play a pivotal role in determining both the performance and longevity of lithium-ion batteries. With the growing demand for extended cycle life, fast charging, and increased driving range in EV applications, mechanical degradation is one of the obstacles that directly regulates the mechanisms of capacity deterioration [17,18]. Therefore, the mutual impacts of electrochemistry and mechanics in LIB cells need to be addressed in the pursuit of electrodes with high energy density. In this regard, the study of electrochemical–mechanical interactions within an electrode structure can be effectively conducted through the application of heterogeneous models. This encompasses the examination of the impact of lithiation and delithiation on the evolution of mechanical strain and stress within active material particles [19,20], the influence of the mechanics on electrochemical response [21,22,23,24], and the exploration of potential fracture nucleation and propagation phenomena within the electrode structure [25,26,27]. Accordingly, the model enables the development of measures to suppress mechanical degradation, ranging from the electrode’s geometrical features to mechanical properties to operational considerations such as cycling voltage windows.

In addition to the fracture inside active material as one degradation mechanism, the mismatch between the mechanical behaviors of AM particles and inactive materials may cause decohesion at the interfaces between them. The conductive additive, which is typically carbon black, provides the electronic conduction pathways between the active particles and the current collector (CC). The contact between the AM particle and the conductive additive could be lost during cycling, partially hindering the electronic path. As a result, the impedance of the electrode increases, and the capacity retention is reduced, leading to electrode performance degradation. In contrast to the particle fracture, the investigation of mechanical failure at the interface has received comparatively less attention. This disparity arises from the difficulty in experimentally resolving the intricate interfacial interactions between the active materials and the conductive agents. Nevertheless, the scanning electron microscope (SEM) images presented by some researchers confirmed that the continuing deformation of AM particles during cycling can result in decohesion at the interface between the active material and the conductive agents [28]. Figure 1 shows cross-sectional SEM images of an NMC cathode electrode before and after 10 cycles at 5 C-rate to visualize the interfacial debonding. While Figure 1a illustrates the existing cohesion between the particle and surrounding inactive matrix, the incurred interfacial decohesion in the cycled state is highlighted by the red dashed lines in Figure 1b. In the rapidly advancing field of solid-state batteries, the occurrence of delamination between particles and the solid electrolyte (SE) can be more pronounced compared to conventional batteries employing liquid electrolytes. Figure 2 contains SEM images showing the developing bond at the interface between NMC particles and solid electrolyte in three states: the pristine state, after the first cycle, and after 50 cycles. The SEM images distinctly reveal the NMC particles due to their characteristic, nearly spherical morphology. While the NMC particles and the solid electrolyte form a densely packed composite electrode in the pristine state, (Figure 2a,b), the NMC particles are surrounded by gaps after the first charging or after 50 cycles, (Figure 2c–f).

Figure 1. Cross-sectional SEM images of NMC particles and surrounding inactive matrix in (a) pristine state and (b) cycled state after 10 cycles at 5 C-rate of a lithium-ion cell cathode reported by Xu et al. [28] (shown with publisher’s permission).

Figure 1. Cross-sectional SEM images of NMC particles and surrounding inactive matrix in (a) pristine state and (b) cycled state after 10 cycles at 5 C-rate of a lithium-ion cell cathode reported by Xu et al. [28] (shown with publisher’s permission).

The mechanistic comprehension of LIB cell degradation faces a critical challenge due to the intrinsic heterogeneity in the structural configuration and the composition of the electrode. While Yang et al. [29] have demonstrated a noteworthy spatial and temporal dependency of damage in a commercial electrode, many prior experimental studies primarily concentrated on an idealized structure, such as thin films and individual particles. Others focused on localized regions within the electrode, characterized by limited dimensions, thereby being unable to provide a statistical representation of the overall chemomechanical behavior exhibited by the electrodes at a large scale. Hence, the exploration of heterogeneous chemomechanics and damage within lithium-ion batteries remains at its early stages of development, with a multitude of questions yet to be elucidated. The above-mentioned experimental challenges also exist in numerical investigations. The common continuum modeling approach and the extended single-particle approach oversimplify and assume free-standing spherical particles of the same size that are not subject to mechanical or electrochemical constraints from the adjacent environment. Recently, Baboo et al. [30] employed a simple single-particle model (SPM) and adjusted the solid diffusivity and specific active surface area to predict the effects of the formation of solid electrolyte interphase (SEI), side reactions, and AM particle fragmentation on the degradation of a LIB full cell performance during cycling. Their investigation was conducted for three different binders with diverse adhesion strengths. However, the simple SPM was unable to replicate the decreasing capacity during cycling, which was observed in the experiment. It can be attributed to the existing simplifications in SPM to express the complicated and spatially-resolved physicochemical interactions taking place inside the LIB cell during operation. Nonetheless, in commercial batteries, both the cathode and anode exhibit a composite nature characterized by substantial heterogeneity on the nano- to microscale, comprising irregularly shaped and sized active particles within a matrix of polymeric binders and conductive additives, along with pores containing electrolyte.

Figure 2. SEM images of the cathode composite of NMC particles embedded in the solid electrolyte (${\mathrm{Li}}_3{\mathrm{PS}}_4$) (a,b) as prepared in pristine state, (c,d) after single charging at 0.1 C-rate and, (e,f) after 50 full cycles in the discharged state provided, by Koerver et al. [31] (shown with publisher’s permission).

Figure 2. SEM images of the cathode composite of NMC particles embedded in the solid electrolyte (${\mathrm{Li}}_3{\mathrm{PS}}_4$) (a,b) as prepared in pristine state, (c,d) after single charging at 0.1 C-rate and, (e,f) after 50 full cycles in the discharged state provided, by Koerver et al. [31] (shown with publisher’s permission).

Recently, a few researchers have developed experimental setups and fully coupled electro-chemo-mechanical heterogeneous models to investigate the interface of AM particles and inactive matrix in electrodes. Xu et al. [28] employed experimental tomography to visualize nanoscale interfacial debonding between NMC particles and the conductive matrix. Their observations revealed that the interfacial detachment during cycling primarily occurs in the electrode region near the separator, rather than in the vicinity of the current collector. Moreover, Xu et al. built a microstructural-resolved model utilizing tomographic data to examine the heterogeneous damage within composite cathode particles situated within the binder matrix. However, the volumetric fraction of AM particles in the reconstructed virtual morphology by Xu et al. equates ${\epsilon }_{\mathrm{AM}}^{\mathrm{ave}}=40\%$ and is thus not as dense as the current commercial electrodes with high energy densities. Moreover, the pore and carbon-binder domains in the re-built morphology in their work are simplified and represented by a single integrated composite domain. Later, Liu et al. [32] also investigated the degradation of contacts between the AM particles and the inactive matrix in an NMC half-cell electrode using a reconstructed morphology with similar simplified geometrical features as developed in the work of Xu et al. [28]. Accordingly, Liu et al. adopted a phenomenological approach to reproduce the mechanical fatigue of cohesion along the AM particles and the inactive matrix during cycling. It was observed that the steady decay of the interfacial strength causes a growing interfacial debonding at cycling.

In this work, a 3D microstructural-resolved model is developed to study the impact of degradation of contact strength between the AM particles and CBD in the electrode in LIB cells during cycling. In contrast to previous models, the constructed virtual morphology of the electrode encompasses distinguished domains to represent the individual existing components inside the real electrode. Such a heterogeneous model enables predicting the electrochemical and mechanical behavior of the electrode’s constituents more reliably. In the end, the influence of interfacial decohesion upon cycling is studied and discussed in terms of capacity retention and electrochemical impedance spectroscopy (EIS) in the time and frequency domains, respectively.

2. Methods

2.1. Electrode Microstructure

In this work, the heterogeneous model, introduced by the authors in their previous work [2], is employed to build a fully coupled electro-chemo-mechanical model and study the interfacial bonding of NMC622 AM particles and carbon-binder domain in a cathode electrode with conventional liquid electrolyte. The presented model allows to observe the heterogeneous electrochemistry, stress, and interfacial detachment in a commercial cathode electrode. Therefore, the spatial and temporal variations on different levels can be investigated. The offered modeling approach in this work can be further utilized to study the mechanical degradation in solid-state batteries as well.

As opposed to the developed cathode electrode’s virtual morphology with low AM domain volumetric fraction by Xu et al. [28] and Liu et al. [32], a denser cathode electrode with ${\epsilon }_{\mathrm{AM}}^{\mathrm{ave}}=62.83\%$ has been built in this study to represent an electrode with higher energy density. In addition, the generated morphology contains separate distinguished domains of pore and CBD that increase the reliability of simulation results as compared to the assumed simplification made by the previous research to consider these two domains as one integrated domain.

Subsequently, a half-cell representative volume element (RVE) is assembled by incorporating the developed cathode electrode microstructure and two cubic domains. The resulting model explicitly encompasses the crucial half-cell components, including the separator, current collector, and lithium foil, as vividly illustrated in Figure 3.

Table 1. Half-cell RVE microstructure generation specifications.

Property	Value
Cathode microstructure dimension	$50\times 50\times$ 25 μm3
Average NMC particle radius $r_p$	2 μm and 5 μm
Particle radius standard deviation	0.2 μm
Active material volumetric fraction (${\epsilon }_{\mathrm{AM}}^{\mathrm{ave}}$)	62.83%
Carbon-binder domain volumetric fraction (${\epsilon }_{\mathrm{CBD}}^{\mathrm{ave}}$)	10.55%
Porosity (${\epsilon }_{\mathrm{e}}^{\mathrm{ave}}$)	26.6%
Separator thickness	11 μm
Al current collector thickness	6.5 μm

includes the considered geometrical properties in the generation of the half-cell RVE.

The simulations are carried out in COMSOL Multiphysics^® version 6.0 using an Intel^® Xeon^® CPU @ 3.70 GHz (2 processors) with 128 GB of RAM. Moreover, COMSOL LiveLink^® for MATLAB^® is used in the current work to transfer model’s mesh data from MATLAB ^® version R2020b to COMSOL Multiphysics^® software to create the FEM model.

2.2. Electro-Chemo-Mechanical Model

The governing equations to describe diverse coupled electrochemical and mechanical interactions within various domains, along with their corresponding parameters, are delineated in

Table 2. List of physics implemented in different domains of the half-cell.

Domain/Boundary	Equation	No.
Electrochemistry
AM, CBD, CC	$\nabla ·\left({\sigma }^{\mathrm{eff}}{*}_s{\nabla }_{{\varphi }_s}\right)-j^{{\mathrm{total}}^{**}}=0$	$\left(1\right)$
AM	$\frac{\partial c_s}{\partial t}=\nabla .(D_s\nabla c_s)-(\frac{D_s\Omega }{RT}\nabla c·\nabla {\sigma }_h+c\nabla ·\nabla {\sigma }_h)$	$\left(2\right)$
Pore, CBD, separator	$\frac{\partial \left({\epsilon }_e^{***}{c}_e\right)}{\partial t}=\nabla ·D_e^{\mathrm{eff}}\nabla c_e+\frac{1-t_{+}^0}{F}{j}^{{\mathrm{total}}^{**}}$	$\left(3\right)$
	$\nabla .({\kappa }^{\mathrm{eff}}\nabla {\varphi }_e)+\nabla ·({\kappa }_D^{\mathrm{eff}}\nabla ln\left(c_e\right))+j^{{\mathrm{total}}^{**}}=0$	$\left(4\right)$
CBD, Separator	$D_e^{\mathrm{eff}}={\left({\epsilon }_e\right)}^p·D_e^{\mathrm{bulk}}$	$\left(5\right)$
	${\kappa }_e^{\mathrm{eff}}={\left({\epsilon }_e\right)}^p·{\kappa }_e^{\mathrm{bulk}}$	$\left(6\right)$
CBD	${\sigma }_s^{\mathrm{eff}}={\left({\epsilon }_s\right)}^p·{\sigma }_s^{\mathrm{bulk}}$	$\left(7\right)$
AM-electrolyte interface	$j^f=A_v·j_0\left(exp\left(\frac{{\alpha }_{a}F\eta }{RT}\right)-exp\left(-\frac{{\alpha }_{c}F\eta }{RT}\right)\right)$	$\left(8\right)$
	$j_0=Fk{c}_s^{{\alpha }_c}{\left(c_{s,\max }-c_s\right)}^{{\alpha }_a}{\left(\frac{c_e}{c_{e,\mathrm{ref}}}\right)}^{{\alpha }_a}$	$\left(9\right)$
	$\eta ={\varphi }_s-{\varphi }_e-U-\frac{\Omega {\sigma }_h}{F}$	$\left(10\right)$
	$j^c=A_v·\frac{\partial ({\varphi }_s-{\varphi }_e)}{\partial t}·C_{\mathrm{DL}}$	$\left(11\right)$
	$j^{\mathrm{total}}=j^f+j^c$	$\left(12\right)$
Mechanics
All domains	$\nabla ·\sigma =0$	$\left(13\right)$
	$\sigma =\mathbf{C}:{ϵ}_{\mathbf{e}}$	$\left(14\right)$
	$ϵ=\frac{1}{2}\left[\nabla \mathbf{u}+{(\nabla \mathbf{u})}^T\right]$	$\left(15\right)$
	$ϵ={ϵ}_{\mathrm{e}}+{ϵ}_{\mathrm{Li}}^{****}$	$\left(16\right)$
AM	${ϵ}_{\mathrm{Li}}=\frac{\Omega }{3}(c_s-c_{s,0})\mathbf{I}$	$\left(17\right)$
AM-CBD interface	$K_{\mathrm{int}}=K_{int,0}-k_d\left(n_{cyc}-1\right)$	$\left(18\right)$

* ${\sigma }_s^{\mathrm{eff}}$ for CBD domain while ${\sigma }_s^{\mathrm{bulk}}$ in AM and CC domains. ** $j^{\mathrm{total}}=0$ in CC and separator domains. *** ${\epsilon }_e$ excluded in equation for pore domain. **** ${ϵ}_{\mathrm{Li}}=0$ for all domains except AM.

, Table 3 and Table 4, respectively. In the equations and associated parameters, the subscripts e and s denote the electrolyte and the solid-phase properties, respectively.

2.2.1. Solid Domains

Ohm’s law, Equation (1) in Table 2, is solved in the AM, CBD, and CC domains to describe electronic migration. Current collector domains are impermeable to lithium ions, so normal ion flux is set to zero at their boundaries.

At the interfaces between AM and the electrolyte, the Butler–Volmer (BV) equation is used to describe the electrochemical reaction kinetics or the volumetric faradaic current density $j^f$. In the BV equation, the total overpotential is modified by introducing a stress-biased term to account for the effect of mechanics on the kinetics of the faradic reaction. Accordingly, the fourth term on the right side of Equation (10) in Table 2 represents the stress-biased term when $\Omega$ and ${\sigma }_h$ denote the partial molar volume of NMC particles and the hydrostatic stress, respectively. In addition, to account for the formation of the double layer (DL) at this interface, a homogeneous volumetric capacitive current density represented as $j^c$ is incorporated alongside the faradaic current density. These components collectively contribute to the determination of the total volumetric current density, denoted as $j^{\mathrm{total}}$, at these interfaces.

To describe solid-state diffusion, the formulation of Fick’s second law is modified to account for the stress-induced diffusion flux, Equation (2) in Table 2. The additional terms inside the second bracket on the right side of this equation reflect the stress-driven flux of the lithium inside the AM particles. Then, at the interface between the active material and the electrolyte, the Neumann boundary condition is established to consider the interplay between the faradaic current density and the diffusion processes occurring within the AM particles.

2.2.2. Pore and Separator Domains

It is assumed that the pore domain located between the solid domains is completely saturated with the liquid electrolyte. The separator domain is treated as an electrolyte-permeable media with a porosity of 50%. To account for the separator’s porous nature, the Bruggeman relation is utilized to predict the effective transport coefficients of the separator domain.

Also, a rather small porosity of 27.6% is assigned to the CBD to account for the reported porosity factor by Daemi et al. [33]. Accordingly, the electrolyte’s effective ionic conductivity and diffusivity in CBD are obtained by the corresponding bulk values scaled by a factor of 0.276, following Boyce et al. [34].

The concentrated solution theory is employed to mathematically describe electrolyte transport. A porosity factor ${\epsilon }_e$ is introduced to Equations (3) and (4) in Table 2 within domains of the separator and the CBD to consider their porous nature, which is not explicitly included in their re-built modeled morphology.

2.2.3. Lithium Foil Boundary

The lithium foil is modeled as a boundary, and the BV equation is used to compute the volumetric current density associated with lithium metal deposition and dissolution at the interface between the lithium foil and the electrolyte. Here, the Butler–Volmer equation is formulated with zero overpotential in order to not introduce additional impedance. As a result, the calculated impedances are primarily attributed to contributions originating from the porous separator and cathode microstructure, saturated with electrolyte, and a current collector foil attached.

2.2.4. Mechanics

A mechanical analysis is performed to consider the effects of lithiation-induced deformation on both mechanical and electrochemical responses of the various electrode constituent domains. The heterogeneity of the reaction kinetics and the non-uniform mechanical behavior of the active and inactive materials lead to the development of stresses followed by mechanical degradation in the cell. Due to the rather slow kinetics of solid diffusion within AM particles, the mechanical equilibrium equation is solved in the absence of body forces, Equation (13) in Table 2, where $\sigma$ denotes the Cauchy stress tensor. In the case of the cathodic NMC electrode with limited expansion upon lithiation, Hooke’s law, Equation (14), is employed as the constitutive relation to describe the linear elastic material model. Accordingly, $\mathbf{C}$ and ${ϵ}_{\mathbf{e}}$ refer to the elastic stiffness tensor and the elastic strain, respectively. By assuming small strains and rigid body rotations, the Green-Lagrange strain tensor is reduced by the removal of quadratic terms, resulting in the well-known compatibility relation, Equation (15), which correlates the total strain tensor and the symmetric part of the gradient of the displacement field. The total strain is further decomposed into the elastic and lithiation-induced parts, as expressed in Equation (16). The lithiation-induced strain in the AM particles, ${ϵ}_{\mathbf{Li}}$, is attained by Equation (17), in analogy to the thermal strain in materials, where $\mathbf{I}$ is the identity tensor. An initial active material concentration, $c_{s,0}$, is prescribed whilst the electrode and all constituent domains are assumed to be initially in an unstressed state. Equation (17) formulates the contribution of lithiation-induced deformation by the development of hydrostatic stress. As illustrated in Figure 3, while symmetrical boundary condition is applied on the RVE sides in the x-y and x-z planes, fully-clamped and free boundary conditions are allocated for the remaining two sides in the y-z plane at $x=l_{electrode}$ and $x=0$, respectively.

In order to evaluate the strength along the interfaces between the AM particles and the carbon-binder domain, an interfacial cohesion model is required. Despite the model developed by Müller et al. [35], which considers a serial spring and a dissipative damper to reproduce the mechanical response of the contact between the AM particles and carbon-binder domain, the induced mechanical hysteresis effect is ignored in this work. Accordingly, a spring layer with zero thickness across the electrode represents the interfacial contact between the AM cathode particles and the CBD. The constant value of this spring per unit area in the pristine state is denoted as $K_{int,0}$. When the spring expands along the interface under tensile stress, the contact locally detaches, hindering electron transfer paths and elevated impedance. On the other hand, the bond is maintained for the remaining contacts, which are locally under compression or no stress. The above-mentioned spatial contact loss induces a drop in the local solid potential ${\varphi }_s$ along the debonded contact, which in turn has a local impact on the charge conservation equation in the solid domain. Moreover, due to the introduced drop in solid potential by detachment, the overall overpotential $\eta$, as the driving force for electrochemical reaction, is retarded, and consequently, the local reaction kinetics described via the Butler–Volmer equation is impeded, resulting in degradation of battery cell performance.

Zhu et al. [36] demonstrated that the interfacial strength gradually decreases due to mechanical fatigue during cycling. Similar to the phenomenological approach presented by Liu et al. [32], the fatigue deterioration of contact between AM particles and CBD during cycling is realized in this work by a decreased spring constant along this interface. Equation (18) expresses the considered degradation of the spring stiffness per cycle, where $k_d$ represents the rate at which the spring stiffness decreases per cycle, and $n_{cyc}$ indicates the cycle number.

Table 3. Electrochemical parametrization of the half-cell model.

Parameters	Value	Ref.
NMC622 Particles
AM solid conductivity (${\sigma }_s$)	$1.6\times {10}^{-4}$ $\mathrm{S}·{\mathrm{m}}^{-1}$	[34]
AM solid diffusivity ($D_s$)	$f\left(SoL\right)$ ${\mathrm{m}}^2·{\mathrm{s}}^{-1}$	[34]
max AM solid concentration ($c_{s,\max }$)	48,700 mol · m−3	[34]
Initial AM solid concentration ($c_{s,0}$)	500 mol · m−3
Equilibrium potential (U)	$\begin{array}{l}-324.2·So{L}^8+1034.4·So{L}^7-129.6·So{L}^6-\\ 777.8·So{L}^5-214.5·So{L}^4-9.8·So{L}^3+\\ 8.2·So{L}^2-2.8·SoL+4.4 \mathrm{V}\end{array}$	[34]
Kinetics
Reaction rate constant (k)	$2\times {10}^{-11}$ $\mathrm{m}·{\mathrm{s}}^{-1}$	[34]
Transfer coefficients (${\alpha }_a$, ${\alpha }_c$)	$0.5$
Surface double layer capacitance ($C_{\mathrm{DL}}$)	$0.2$ $\mathrm{F}·{\mathrm{m}}^{-2}$	[28]
Bruggeman exponent (p)	$1.0$	[34]
Current collector
Conductivity (${\sigma }_s$)	$3.7\times {10}^7$ $\mathrm{S}·{\mathrm{m}}^{-1}$	[34]
Carbon-binder domain
Conductivity (${\sigma }_s$)	375 $\mathrm{S}·{\mathrm{m}}^{-1}$	[34]
Electrolyte
Initial electrolyte concentration ($c_{e,0}$)	1000 $\mathrm{mol}·{\mathrm{m}}^{-3}$	[34]
Conductivity (${\kappa }_e$)	$f\left(c_e\right)$ $\mathrm{S}·{\mathrm{m}}^{-1}$	[37]
Diffusivity ($D_e$)	$f\left(c_e\right)$ ${\mathrm{m}}^2·{\mathrm{s}}^{-1}$	[37]
Activity ($\partial lnf/\partial ln{c}_e$)	$0.43$	[37]
Transference ($t_{+}^0$)	$0.37$	[37]

Table 3. Electrochemical parametrization of the half-cell model.

Parameters	Value	Ref.
NMC622 Particles
AM solid conductivity (${\sigma }_s$)	$1.6\times {10}^{-4}$ $\mathrm{S}·{\mathrm{m}}^{-1}$	[34]
AM solid diffusivity ($D_s$)	$f\left(SoL\right)$ ${\mathrm{m}}^2·{\mathrm{s}}^{-1}$	[34]
max AM solid concentration ($c_{s,\max }$)	48,700 mol · m−3	[34]
Initial AM solid concentration ($c_{s,0}$)	500 mol · m−3
Equilibrium potential (U)	$\begin{array}{l}-324.2·So{L}^8+1034.4·So{L}^7-129.6·So{L}^6-\\ 777.8·So{L}^5-214.5·So{L}^4-9.8·So{L}^3+\\ 8.2·So{L}^2-2.8·SoL+4.4 \mathrm{V}\end{array}$	[34]
Kinetics
Reaction rate constant (k)	$2\times {10}^{-11}$ $\mathrm{m}·{\mathrm{s}}^{-1}$	[34]
Transfer coefficients (${\alpha }_a$, ${\alpha }_c$)	$0.5$
Surface double layer capacitance ($C_{\mathrm{DL}}$)	$0.2$ $\mathrm{F}·{\mathrm{m}}^{-2}$	[28]
Bruggeman exponent (p)	$1.0$	[34]
Current collector
Conductivity (${\sigma }_s$)	$3.7\times {10}^7$ $\mathrm{S}·{\mathrm{m}}^{-1}$	[34]
Carbon-binder domain
Conductivity (${\sigma }_s$)	375 $\mathrm{S}·{\mathrm{m}}^{-1}$	[34]
Electrolyte
Initial electrolyte concentration ($c_{e,0}$)	1000 $\mathrm{mol}·{\mathrm{m}}^{-3}$	[34]
Conductivity (${\kappa }_e$)	$f\left(c_e\right)$ $\mathrm{S}·{\mathrm{m}}^{-1}$	[37]
Diffusivity ($D_e$)	$f\left(c_e\right)$ ${\mathrm{m}}^2·{\mathrm{s}}^{-1}$	[37]
Activity ($\partial lnf/\partial ln{c}_e$)	$0.43$	[37]
Transference ($t_{+}^0$)	$0.37$	[37]

Table 4. Mechanical parameterization of the half-cell model.

Parameter	Value	Ref.
NMC622 particles
Partial molar volume ($\Omega$)	$1.8\times {10}^{-6}$ ${\mathrm{m}}^3·{\mathrm{mol}}^{-1}$	[34]
Young’s modulus (E)	140 GPa	[34]
Current collector
Young’s modulus (E)	70 GPa	[34]
Carbon-binder domain
Young’s modulus (E)	$0.3$ GPa	[34]
Electrolyte
Instantaneous shear modulus (G)	$0.3$ MPa	[38]
AM-CBD interface
Spring constant per unit area in pristine state ($K_{int,0}$)	$2\times {10}^{15}$ $\mathrm{N}·{\mathrm{m}}^{-3}$	[28]
Decreasing rate of spring constant per cycle ($k_d$)	$0.375\times {10}^{15}$ $\mathrm{N}·{\mathrm{m}}^{-3}$	[28]

Table 4. Mechanical parameterization of the half-cell model.

Parameter	Value	Ref.
NMC622 particles
Partial molar volume ($\Omega$)	$1.8\times {10}^{-6}$ ${\mathrm{m}}^3·{\mathrm{mol}}^{-1}$	[34]
Young’s modulus (E)	140 GPa	[34]
Current collector
Young’s modulus (E)	70 GPa	[34]
Carbon-binder domain
Young’s modulus (E)	$0.3$ GPa	[34]
Electrolyte
Instantaneous shear modulus (G)	$0.3$ MPa	[38]
AM-CBD interface
Spring constant per unit area in pristine state ($K_{int,0}$)	$2\times {10}^{15}$ $\mathrm{N}·{\mathrm{m}}^{-3}$	[28]
Decreasing rate of spring constant per cycle ($k_d$)	$0.375\times {10}^{15}$ $\mathrm{N}·{\mathrm{m}}^{-3}$	[28]

3. Results and Discussion

To visualize the expansion of NMC particles upon lithiation and the subsequent temporally and spatially resolved detachment along the contacts between the AM particles and the CBD, the half-cell’s first discharge is simulated in the pristine state. The half-cell is initially in the fully charged state and is then galvanostatically discharged at 1 C-rate to 3 V. The applied current at different C-rates can be calculated as expressed below.

$$I_{\mathrm{app}}=\frac{c_{s,max}·{\mathrm{V}}_{\mathrm{particles}}^{\mathrm{total}}·\mathrm{F}·{\mathrm{C}}_{\mathrm{rate}}}{3600\mathrm{s}}$$

where ${\mathrm{V}}_{\mathrm{particles}}^{\mathrm{total}}$ represents the total volume of the NMC particles in the reconstructed cathode electrode. While Figure 4a,b illustrate the state of lithiation (SoL) in the NMC particles across the cathode electrode at the beginning and end of the discharge process, their corresponding deformations are depicted in Figure 4c,d. At the beginning of discharge, SoL looks significantly uniform in the particles across the electrode, and as the discharge continues, non-uniformity becomes more noticeable. As can be seen in Figure 4b, some particles are still not highly lithiated at the end of the discharge process, indicating partially isolated particles, as also observed in previous experimental works on real electrodes [39,40,41]. Such heterogeneous electrochemical response in the electrode is due to the spatial arrangement and orientation of the different constituent domains in the electrode, resulting in different reaction kinetics.

Figure 4d shows the deformed electrode at the end of discharge compared to the almost undeformed electrode at the beginning of the lithiation process, Figure 4c. The applied free boundary condition in the model, as shown in Figure 3, results in the expansion and contraction of the electrode in the thickness direction. Such displacement can form tension stress perpendicular to the contact boundaries between AM particles and CBD, depending on the spatial positioning and arrangement of the constituent domains. Accordingly, Figure 5a exhibits the interfacial area, which undergoes detachment at SoL $=50\%$ upon the first discharge process. In agreement with what has been observed by means of the X-ray contrast tomography technique conducted by Xu et al. [28] and X-ray computed tomography (XCT) by Parks et al. [42], both the fractured AM particles and the detached contacts take place mostly close to the separator upon cycling. Moreover, the portion of contacts experiencing debonding directly relies on the intensity of the spatial mismatch between the mechanical response of the variously positioned and orientated AM particles and the CBD at each lithiated state of NMC particles.

To study the effect of the degrading contacts between AM particles and CBD on electrochemical performance during cycling, the fully charged half-cell is galvanostatically cycled at 1 C-rate for five sequential times in a defined cell voltage window of 3–4.3 V. Figure 5b exhibits the contact area experiencing decohesion at SoL = 50% during discharge at the fifth cycle. As noted, compared to the pristine state, a larger area undergoes detachment due to the degraded contact strength between AM particles and CBD during cycling. This outcome is in agreement with what was observed in SEM images conducted by Xu et al. [28], shown already in Figure 1, where the gap development between active and inactive materials occurs in the electrode during cycling. Moreover, Figure 6 illustrates the evolution of the cell characteristics during cycling. While Figure 6a exhibits the cell voltage and OCV, Figure 6b,c display the evolution of the averaged values for the state of lithiation and volumetric strain over the whole NMC particles during the cycling. In Figure 6d, the averaged interfacial gap associated with the detached contacts at the interfaces between AM particles and CBD is shown. As can be observed, the degrading contact during cycling, which is defined by Equation (18), results in a growing detachment as cycling proceeds. In other words, the diminishing spring constant through cycling leads to not only a larger detached area (Figure 5), but also to the formation of a larger decohesion along these contacts. Furthermore, the interfacial gap in Figure 6d shows an irreversible growth of detachment during cycling. This can be explained by the degradation of interfacial strength during cycling. The calculated interfacial gap maxima in Figure 6d varies from almost 10 nm in the first cycle to 20 nm in the fifth cycle. Liu et al. [32] performed a numerical study and attained a range of 50–65 nm for an NMC electrode cycled at 1 C-rate for five times. This discrepancy primarily originates from two different assumptions in the model built by Liu et al.: a fully clamped boundary condition was assigned to all RVE sides, and the domains of pore and CBD were simplified and represented by a single integrated composite domain with considerably high mechanical stiffness. While the former results in higher absolute values of interfacial gap during cycling, the latter leads to a narrower interfacial gap range compared to the computed values in this work.

The effect of interfacial detachment on the polarization of the cell, and hence on the capacity retention, is better illustrated in Figure 7, which shows the delivered discharge capacities, normalized by the theoretical capacity, during the cycling. The gradual decrease in discharge capacity during cycling is due to increasing impedance caused by increasing detachment. However, the more significant decrease in capacity from the first to the second discharges is mainly due to the higher initial cell voltage and more complete lithium insertion in the first discharge.

Electrochemical Impedance Spectroscopy (EIS)

As noted earlier, the development of such spatially resolved gaps along the contacts hinders the electron transfer paths and impedes the reaction kinetics. The resulting elevated overpotential caused by the interfacial debonding after only five cycles and its subsequent capacity retention, depicted in Figure 7, indicates that the delivered discharge capacity drops smoothly over the successive cycles. On the other hand, since electrochemical interactions are biased by various degradation mechanisms, electrochemical impedance spectroscopy is often used to characterize the deterioration of the battery’s performance. Therefore, to quantify the effect of interfacial decohesion of AM particles and CBD on the performance of the half-cell, the impedance spectra are also numerically computed in different cycles. The EIS simulation is conducted around SoL = 50% as the equilibrium state. Then, a harmonically oscillating voltage signal with varying frequency and small magnitude is introduced to perturb the half-cell. While the input signal magnitude is set to 10 mV, its frequency is varied from 1 mHz to 100 kHz.

Figure 8 illustrates the computed impedance spectra for different studied cases. The calculated impedance spectrum is composed of two semicircles in the high- and mid-frequency ranges, followed by a sloping tail in the low frequencies. To quantify the effects of interfacial detachment on the dynamic response, the EIS simulation was initially carried out under the two conditions “with interfacial gap” and “without interfacial gap” (Figure 8a). As can be observed, when interfacial detachment is accounted for, the spectrum differs from the case with well-bonded contact. This variation primarily includes a shift of the spectrum to the right and also a larger dimension of the high-frequency semicircle, both due to the existing debonded area and hindered electronic transport. Moreover, the mid-frequency semicircle enlarges slightly, which is explained by the impeded reaction kinetics caused by the introduced local drop in the solid-phase potential at the particle surfaces detaching from the CBD. On the other hand, Figure 8b exhibits a similar but significantly more gradual change in the impedance spectra during cycling, which is due to the mechanical deterioration of the interfacial cohesion. Similar to the observed slight capacity decrease in Figure 7, the impedance spectrum slightly alters while cycling. It is certain that the strength of the interface would continue to decrease by further cycling after the fifth cycle, resulting in greater detachment at the interface and thus higher impedance.

4. Conclusions

Mechanical degradation is an obstacle on the way to high-energy-density electrodes. Different types of mechanisms lead to mechanical degradation of the LIB cell by contributing to the capacity fade and impedance increase. The numerical methods provide a low-cost and fast alternative to investigate mechanical degradation as compared to experimental tomography techniques. Such methods play an even more important role in the case of cyclic mechanical fatigue, which requires more resources for experimental testing.

In this context, a 3D half-cell’s representative volume element (RVE) was built using the virtual morphology of the NMC622 cathode electrode to allow for exploration of the degrading contact of the AM particles and CBD in the electrodes. The half-cell was cycled and its characteristics during cycling were visualized. The degrading interfacial strength over five cycles was quantified in terms of the evolution of the interfacial gap in the time domain and the increased impedance in the frequency domain by means of numerical electrochemical impedance spectroscopy. In agreement with previous experimental tomography images, it was observed that most of the detachments occur in the regions close to the separator. Moreover, it was found that by cycling and degradation of contact between the AM particles and CBD, the impedance spectrum shifts to the right and the high-frequency semicircle enlarges. Furthermore, the mid-frequency semicircle’s dimension expands as well. The former changes in the spectrum are due to the more retarded electronic transport paths, while the latter reflects the introduced drop in the solid-phase potential, hindering the reaction kinetics.

The framework presented in this work can be further employed to investigate different interfacial characteristics in LIBs. For instance, it enables determining the minimum required CBD cohesion strength, the optimized CBD volumetric content, and the optimized CBD morphology. Additionally, the framework can be used to study the critical role of delamination of the solid electrolyte (SE) in solid-state batteries. Such investigations provide a conceptual design in the early stages of a research project that allows optimization concerning desired functionality. In other words, virtual material testing and experimenting with various properties of the electrode constituents is possible in short time and with little cost.