Electrode Mesoscale as a Collection of Particles: Coupled Electrochemical and Mechanical Analysis of NMC Cathodes

NMC active material particles for a 92 wt% NMC333, 4 wt% PVDF, and 4 wt% carbon black cathode calendered under a 300 bar pressure are derived from available XCT image data.²⁷ Ebner et al.²⁷ provide a set of greyscale images where each NMC particle is separately labeled by color. The Avizo 9 software (FEI; Hillsboro, OR) was used to first expand each particle by approximately one voxel to account for watershed errors and to match binary segmented volume fractions. Then, each particle was separately segmented in Avizo using the greyscale labels. A surface mesh was created on each particle, smoothed to remove unphysical voxellation stair-steps, and exported in a standard tessellation library (STL) format.

A simulation domain of 50 × 50 × 50 μm from within the set of imaged particles was selected and meshed with tetrahedral elements using Cubit.⁸¹ Additional mesh domains representing the current collector and separator were also added, as shown in Fig. 1c. This mesh of the domain is called the background mesh. A mesh size of 0.5 μm is used to be consistent with previous solution verification work.²⁰ Roberts et al.²⁰ do recommend a larger domain than (50 μm)³, however we plan to study the variability of domain size and selection in our future work.

The background mesh is decomposed into particle and non-particle phases using the Conformal Decomposition Finite Element Method (CDFEM).^20,82,83 Each particle STL file is used to initialize a level-set field on the background mesh, upon which the CDFEM algorithm decomposes elements containing level-set zero crossings into phase-conformal child elements. This approach has been used heavily in previous publications.^50,53,54

Subsequently, the CBD phase is added to the domain using the "binder bridge" approach proposed by Trembacki et al.^19,21,22 This algorithm places CBD into porous, consolidated domains that are concentrated near particle-particle contacts, attempting to replicate the morphologies observed by Jaiser et al.⁸⁴ The total volume of composite binder was set to match the reported weight fractions from the image data,²⁷ additionally accounting for 47% nanoporosity within the CBD phase.^52,85

The final computational domain is composed of five material phases (illustrated in Fig. 2): a current collector, active material particles (NMC333), composite binder domain (CBD), electrolyte, and a separator. Properties for each constitutive material are displayed in Table I. The separator phase is assumed to have the same intrinsic transport properties as the electrolyte, as they are typically very porous to facilitate ion transport.

Zoom In Zoom Out Reset image size

Table I.  Material properties used in the mathematical model.

Property	Symbol	Value	Units	Source
Li Diffusivity	DLi	2.51 × 10−10	cm2 s−1	86
Li+ Diffusivity	${D}_{{\mathrm{Li}}^{+}}$	2.6 × 10−6	cm2 s−1	87, 88
${\mathrm{PF}}_{6}^{-}$ Diffusivity	${D}_{{\mathrm{PF}}_{6}^{-}}$	4.88 × 10−6	cm2 s−1	87
Charge Transfer Coefficient	αi	0.5	—
Reaction Rate Constant	k	4.38 × 10−11	m2.5 s−1-mol−0.5	89
Maximum Li Concentration	Cmax,Li	3.61 × 104	mol m−3	89
NMC Young's Modulus	ENMC	138	GPa	90
NMC Poisson's Ratio	νNMC	0.3	—	91
Volumetric Expansion	β	4.97 × 104	Pa mol−1-m−3	92
CBD Young's Modulus	ECBD	0.30	GPa	91
CBD Poisson's Ratio	νCBD	0.4	—	91
NMC Conductivity	κNMC	f(CLi) Eq. 2	S cm−1	86
CBD Electrical Conductivity	κCBD	f() Eq. 3	S cm−1	19
CBD Porosity	φ	0.47	—	52
CBD Tortuosity	τ	9.4	—	52, 93
Equilibrium Potential	ϕeq	f(CLi) Eq. 6	V	94
Faraday Constant	F	96485	s-A mol−1
Gas Constant	R	8.314 5	J mol−1-K−1
Temperature	T	298.0	K

Electron transport in the particles and CBD is governed by Ohm's Law for current conservation,

$$\begin{eqnarray}&&{\boldsymbol{\nabla }}\cdot ({\kappa }_{i}{\boldsymbol{\nabla }}{\phi }_{s})=0,\end{eqnarray} \tag{ 1 }$$

where κ_i is the electrical conductivity of material i and the solution variable ϕ_s is the solid phase voltage. Amin and Chiang⁸⁶ report that MC conductivity has a heavy dependence on lithiation, ranging over multiple orders of magnitude through the course of discharge. Thus, we use a fit for conductivity (κ_NMC):

$$\begin{eqnarray}&&\begin{array}{l}\mathrm{log}[{\kappa }_{\mathrm{NMC}}/({\rm{S}}\,{\mathrm{cm}}^{-1})]=-354.49{y}^{6}+1061.8{y}^{5}-1261.8{y}^{4}\\ \quad \ +\ 757.19{y}^{3}-242.81{y}^{2}+43.668y-7.3248,\end{array}\end{eqnarray} \tag{ 2 }$$

where y = 1 − C_Li/C_Li,max. Additionally, binder conductivity is dependent on strain (), and we use the "fresh binder" model with a 10× multiplier published by Trembacki et al.¹⁹

$$\begin{eqnarray}&&{\kappa }_{b}=0.1593-173.97\epsilon \end{eqnarray} \tag{ 3 }$$

Lithium transport within the particles is governed by,

$$\begin{eqnarray}&&\displaystyle \frac{\partial {C}_{\mathrm{Li}}}{\partial t}+{\boldsymbol{\nabla }}\cdot {{\bf{J}}}_{\mathrm{Li}}=0,\end{eqnarray} \tag{ 4 }$$

where C_Li is the lithium concentration and ${{\bf{J}}}_{\mathrm{Li}}$ is the lithium flux. Lithium transport within NMC follows non-ideal solution theory due to preference of lithium insertion at specific intercalation sites^53,54,95:

$$\begin{eqnarray}&&{{\bf{J}}}_{\mathrm{Li}}={{MFC}}_{\mathrm{Li}}{\boldsymbol{\nabla }}{\phi }_{\mathrm{eq}}.\end{eqnarray} \tag{ 5 }$$

Here, M = D_i/RT is the mobility of lithium, given by the Nernst-Einstein relation, ϕ is the equilibrium potential, F is Faraday's constant, D_Li is the lithium diffusivity, R is the gas constant, and T is the temperature. The equilibrium potential is from data published by Wu et al.⁹⁴ for a NMC half-cell discharged at a C/25 rate, to which we apply the following polynomial fit with an added exponential tail to represent the portion of the curve near full lithiation

$$\begin{eqnarray}\begin{array}{rcl}{\phi }_{\mathrm{eq}}(x) & = & -2960.98{x}^{7}+14272.3{x}^{6}-29127.0{x}^{5}\\ & & +\,32600.0{x}^{4}-21599.4{x}^{3}\\ & & +\,8471.45{x}^{2}\,-\,1823.91x\,+\,170.967-\exp \left[250.0(x\,-\,1)\right],\end{array}\end{eqnarray} \tag{ 6 }$$

where x is the local lithiation fraction C_Li/C_Li,max.

Current and ionic species conservation are conserved in the porous CBD and electrolyte phases through Eqs. 7 and 8, respectively,

$$\begin{eqnarray}&&{\boldsymbol{\nabla }}\cdot ({{\bf{J}}}_{{\mathrm{Li}}^{+}}-{{\bf{J}}}_{{\mathrm{PF}}_{6}^{-}})=0,\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray}&&\varphi \displaystyle \frac{\partial {C}_{{\mathrm{Li}}^{+}}}{\partial t}+{\boldsymbol{\nabla }}\cdot {{\bf{J}}}_{{\mathrm{Li}}^{+}}=0.\end{eqnarray} \tag{ 8 }$$

Here, ${C}_{{\mathrm{Li}}^{+}}$ is the lithium-ion concentration and φ is the porosity of the medium (1 for electrolyte, 0.47 for CBD). The current conservation, Eq. 7, contains contributions from both Li⁺  and ${\mathrm{PF}}_{6}^{-}$ whose fluxes are governed by the Nernst–Planck model,

$$\begin{eqnarray}&&{{\bf{J}}}_{i}=-\displaystyle \frac{\varphi {D}_{i}}{\tau }\left[{z}_{i}{C}_{i}\displaystyle \frac{F}{{RT}}{\rm{\nabla }}{\phi }_{l}+{\rm{\nabla }}{C}_{i}\right].\end{eqnarray} \tag{ 9 }$$

Here, τ is the tortuosity, D_i is the ionic diffusivity, and z_i is the species valence. Finally, we assume electroneutrality for the ${\mathrm{PF}}_{6}^{-}$ concentration, ${C}_{{\mathrm{Li}}^{+}}={C}_{{\mathrm{PF}}_{6}^{-}}$.

Reaction kinetics are modeled by applying a Butler-Volmer current density (j) as a normal flux boundary condition for potential on both the solid and electrolyte sides of the interfaces between particles and CBD/electrolyte:

$$\begin{eqnarray}&&j={i}_{0}F\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{ 10 }$$

Here, the exchange current density is given by ${i}_{0}=\varphi k{\left({C}_{{\mathrm{Li}}^{+}}\right)}^{{\alpha }_{a}}{\left({C}_{\max ,\mathrm{Li}}-{C}_{s}\right)}^{{\alpha }_{c}}{C}_{s}^{{\alpha }_{c}}$, η = ϕ_s − ϕ_l − ϕ_eq is the overpotential, α_i are the charge transfer coefficients, and k is the rate constant. Concentration flux (r) is tied to the current density through reaction stoichiometry, meaning that r = j/F is applied as a normal flux condition for concentration on the same surfaces. The system is fully defined by a Dirichlet boundary condition on the liquid phase voltage at the free end of the separator (ϕ_l = 0) and a current flux boundary condition at the free end of the current collector to match the desired C-rate. A bulk boundary condition at the separator interface (${C}_{{\mathrm{Li}}^{+}}={C}_{{\mathrm{Li}}^{+},0}$) effectively provides an infinite source of Li⁺ . While we recognize that there could be concentration gradients in Li⁺  within the separator that impact this bulk boundary condition,⁹⁶ because of the low concentration gradients we see within the electrode itself, we choose this simpler approach. The cathode is considered charged in its initial state, set to a spatially uniform C_{Li, 0} = 0.42C_Li,max. Initial electrolyte concentration is defined to be 1.2M LiPF₆.⁹⁷

Mechanical deformation of the cathode particles is governed by the equation

$$\begin{eqnarray}&&{\boldsymbol{\nabla }}\cdot {\boldsymbol{\sigma }}=0\end{eqnarray} \tag{ 11 }$$

where ${\boldsymbol{\sigma }}$ is the Cauchy stress tensor. The stress includes contributions from linear elastic mechanical deformation as well as charge-induced swelling,

$$\begin{eqnarray}&&{\boldsymbol{\sigma }}=\lambda \mathrm{tr}({\boldsymbol{\epsilon }}){\bf{I}}+2\mu {\boldsymbol{\epsilon }}-\beta {\rm{\Delta }}{C}_{\mathrm{Li}}{\bf{I}}.\end{eqnarray} \tag{ 12 }$$

Here, λ and μ are Lamé parameters, ${\boldsymbol{\epsilon }}=[{\boldsymbol{\nabla }}{\bf{d}}+{\left({\boldsymbol{\nabla }}{\bf{d}}\right)}^{{\rm{T}}}]/2$ is the small strain deformation tensor, and ${\bf{d}}$ is the mesh displacement vector.

The term βΔC_Li in Eq. 12 represents lithiation-induced swelling, following Vegard's law,⁹⁸ where β assumes isotropic swelling of the NMC secondary particle structures. This parameter is obtained through a process outlined by Malavé et al.⁹⁹ which utilizes data of unit cell volume as a function of lithiation, $\beta =\tfrac{\chi K}{{C}_{\max ,\mathrm{Li}}}$. Here, K is the bulk modulus, C_max,Li is the maximum concentration in the particle phase, and χ is the slope of the strain vs lithium fraction curve. A linear fit for the unit cell volume was taken between lithium fraction x = 0.253 and x = 1.0, with x = 0.253 taken to be the fully charged and thus zero stress state.

Identical mechanical deformation equations were applied to the CBD phase, but with β = 0. Mechanical behavior in this phase is especially important due to the strain-dependence of electrical conductivity.

In all cases, no-penetration boundary conditions, n · d = 0 are applied to the edges of the computational domain. This fixes domain boundaries to allow no macroscopic swelling of the cathode, ensuring that all swelling of NMC particles results in compression of the CBD/electrolyte phases. While this restriction can help isolate the influence of mesostructure on local stresses, we acknowledge that real cathodes do exhibit macroscopic swelling over the course of discharge, an effect we investigate later.

The mathematical model is implemented and solved using Sandia's Sierra/Aria Galerkin Finite Element Model code.¹⁰⁰ The model is formulated into two separate equation systems; one containing the electrochemical model and the other the mechanics model. Each equation system is solved monolithically using a full-Newton method to compile the linear system and iteratively solve the nonlinear system of equations. The linear system of equations is solved iteratively using GMRES with a DD-ILU(T) preconditioner. The transient problem is incremented using a second-order backward differentiation formula time scheme (BDF2) with an adaptive time stepping algorithm to control the error in the predictor portion of the time integration algorithm.

The simulation domain contains 2.4 million tetrahedral elements and 0.4 million nodes, resulting in 1.5–1.6 million degrees of freedom in each linear system. Each simulation was run on 30 processors and the simulation was discretized into 500 time steps.

While visually electrodes at the mesoscale provide significant insight, it is often the macroscale behavior of the electrode that is the most relevant to practitioners. Therefore, we first predict some macroscale quantities, shown in Fig. 3.

Zoom In Zoom Out Reset image size

The most obvious quantity to study is the discharge curve; voltage as a function of state of charge compared to the open circuit voltage, shown in Fig. 3a. Two features are immediately noticeable. First, the cell shows polarization losses on the order of 100 mV–200 mV throughout the entire discharge, primarily due to the large current density. Based on the various phenomena shown in Fig. 5, these losses are due to a combination of charge transfer kinetics, electrical transport, and ionic transport. Electrical transport losses become more important later in the discharge as fewer particles are reacting. Second is capacity loss; the cell reaches a cutoff voltage of 3.4V at approximately 90% of full discharge, representing a noticeable (although reversible) loss in capacity.

Lithiation-induced swelling generates local stresses throughout the cathode domain, which are trasmitted through the mesostructure to impose a macroscale pressure on the current collector and separator (Fig. 3b). As the stress generated from lithiation is linear with state of charge, so is the resulting macroscale pressure, reaching a maximum of nearly 80 MPa by the end of discharge. This mean pressure is significantly lower than some of the localized stresses that will be discussed later. We note that all external boundaries of the simulation domain are fixed, preventing the stress from relieving itself by expanding into the separator or macroscopically deforming the cell, an effect that will be addressed later.

Mesoscale modeling is commonly used to predict the effective electrical conductivity of the electrode for macroscale modeling,⁶³ treating the effective conductivity as a static parameter. Recent works^19,21,22 have recognized that the CBD conductivity changes as a function of lithiation due to the mechanical forces generated from lithiation. Absent a full electrochemical model, however, these works must assume a constant state of lithiation in the particles. In Fig. 3c, we instead present the instantaneous effective electrical conductivity as a function of state of charge, fully considering the inhomogeneous state of charge of each particle. The effective electrical conductivity is calculated by separately applying a voltage drop across the thickness of the electrode, calculating the current flux into the collector, and then using ${\kappa }^{\mathrm{eff}}=-{\bar{J}}_{\mathrm{elec}}L/{\rm{\Delta }}V$, where ${\bar{J}}_{\mathrm{elec}}$ is the mean current density, L is the thickness of the electrode, and ΔV is the applied voltage difference.

Throughout discharge, the cathode's effective conductivity demonstrates a nearly linear increase with time. This is largely due to the linear increase in conductivity of the composite binder, coupled with the linear stress increase shown in Fig. 3b, as it is compressed by the surrounding swelling NMC particles. Deviation from a perfectly linear relationship occurs because our expression for strain-dependent binder conductivity is capped at 5 S cm⁻¹ to match the range of values reported by Trembacki et al.,¹⁹ thus causing the rate of increase to slow as more local hotspots reach this maximum value.

Impact of mechanical boundary conditions

Results in Fig. 3 can be affected by the application of mechanically rigid boundary conditions on all boundaries, most importantly for the separator and current collector boundaries, effectively limiting swelling in the z-direction. Allowing compression of the separator can introduce a flexibility that mimics actual cell stiffness. Figure 4 shows the effect of introducing a stiffness parameter k = EA/L to the separator phase, where E is the modulus, A is the cross-sectional area, and L is the thickness of the phase. The vertical dotted line represents the equivalent stiffness of a swollen polymeric separator, where k = 2.04 × 10⁴ N m⁻¹ is calculated from mechanical properties reported by Gor et al.¹⁰¹ The vertical, dashed line is the effective stiffness of the full cathode (k = 9.66 × 10⁶ N m⁻¹).

Zoom In Zoom Out Reset image size

Figure 4a displays the volume-average von Mises stress calculated over all particles in the domain. As this value reaches a plateau, the system begins to resemble the rigid boundary conditions used previously. The mean von Mises stress at high stiffness values is over 110 MPa. However, at low stiffness values, the mean stress decreases to around 70 MPa, a 35% decrease. At low k, the separator stiffness does not contribute, and the average stress is primarily due to particle contact within the mesostructure and the impact of the in-plane no-displacement constraints. The transition between these two plateaus occurs near where the imposed stiffness is between the stiffness of the separator (the softest material in the system) and the stiffness of the cathode itself. These pressures and stresses are transmitted through the separator to the anode, where they could have important implications on phenomena such as dendrite growth and lithium plating.^102,103

Figure 4b provides a detailed breakdown of where the volume changes primarily occur. For the rigid boundary case (high stiffness), the particle phase increases in volume around 0.35% halfway through discharge. This projects to be significantly less than the 3.2% volume change at full discharge shown by Kam et al.,⁹² which was used to calibrate the swelling of a single particle in our simulations. However, results are similar to those shown by Mendoza et al.⁵³ This is due to the combination of mesoscale and macroscale confinement of these elastic particles (illustrated by the high stresses in Fig. 5) and highlights the importance of mesoscale modeling over single-particle modeling. In this high stiffness limit, all of the volume increase of particles goes into decreasing the volume of the non-particle phases. While the CBD phase does exhibit a small amount of compression (up to 0.1%), the liquid electrolyte phase contributes the majority of this volume change.

Zoom In Zoom Out Reset image size

As the boundary stiffness decreases, a number of changes are observed. First, the overall volume change of the particles increases slightly, to 0.37%. While the z direction is not constrained, there is still a significant mesoscale confinement effect in the particle network, manifesting itself by transmitting stress to the x and y boundaries. Also, as stiffness decreases, the volume change of the non-particle phases (CBD and electrolyte) actually decrease, even though the overall particle expansion is increasing. This is because the cathode has room to grow in the z direction, evidenced by a 0.28% increase in cathode volume.

Mesoscale simulation allows us to probe a number of different quantities relevant to cell performance. While we will quantify many of these quantities later, it is instructive to first visualize the relevant fields, as is illustrated in Fig. 5. These plots represent planar cuts through the 3D mesostructure taken at different times (25%, 50%, 75%, and 92% of theoretical total discharge) over the course of a single discharge simulation.

Figure 5a shows the variation of lithium concentration within the particles of the mesostructure. To highlight concentration differences, we show this as a scaled concentration difference from the mean value, ${\rm{\Delta }}{C}_{\mathrm{Li}}=({C}_{\mathrm{Li}}-{\bar{C}}_{\mathrm{Li}})/{C}_{\mathrm{Li},\max }$, where ${\bar{C}}_{\mathrm{Li}}$ is the average lithium concentration at each given time. Smaller particles tend to discharge faster, meaning they display a concentration higher than the bulk average, while larger particles tend to respond slower to the applied current. This phenomenon is due to a larger surface area to volume ratio in the small particles, allowing the reaction to add more lithium per particle volume. This suggests that the reactions may be kinetically limited, which will be discussed in more detail later.

Additionally, we can see much more variance in particle concentration at 50% and 75% SOC than in the final snapshot. Small particles continue to more fully lithiate throughout the middle portion of the discharge, exacerbating the difference at those snapshots. However, toward the end of discharge, these particle become fully lithiated and no longer participate in the reaction, allowing the larger particles to "catch up" to their state of lithiation.

Finally, the largest particles show a measurable internal concentration gradient, particularly in moderate to late times. Lithium diffusion is driven by gradients in the electrochemical potential Eq. 5. In the middle portion, the equilibrium voltage curve becomes relatively flat (i.e. doesn't change much with concentration), making it harder to diffuse lithium and allowing gradients to build up.

Figure 5b displays the concentration of lithium ions in the electrolyte and porous CBD phases. The CBD/electrolyte boundary is shown as a black outline. This illustration shows the expected trend of Li⁺ diffusion, as there is a general decrease in concentration further away from the separator, exhibiting a 25% decrease across the thickness of the electrode. Furthermore, it is evident that areas with larger masses of binder display slightly lower ionic concentration, signifying more difficult transport through this porous phase (as it has only 5% of the diffusivity of pure electrolyte). However, subsequent sections will show that ionic gradients both across the cathode and within the porous CBD do not pose a limitation to the electrochemical reaction and thus do not measurably affect the macroscopic discharge performance. Note that literature in this area has shown that electrode thickness can have a significant effect on cathode utilization, particularly for thicknesses greater than 60 μm.^14,78 While ionic gradients are not impactful in the domains studied here, mesoscale transport effects in thick electrodes warrant further study.

Over the course of discharge, this profile remains fairly consistent with no noticeable differences from 25%–75% SOC. However, near the end of the simulation, the gradient becomes much more pronounced near the largest particle in the domain. This is due to the variation in the state of lithiation for differently sized particles, as previously discussed. Near the end of the discharge, the smaller particles have already been fully lithiated. Thus, the larger particles assume all of the remaining reactive flux from the applied current. This forces a larger reaction rate at the surface of the few remaining un-lithiated particles, consuming Li⁺ faster than it can diffuse. In the 92% discharged snapshot, it is evident that the Li⁺ concentration is smaller near the largest particle due to this effect. This will, in effect, increase the local overpotential of the few particles continuing to contribute to the reaction and lower the macroscopically observed voltage (as will be shown later). Ultimately, this is evidence of local mesostructure influence on cell performance and highlights the need for simulation at this scale.

Mesoscale resolution also allows us to see the electrical current pathways through the system over the course of discharge, as illustrated by the magnitude of the current density in Fig. 5c. As the composite binder is much more conductive than NMC, there is a much stronger presence of current through the binder than the active material. Thus, the most efficient pathway to carry current from the collector to a reactive surface will be one that includes the most binder. However, there is still a large driving force for the electrochemical reaction at all NMC surfaces, whether coated with CBD or simply in contact with electrolyte. For electrons to reach the particle/electrolyte surfaces, they must pass through the low conductivity NMC, which can be seen in each snapshot as diffusive current regions near electrolyte interfaces in the larger particles.

Furthermore, this affects the voltage drop seen in Fig. 5d, which is displayed as the local solid or liquid phase voltage compared to the voltage at the relevant boundary (ϕ_s at the current collector, ϕ_l at the separator). As the current is directed through the active material, the increased resistance creates a significant voltage drop (as high as 100 mV) that can be seen by much higher solid phase potentials on particle/electrolyte surfaces. At late times, the largest particle has substantial current throughout, as it does not have enough CBD pathways to support the large reaction rate it sustains late in the discharge, also contributing to a drop in the observed macroscale voltage.

CBD has a strain-dependent electrical conductivity, which also plays a significant role in electrical current transport, as shown in Fig. 5e. As the NMC particles lithiate, they will swell, imposing a mechanical stress on the composite binder and deforming it. When this binder phase compresses, the carbon black becomes more densely packed, thus increasing conductivity. This phase is already much more conductive than the active NMC particles, but this difference is exaggerated as the cell is discharged. NMC lithiation decreases its conductivity, and its lithiation-induced swelling increases conductivity in the binder phase. By the end of the simulation, electron transport is extremely favored in the binder phase. It is also worth noting that the snapshots in line four are also a direct representation of binder strain, as it shares a linear relationship with conductivity.

This binder strain is a direct manifestation of the mechanical stresses generated by lithiation, shown using the Von Mises stress invariant in Fig. 5f. As was demonstrated by Mendoza et al.⁵³ and Roberts et al.,⁵⁴ even, though the articles are nearly spherical, their internal stresses are not uniform. Instead, is is the arrangement of particles at the mesoscale (combined with the mechanical boundary conditions applied) that constrains particle motions and leads to strong stress-supporting chains throughout the percolated particle network, highlighted by the stress concentrations near particle-particle contacts. These mesoscale stress effects strongly outweigh any stresses that would be observed for isolated particles. It is also the inter-particle effects that lead to the binder strain and increased binder conductivity, highlighting the importance of image-based mesoscale simulations for truly understanding performance of these battery materials.

Perhaps the most significant advantage of the present mesoscale approach enabled by our meshing technique is the ability to extract data on a per-particle basis. This allows exploration of how particle attributes such as size or location can impact the cell characteristics during discharge.

In Fig. 6a, individual mean particle concentration is plotted vs particle size (volume to surface area ratio) at a single timestep. Here, it is clear that there is a general trend suggesting that smaller particles are at a higher concentration, which is consistent with results illustrated in Fig. 5. This apparent exponential decay would be expected for kinetically limited problems, where there are not strong concentration gradients within the particles and the concentration is limited by reaction rate at the surface.

Zoom In Zoom Out Reset image size

However, there is significant scatter apparent in Fig. 6a with a large number of small particles displaying a lower concentration than the trend would suggest. We hypothesize two possible explanations or this scatter. First, if the reactions are primarily kinetically limited, then the surface area may not be the most accurate measure of reaction area, as some parts of the surface area are occluded by contact with other particles. Therefore, Fig. 6b has an x-axis where the volume is scaled by the particle active area (particle surface touching either electrolyte or porous CBD phases). While this somewhat reduces the scatter and makes the trend more pronounced, the scatter is still significant.

The second hypothesis is that particles that are touching allow lithium diffusion between them, making them act as a single apparent particle in terms of lithium volume for diffusion, while having a drastically different apparent area. Our mathematical model does not have any resistance to lithium diffusion between particles, making this plausible. The impact of a transport resistance between particles is an area of future study. To explore this, Fig. 6b shows particles in two colors. The blue markers are single particles; those not directly touching any other particle. The red markers are particles that are touching at least one other particle. Clearly the large majority of particles that do not fall on the expected volume to area trend line are contacting particles.

To further investigate this, we group particles together if they share any contact. If there was an unbroken diffusive pathway from one particle to another, they were considered part of the same "group" and displayed as such on Fig. 6c. With all grouped particles collapsed into a single marker, the trend is very clear. This suggests that the most important factor in a particle's discharge is not a transport property such as diffusivity or conductivity, but in fact is simply due to its surface area to volume ratio. Ultimately, it is the surface reaction rate which governs how quickly a particle will fill with lithium, and its volume determines the capacity to which it is filled.

Figures 6d–6f show that this trend remains throughout the course of the simulation, but undergoes an inflection in its concavity as time progresses. This phenomena was similar to the previous discussion about concentration in Fig. 5 and is due to smaller particles reaching capacity while larger particles still are significantly de-lithiated.

The process by which the per-particle concentration curves smooth out is further analyzed in Fig. 7. Subfigures b and c show per-particle and per-group concentrations, respectively. Each line in these subfigures is on a color scale determined by the particle (or group) volume. It is evident that smaller, disconnected particles discharge more quickly than larger particles. The value of ∂C/∂t in a given particle is also linked to the derivative of the OCV curve, shown in Fig. 7a. Ultimately, this derivative is defined by ∂C/∂t = r · (A/V), where r is the reaction rate, A is the particle surface area, and V is the particle volume. Here, r is determined by the Butler–Volmer boundary condition Eq. 10, which is exponentially dependent on the local overpotential. As we will show later, the voltage drop on binder/particle surfaces is fairly uniform throughout the electrode, so ϕ_s is fairly consistent for each particle. However, there is significant variation in ϕ_eq and the difference between ϕ_s and ϕ_eq is what ultimately drives the electrochemical reaction. In the flat regions of the OCV curve, the overpotential will increase in magnitude because ϕ_s drops while ϕ_eq stays close to the same value. This pushes the reaction forward and allows the particle to discharge faster. This change in slope corresponds between Figs. 7a and 7b. Note that the other term contributing to ∂C/∂t is the surface area to volume ratio, which scales as 1/R, providing further explanation why smaller particles discharge faster. As illustrated in Fig. 7c, grouping contacting particles together results in a more clear trend highlighting the importance of particle capacity.

Zoom In Zoom Out Reset image size

Use of this mesoscale approach not only provides per-particle data, but it can isolate per-surface data as well. The CBD placement algorithm selectively places the CBD phase in areas near particle-particle contact. This means that there will be some particle surfaces covered with the polymeric binder and some surfaces in direct contact with electrolyte. Through our post-processing, we are able to show differences in transport on these surfaces and suggest how they affect global discharge.

Figure 8 illustrates the difference in reaction rate at each particle surface at two separate states of discharge. For a large portion of the discharge time, the reactive flux on the electrolyte surfaces is higher than that of binder surfaces (Fig. 8a). However, as the simulation nears completion, this is no longer the case as binder surfaces begin to react as fast or faster than the electrolyte area (Fig. 8b). This is due to the behavior we discussed for Fig. 5. Later in the simulation when NMC particles become lithiated, their conductivity decreases. This results in a polarization loss that increases ϕ_s on electrolyte surfaces. As these interfaces depend on electron transport through the active material for the reaction, they will ultimately see a smaller magnitude overpotential as the solid voltage gets closer to the OCV. This slows down lithium insertion from these surfaces, suggesting that particles not covered in enough binder could contribute to capacity loss at high C-rate.

Zoom In Zoom Out Reset image size

Augmenting this phenomenon is the increased conductivity of the CBD phase late in discharge due to the compressive forces it is subjected to. While it becomes harder to get electrons to the electrolyte-exposed surfaces to react, it simultaneously becomes easier to get electrons to the binder-covered surfaces. Even though there is a 95% diffusion penalty for lithium ions through CBD, Fig. 5 showed very little effect of this penalty, as the lithium ion diffusivity is sufficiently high.

Binder surfaces only react at locations in contact with pore space (half of the total surface). When correcting the active area to provide a mol s⁻¹ rate, as shown in Fig. 8c, it is apparent that the binder surfaces are comparable to the electrolyte surfaces throughout discharge and exceed their performance late in discharge. This means that there is not a transport limitation of Li⁺ to the particle surface, and the reaction is instead governed by the overpotential. This information can be useful when designing a cathode, as you simply need to consider the amount of available area for reaction when adjusting the weight fraction of porous binder during electrode manufacture.

Through our per-surface analysis, we can also look at the distribution of surface overpotentials and loss mechanisms which influence discharge characteristics. In Fig. 9a, total surface overpotential is plotted against particle size. As previously mentioned, this overpotential is fairly comparable for binder and electrolyte surfaces early on, but begins to deviate toward the end of the simulation with binder surfaces exhibiting higher overpotential and faster reaction rates. At 90% discharge (Fig. 9d), the discrepancy between surfaces is very apparent, particularly for the larger particles.

Zoom In Zoom Out Reset image size

Figures 9b, 9e show the electrical transport loss to each surface. Early in the discharge, when NMC conductivity is relatively high, the losses to both exposed and binder-covered surfaces is negligible (Fig. 9b). However, later in the discharge (Fig. 9e), NMC conductivity is sufficiently low as to lead to significant losses transporting electrons from the CBD to electrolyte surfaces to react. On the largest particles, the electrolyte exposed surfaces can experience as much as a 100 mV penalty. This is the primary contributor to the lower reaction rates for electrolyte surfaces shown in Fig. 9d.

Figure 9c shows polarization losses through ionic transport (blue squares) and composite binder conduction (red circles). These voltages are very small, and thus do not significantly contribute to polarization throughout the cell. However, we do see the expected trends as ϕ_liq increases toward the current collector and ϕ_s increases toward the separator. Here we should note that our analysis did not show any significant dependence for particle discharge on distance from the current collector. This is likely due to the very small values we observe for the voltage drop in the liquid phase, thus suggesting very good ionic transport.

The results in this section have thus far suggested that neither species nor electrical transport place a significant transport limitation on cathode performance over the length scales studied here. The most successful predictor of a single particle's discharge state is its size, and thus, capacity. If this is the case and no transport limitations are strongly present, the volumetric reaction rate for each particle should be constant. Figure 10 illustrates this volumetric rate while utilizing the same particle grouping method from Fig. 6. Figure 10a shows the expected trend for volumetric reaction rate, as it is nearly constant for all particle groups. Only the largest diffusive network exhibits a slight decrease from this value, suggesting that extremely large particles can experience some transport limitation (which is illustrated in Fig. 5a). This Behavior begins to change at 60% discharge, as the smaller particles have become so lithiated that their overpotential drops and rate of reaction begins to slow, causing the left portion of the curve to significantly decrease. Near the end of discharge, Fig. 10c, only the largest particles have the capacity to add more lithium, and the curve shifts to suggest the majority of reaction is occurring on these large groups. While this demonstrates that transport limitations can play a role late in discharge, the most significant contributor to particle discharge performance is simply the particle's size.

Zoom In Zoom Out Reset image size