Rapid Simulation of Electro-Chemo-Mechanical Deformation of Li-ion Batteries Based On Porous Electrode Theory

Lithium-ion batteries change their geometric dimensions during cycling as a macroscopic result of a series of microscale mechanisms, including but not limited to diffusion-induced expansion/shrinkage, gas evolution, growth of solid-electrolyte interphase, and particle cracking. Predicting the nonlinear dimensional changes with mathematical models is critical to the lifetime prediction, health management, and non-destructive assessment of batteries. In this study, we present an approach to implement an elastoplasticity model for powder materials into the porous electrode theory (PET). By decomposing the overall deformation into elastic, plastic, and diffusion-induced portions and using the powder plasticity model to describe the plastic portion, the model can capture the reversible thickness change caused by Li-ion (de-)intercalation as well as the irreversible thickness change due to the rearrangement and consolidation of particles. For real-world applications of the model to predict battery health and safety, the key lies in solving the mathematical equations rapidly. Here, we implemented the coupled model into the open-source software PETLION for millisecond-scale simulation. The computational model is parameterized using values gathered from literature, tested under varying conditions, brie ﬂ y compared to real-world observations, and qualitatively analyzed to ﬁ nd parameter-output relations. ©2024The

Lithium-ion batteries (LIBs) are widely used in electronics and road transportation.Over the past decades, extensive research on material innovation and engineering has led to many mature chemical systems that have been successfully commercialized.One of the important unsettled challenges facing the battery community is to reduce battery aging and its effects on safety and reliability.][3][4] Experimental evidence has shown that LIBs change their geometric dimensions during charge/discharge, referred to as "cyclinginduced deformation" hereinafter.6][7][8] Since some early experimental studies on the macroscopic measurement of battery cells, 9 many publications have investigated this important topic with microscopic tests and non-destructive assessment techniques. 10It is understood that the total change of dimensions consists of two parts, namely reversible and irreversible deformation. 11,12The former is fully recovered after a charge-discharge cycle and mainly comprises the deformation caused by the insertion and extraction of Li-ions in the solid active materials as well as the thermally induced deformation before any fracture.The latter is the unrecoverable deformation during cycling, which usually accumulates and becomes significant as the cycle number increases.The main sources of irreversible dimensional change include side reactions (e.g., Li plating, the formation of solid-electrolyte interphase (SEI), gas evolution), plastic mechanical deformation, and mechanical failure (e.g., cracks, fractures, delamination). 10,13Compared with the reversible part, the irreversible dimensional change is a more crucial factor of battery health because it is deeply two-way coupled with various aging mechanisms.In other words, aging induces deformation, and deformation, in turn, aggravates aging.
Like many engineering problems, modeling the cycling-induced deformation of battery cells is challenging due to the trade-off between accuracy and efficiency.For accuracy, it is preferable to implement the physical mechanisms at the micro-scale particle level.The finite element method (FEM) is often used, which requires a daunting computation cost due to the small scale of interest and the large number of elements (usually at the scale of millions). 5Thus, FEM cannot be applied for optimization-focused or control-oriented purposes.The discrete element method (DEM) is another useful tool for particle-level modeling. 14,15DEM treats particles as rigid spheres, which greatly improves computational efficiency.The problem with DEM-based electro-chemo-mechanical models is that it is difficult to implement electrochemical processes as the spheres are already greatly simplified for contact mechanics.To combine efficiency and accuracy, a practical approach is to develop an electrode-level model based on the Newman-type pseudo-2-dimensional (P2D) porous electrode theory (PET) model. 16,17Many studies have simulated the battery expansion/shrinkage by introducing a linear elasticity model with the diffusion-induced strain into the P2D model and have successfully predicted the expansion/shrinking over a few cycles. 8he fundamental cycling-induced deformation mechanisms during multiple cycles, however, have not been fully understood or modeled.
In real-world applications, the most critical requirement for computational efficiency of battery models is usually from battery management systems (BMSs). 3,18Recently, Berliner et al. 19 developed open-source software-PETLION-for millisecond-scale porous electrode theory-based lithium-ion battery simulations in the Julia language.A typical runtime for a dynamic simulation of full charge or discharge is only three milliseconds on a laptop.New features describing the aging effects such as the growth of SEI have been implemented into PETLION for lifetime prediction and management of batteries. 20Other widely used computational platforms include LIONSIMBA, 21 PyBaMM, 22 DUALFOIL, 23 and COMSOL MultiPhysics.
So far, these platforms have not been enabled to characterize cycling-induced deformation or swelling because the complex interplay between electrochemistry and mechanics has not been fully understood.There were promising computational attempts to tackle this problem.Garrick, Weidner, and co-workers [24][25][26] developed a methodology of implementing rock plasticity into PET to predict the volume change of porous electrodes during intercalation, by effectively introducing a compressibility coefficient.Thus, the complete traditional rock plasticity theory, which has a set of equations governing yield, hardening, plastic flow, and other effects, is greatly simplified.This inspiring work is worth more attention from the general battery community to further the electro-chemomechanical modeling of porous electrodes.For example, the compressibility coefficient that the authors developed is hard to determine experimentally or theoretically, making the method difficult to use.Furthermore, by nature, battery porous electrodes are closer to powders rather than rocks.Modeling porous electrodes using powder plasticity 14,[27][28][29][30] would enable the considerations of softer and more ductile particles, larger porosity, and the change of porosity during compaction (e.g., calendering an electrode).On the powder plasticity aspect, existing knowledge in the pharmaceutical industry of mechanically making pills 31,32 can be potentially used.In this study, we combine material-level theories of diffusion-induced volume strain with an electrode-level, elastoplastic continuum model, coupled tightly with a PET model for a continuous feedback loop.For this, a powder plasticity approach characterized by the Drucker-Prager yield criterion and flow potential is used.The overall computational framework is illustrated in Fig. 1.

Theory
Mass conservation, dimensional change, and porosity change.-Tomodel the deformation of the electrode, we start with the conservation of active materials in a control volume Ω t of the porous electrode, where e is the porosity of the electrode, and R i and V i are the reaction rate and molar volume of species i, respectively.Expanding the lefthand-side (LHS) of Eq. 1 to partial derivatives yields the continuity equation of the solid phases, where u is the displacement field and u is the velocity field.As pointed out by Garrick et al., [24][25][26] the first term of Eq. 2 is the change of porosity over time, and the second term is the dimensional change.Neglecting the non-uniformity of the porosity in the porous electrode, the second term on the LHS can be simplified as where ϵ v is the volumetric strain of the electrode.The continuity equation lays the foundation of the overall electro-chemo-mechanical framework by bridging the volumetric deformation, change of porosity, and reaction rate.In the following subsections, we further elaborate on every aspect of the governing equations.
General analysis of the mechanical deformation.-Fora deformable domain Ω, the kinematics is given by where ϵ is the strain tensor.The Dirichlet boundary condition of Eq. 4 is described as = û u on Ω u .The force equilibrium condition without body force is given by where σ is the Cauchy stress tensor.The Neumann boundary condition of Eq. 5 is described as σ• = n t on Ω t .The total strain ϵ is decomposed into three components, i.e., the diffusion-induced strain (also known as the eigen strain) ϵ d , the elastic strain ϵ e , and the plastic strain ϵ p , Figure 2 shows the geometry of the porous electrode of interest in the coordinates {x 1 , x 2 , x 3 }.In the Voigt notation, the stress tensor could be expressed as T , where σ ( = ) i j ij are the normal components and σ ( ≠ ) i j ij are the shear components.The same notation rule also applies to the strain components.In porous electrodes, the active powder materials are in a special condition that we can greatly simplify the stress and strain components.The thickness of the coating (h ≈ 60-100 μm) is significantly smaller than the width of the electrode disk (W ≈ 1-100 cm).Moreover, the coatings are firmly attached to a layer of the metal current collector.So, the powders far from the edge of the electrode are in an ideally confined condition.It is, Journal of The Electrochemical Society, 2024 171 050557 therefore, safe to neglect the "in-plane" and the shear components of the strain, i.e., ϵ 11 = ϵ 22 = ϵ 12 = ϵ 23 = ϵ 13 = 0.The only strain component of interest is ϵ 33 in the thickness direction.In other words, the deformation is simplified to one-dimensional.On the other hand, the low shearing resistance nature of powders enables us to neglect the shear components of the Cauchy stress in this case, σ 12 = σ 23 = σ 13 = 0. Given the large W/h ratio of the layer and the "flowing" nature of the powders, it is also reasonable to assume that the stress components in the x 1 and x 2 directions are equal: σ 11 ≡ σ 22 .In the 1D deformation assumption, these two stress components only vary in the thickness direction x 3 , and the remaining non-zero stress component σ 33 is constant as −p c , where p c is the stack pressure applied externally.This stress tensor automatically satisfies Eq. 5 for the force equilibrium.Based on these simplifications, the strain decomposition equation (Eq.6) can be rewritten as The governing equation for the x 2 direction is the same as that for the x 1 direction, and we only list the x 1 equation hereafter.These two equations clarified the main tasks of this article.The goal is to predict the deformation in the thickness direction, ϵ 33 .To achieve this goal, we build the relation between the diffusioninduced strain and the concentration field in Diffusion-induced strain section.Furthermore, we develop an elastoplasticity constitutive model to determine the components of elastic and plastic deformation in Powder plasticity for the active coatings section.Furthermore, assuming isotropic linear elasticity, the constitutive equations for stress, are given by the relations: where E is the elastic modulus, and ν is Poissons ratio.The porouselectrode-level values for these two parameters were gathered from literature 14,27 and are shown in Table I.Considering that the binder material could also influence the elastic properties of the electrode even at small percentage of the binder loading, 33 we used experimentally measured elastic modulus (∼500 MPa) of the coatings from literature in the present study, which already considers the contribution of the binder; the elastic modulus of active particles is much higher (tens to hundreds GPa). 34ffusion-induced strain.-Thedimensional change of a battery cell is a combined result of the intercalation and de-intercalation of Li ions and several complex degradation mechanisms as described in Introduction Section.In this article, we focus on the dimensional change in less than 50 cycles after formation.During this stage, it is reasonable to assume that the major degradation mechanisms such as gas generation and Li plating are not dominant. 10,35,36Our hypothesis is that the dramatic deformation at this early stage of the battery life cycle is due to the intercalation and de-intercalation of Li ions, which causes the volumetric expansion and extraction of active particles.In the fully confined condition as shown in Fig. 2, the particle-level volumetric change leads to the elasto-plastic deformation at the electrode level, which is elaborated on in the following section.In the literature, the diffusion-induced strain (eigenstrain) of the active material is often related to the concentration c s through Vegard's law, 11,37 , where c ref is a reference strainfree concentration and β is the Vegard coefficient tensor.For isotropic materials, this expression can be written as where Ω is the partial molar volume of lithium, 38 and I is the identity tensor.The resulting diffusioninduced volumetric strain ϵ d,v = Ω(c s − c ref ).Here, we choose the initial condition as the reference c ref = c s (t = 0) and normalize the concentration with the max concentration c max .A linear relation between the state-of-charge (SOC) and the volume change is obtained, where Ω * is the linear coefficient in terms of SOC.This simplified equation effectively captures the SOC-dependent volume change for many electrode materials such as NMC.During charging, the graphite electrode with intercalation shows volume expansion while in NMC electrode, it shows an approximately linear trend of volume shrinkage of 3%-5% 39 which is illustrated in Fig. 3. NMC particles exhibit distinctly anisotropic deformation induced by diffusion, 39 which is captured by implementing different coefficients for  different directions into the Vegard tensor β. 40 In this study, we neglect the anisotropy to simplify the model.On the other hand, realworld experimental data has shown highly nonlinear dependencies on ion concentration for graphite anode 39 because of its complex multi-stage phase transition process, 41 as illustrated in Fig. 3. Physically modeling this trend requires an electrochemical model that includes phase transition. 42Here, this trend is represented by a simple polynomial fit, where Ω * ( = ) i 0, 1, 2, 3 i are coefficients that are fit to data.
Powder plasticity for the active coatings.-Theelastoplastic behavior of the electrodes is modeled using the Drucker-Prager/Cap (DPC) model, as identified for granular electrode coatings by Zhu et al. 27 The DPC model describes plasticity in powder materials using a yield surface and flow potential. 32Its yield surface consists of two parts, a straight shear failure (sliding) surface after Drucker and Prager and an elliptical cap reminiscent of the crushable foam model developed by Deshpande and Fleck 44,45 In plasticity theory, the concept of yield surface is used to graphically describe the status of stress.While the stress state is within the yield surface, the material is elastic, and no plastic deformation takes place.When the stress state reaches the yield surface, the material undergoes plastic deformation such that the stress state does not surpass the yield surface.The DPC yield surface employed in this work combines the sliding behavior of powder materials and the consolidation behavior observed in crushable foams. 28,30This yield surface is defined in a space of the properties hydrostatic pressure, (the first invariant of the stress tensor) and equivalent stress (the second invariant of the deviatoric stress tensor).The yield surface is mathematically described by the yield function φ(p, q).Its isoline at φ = 0 is the yield surface, while φ < 0 indicates the stress state is within the yield surface which indicates the elastic region.
When the stress state reaches the yield surface, it will go to the plastic zone and flow along the surface.The yield function is defined as where R, β, d and p a are parameters determining the shape of the yield function, as illustrated in Fig. 4. The friction angle β describes the steepness of the shear failure (sliding) surface and the cohesion d its intersection with the q-axis.R influences the eccentricity of the elliptical consolidation surface beyond its intersection point with the sliding surface, which is given by p a , the hardening pressure at which the sliding mechanism transitions to the consolidation mechanism.
The hardening pressure p a is defined as a function of the volumetric plastic strain ϵ p,v , which allows the elliptical yield surface cap to grow or shrink, representing the hardening behavior of compressed foams.Compressive volumetric plastic strain, which is accumulated on the elliptical consolidation surface, causes hardening, expanding the cap, whereas volumetric plastic dilation, exhibited on the sliding surface, causes softening, shrinking the cap. 31Zhu et al. 27 have identified the power law, for the hardening pressure as a response to increasing compressive plastic volumetric compaction through indentation tests, where A and n are material parameters, A is a pressure in Pa, and n > 1 is a unitless exponent.The volumetric plastic strain ϵ p,v appears with a negative sign, so that compressive plastic strain leads to hardening.It is reasonable to assume that p a ≠ 0 at the initial time of the model, requiring the inclusion of a term for its initial value, p a,0 .As strains represent a change relative to an initial length, it is difficult to define the point at which ϵ p,v = 0, therefore, the ϵ a,0 was added as an offset, which can be directly determined from p a,0 .The power law used in this work is, therefore, where p a < 0 was ensured using a min function.
The parameters of the DPC yield surfaces for NMC and graphite used in this work were identified by Zhu et al. 27 through indentation tests and are shown in Fig. 5. Similarly, the hardening law was parameterized to roughly resemble the experimental relation between volumetric plastic strain and hardening pressure identified in the same work. 27nother important aspect of the plasticity model is the plastic flow potential G, which determines the direction of the plastic strain tensor, , . 19 Like the yield function, the flow potential was defined as a combination of the sliding behavior of powders for p < p a and consolidation behavior for p ⩾ p a .In particular, the right side of the flow potential has the same elliptical cap as the yield function following an associative flow rule, whereas a non-associative flow rule as used in the first case is common for powder materials. 27The above expression is equivalent to the DPC model in Abaqus without shear-induced dilatancy. 46Equation 18 is a tensor equation consisting of two directions of interest, i.e., x 1 and x 3 .To eliminate the additional unknown quantity ϵ | ¯̇| p , we define the ratio of plastic strain increments, ,11 The expression of Γ can be uniquely determined based on the choice of G.
For any condition, the six main equations derived above, i.e., (7-10, 15, 20), can adequately determine the six unknown variables {ϵ e,11 , ϵ e,33 , ϵ p,11 , ϵ p,33 , σ 11 , ϵ 33 }.In finite element-based computation, this set of equations is usually solved by the return mapping algorithm. 47At a time step, it is first assumed that the total strain increment is purely elastic to update the stress components using equations 9 and 10 as "trial stresses."The trial stress components are then plugged into the yield function Eq. 15 for checking the yield condition.If φ < 0, the step is elastic and the trial stress tensor should be used as the final stress.Otherwise, the step is plastic.In the case of φ > 0, the values of plastic strain should be increased until φ = 0.This process is commonly performed using the Newton-Raphson method, a widely used robust algorithm.The goal of this study is to implement powder plasticity into P2D model and solve the combined electro-chemo-mechanical model rapidly, at the millisecond level.The method of solving this set of mechanical equations is introduced in Computational Methods section.
Electrochemical model.-There is vast literature on coupled electrochemical transport models for porous electrodes.Typical modeling efforts leverage the separation between length scales at the particle level and the electrode level, leading to a variety of P2D models with one spatial dimension (x) along the electrode thickness and another within active particles (r) of the porous electrode.Early models include the Doyle-Fuller-Newman (DFN) 16 model and single-particle model (SPM). 48The DFN model is widely used in studies of the analysis, control, and management of battery systems.The DFN model is a system of nonlinear partial differential equations, which are more expensive to solve than the SPM which is much faster owing to major simplifications like quasi-steady fast diffusion in particles and negligible electrolyte potential gradient.More recently, multiphase porous electrode theory (MPET) by Smith and Bazant in 2017 42 was shown to capture essential effects arising out of phase separation phenomena in a distribution of particle sizes.
In this work, we extend the electrochemical P2D model in Torchio et al. 21This model is chosen since the electrochemical part has been implemented in PETLION for comparison against LIONSIMBA, another FVM-based porous electrode simulation framework on MATLAB.Additionally, a simple P2D model is appropriate for the purpose of early demonstration of powder plasticity in the P2D framework.The electrochemical model is characterized by the system of electrochemical reaction and transport equations summarized in Table II.The diffusion equation in the solid phase is modified by introducing a term accounting for temporal variation in the porosity.This modified P2D model is numerically solved in PETLION 19 which uses automatic differentiation (AD) for fast solution of the Partial Differential Algebraic Equation (PDAE) system (Table II) using a Finite Volume Method (FVM).For all simulations in this work, the same set of basic, electrochemical parameters from the NMC_LGM50 default dataset in PETLION 19 were used.

Computational Methods
Porous electrode theory presents a coupled system of partial differential algebraic equations.The finite element method (FEM) and finite volume method (FVM) are two commonly used routes for spatial discretization of the governing equation which results in a system of differential-algebraic system of equations (DAEs) that govern the time evolution of the vector of states, which are differential states if a time derivative appears in the DAE system and algebraic if entirely defined algebraically.PETLION uses the FVM for spatial discretization owing to its flux-conserving properties and ease of implementation/modification.The states in the resulting DAE system are implemented using the symbolics.jlJulia package whereas the equations themselves are represented by residuals that are implemented within ModelingToolkit.jlJournal of The Electrochemical Society, 2024 171 050557 which is a high-performance symbolic numeric computation package.Initialization is performed using initial values for differential states and Newton's method to solve for the initial value of algebraic states.The IDA method within Sundials.jl is used to solve the DAE system and is supplied with its jacobian which is calculated either symbolically or using automatic differentiation.
In this work, the system of partial differential algebraic equations represented by the combined electro-chemo-mechanical model was solved within the computational framework of PETLION.Additional states and corresponding residuals were added for the newly introduced six mechanical variables {ϵ e,11 , ϵ e,33 , ϵ p,11 , ϵ p,33 , σ 11 , ϵ 33 } and the yield function φ.Active material fraction and filler fraction were re-defined as differential variables.To solve the elastoplastic model, we introduced an intermediate variable namely the elastoplastic ratio, , which was added as an algebraic state in the PETLION code.c Numerically, the value of X falls between zero and one, 0 ⩽ X ⩽ 1.
The residual linked to X is unique due to the switching between the elastic and plastic step discussed in Powder plasticity for the active coatings section Section.The implementation of this switching was done by tracking the location of the stress state for each (i th ) volume element using a boolean (p i ) to indicate elastic (p i = 0) or plastic (p i = 1) deformation.Depending on the value of this boolean (p i ) the residual for X was continuously switched between ϕ ϵ ̇ḋ for elastoplastic deformation (p i = 1) and as X for purely elastic deformation (p i = 0).The boolean was dynamically updated using the continuous callback mechanism of PETLION with the conditions added to checks.jl.Aside from the φ > 0 condition for the elastic-to-plastic switch, additional conditions ensuring an increasing trend in φ were added to avoid back-and-forth switching due to numerical inaccuracies.

Results
In this section, we exhibit the performance of the electro-chemomechanical model in PETLION and compare with real-world data.
In practice, Li-ion batteries are cycled under conditions guided by a protocol chosen to limit undesirable degradation without making significant compromises in charging duration.The most common charging protocol is the Constant Current-Constant Voltage (CC-CV) protocol which involves first charging at a constant current until a cutoff voltage is reached and then a constant voltage hold till the current drops below a cutoff.The CC-CV protocol also serves as a baseline for comparison against newly designed protocols aimed at minimizing degradation and plating.Therefore, a CC-CV protocol was chosen to analyze the trends predicted by the simulation, as shown in Fig. 6.Evolution of the state of stress was explored by tracing the state of stress in Fig. 7 in the coordinates of equivalent stress versus pressure.The initial yield surface is plotted in solid black lines, and the red-and blue-colored lines stand for the stress trajectory during a stable cycle.Every point on the lines stands for the stress state at a certain moment.Note that the area enclosed by the yield surface and the pressure axis represents the elastic deformation.It is clear from Fig. 7a that the positive electrode (NMC in this study) only deformed elastically due to the high strength of the ceramic-type material.Graphically, the stress trajectory (red colored) never reached the yield surface.On the contrary, large plastic deformation was observed in the graphite negative electrode (Fig. 7b).During the intercalation, the stress reached the "cap" of the yield surface, indicating that the porous material got consolidated.The plasticity model that we used here allowed the material become stronger through the expansion of the cap surface.The dashed lines depict the post-yield surface.It got maximized at the end of the half-cycle of intercalation.When de-intercalation occurred, the stress decreased linearly in the elastic range and reached the p-axis, indicating a zero equivalent stress and a hydrostatic stress condition.Continued deintercalation causes further decreasing in pressure and growth of equivalent stress, as a reloading.The elastic reloading ended up with reaching the sliding branch of the yield surface, indicating that the material started plastic deformation through a particle frictional movement mode.Compared to the cap surface that governs the consolidation behavior, the sliding yield surface does not allow hardening.Graphically, the dots of stress state stayed on the sliding surface.The physical meaning behind it is that the particles moved   relatively and that there was a small frictional force contributing to the stress.

Stress
These two physical mechanisms, inner-particle consolidation and intra-particle sliding, have also been extensively investigated and confirmed for powder materials.Given the structural similarity between electrode coatings and granular powder materials, it is reasonable to expect that those mechanisms and plasticity models developed for powders should also be applicable to coatings.The main difference between the granular powder materials and the electrode coatings is the application of binder materials in the electrode coating, which could form an elastic matrix to constrain the particle movement and deformation.The binder connecting particles will undergo local shear and tension which could lead to inelastic deformation in the binder phase through fracture, delamination, or plasticity.Therefore, the binder would add the extra resistance to the particle sliding and consolidation during the plastic deformation.
From the results in Fig. 7, it is clear that our model can characterize the main deformation elastoplastic mechanisms of the electrodes, shedding light to predicting the diffusion-induced swelling.
Accumulation of plastic strain.-Figure8 shows the history of the strain components of the negative electrode (graphite) in the through-thickness direction during charge/discharge cycles.As the electrode is lithiated, a tensile diffusion-induced strain is opposed by a compressive elastic strain that accumulates until the consolidation part of the yield surface is reached after which, the additional diffusion-induced strain is opposed by lateral plastic strain.This is coupled with strain hardening that leads to enlargement of the yield surface.In the delithiation step, the diffusion-induced strain rate changes sign and this leads to an elastic unloading until the sliding part of the yield surface is reached after which, plastic strain changes and takes over till the next lithiation step in the next cycle.This trend agrees with the observation in Fig. 7.
The most important feature in Fig. 8 is the accumulation of the plastic strain component, ϵ p,33 .In the figure, the shaded ranges indicate the time during which the plastic strain increases.It is found that those ranges are correlated with the intercalation processes.As explained in Fig. 7, the accumulation of plastic deformation is caused by the consolidation of the electrodes, instead of the sliding.It is evident that our electro-chemo-mechanical model can capture well the deformation mechanisms of the porous electrode as expected.Moreover, the irreversible nature of the plastic strain component makes it possible to characterize the irreversible deformation, which is discussed in the next subsection.
Thickness change under different stack pressures and C rates.-Figure 9 shows the thickness change of the stack of two electrodes Δh relative to their total original thickness h 0 .The x-axis is the normalized time, i.e., the real time divided by the time needed for one cycle.It is clear that the model prediction of the thickness change consists of reversible and irreversible portions.Within one cycle, the reversible portion dominates, but over multiple cycles, the irreversible portion keeps accumulating and eventually exceeds the reversible portion.This trend agrees well with experimental observations reported by the literature 9 especially over the first 50-150 cycles before the overall thickness enters a plateau.Our current model is not able to capture the plateau and the long-range trend, which is a result of the complex interplay between mechanics and   Journal of The Electrochemical Society, 2024 171 050557 degradation mechanisms.We will elaborate on this point in the Discussions.
The external pressure p c is key mechanical boundary condition in the electrode model.As indicated by the force equilibrium equation in the thickness direction, Eq. 10, the σ 33 ) is constantly at the magnitude of −p c . Figure 9 also elucidates the impact of changing external pressure on the thickness change over time.Clearly, as the external pressure is increased from 2.25 MPa to 2.85 MPa, the thickness change reduces by about 100% after 20 cycles.This is in line with the physical picture of a higher external pressure hindering positive accumulation of plastic strain and can further be explained by the impact of external pressure on the elastoplastic ratio.
Except for the stack pressure, the C rate could also influence the electrochemo-mechanical interactions within the electrode.In addition to 1 C charging, simulations were performed for 2 C and 4 C charging rates while keeping our CC-CV protocol fixed.The relative thickness change curve shows a similar trend across different C-rates with time lags that are inversely dependent on the C-rate (Fig. 10a).In line with this, we also observe an earlier onset of plastic strain accumulation with increasing C-rate as shown in Fig. 10b.
Extension to considering long-term cycling.-Theaccumulated plastic strain in the electro-chemo-mechanical model is linked to a change in porosity through mass conservation.This is a contribution that would be neglected in the purely electrochemical porous electrode model.To elucidate this, the voltage response curves for a purely electrochemical model and our electro-chemo-mechanical model are shown in Fig. 11.The differences in the dashed and solid curves can be attributed to differences in the voltage polarization due to changes in transport properties of the cell.
In addtion, the plasticity model in this work can also be modified to perform long-term cycling simulations.One of the key long-term predictions would be the saturation value of the accumulated plastic strain in the system.To illustrate this, the thickness change of the cell was simulated again using a hardening law that has a stiffer (exponential) dependence in contrast to the power law in the main text, which is expressed as: The simulated relative thickness change using the standard protocol with 1C charging is shown in Fig. 12.Further comparisons were made to results in Willenberg 50 that reports thickness change over a single cycle for cylindrical Li-ion pure graphite anode Panasonic NCR Cells.Lacking the details about the cell information, it is difficult to compare the normalized thickness change.However, the nonlinear trend of the profile is observed to agree well with each other.In particular, the thickness change during charge and discharge is found to deviate from each other.This interesting feature is captured by our model, highlighting the importance of considering plasticity.

Discussion
Mechanisms of the irreversible thickness increase.-Asmentioned in the Introduction, there are multiple mechanisms responsible for the irreversible thickness change of a battery cell.Among them, the most widely accepted are the degradation-related ones, i.e., Li plating, formation and growth of SEI, cracking, and gas evolution.There is a popular understanding that irreversible thickness change is a direct result of degradation.Based on the experimental data reported in many publications, the irreversible thickness increase of a battery cell accumulates dramatically in the first 150-200 cycles. 9,49However, this understanding is questionable for shortduration cycling.In this work, we demonstrate that the plastic deformation of the porous electrode results in the irreversible deformation in the short-term cycling.A schematic illustration of the plastic deformation arised from the particle rearrangement and consolidation is shown in Fig. 15.
It should be very carefully stated that here the plastic deformation is defined on the electrode level, instead of the particle level.The irreversible thickness change predicted in this model is from the electrode scale, since the particles were treated as elastic, and no  degradation effect such as SEI growth or cracking was considered at the particle level.This is because the mechanical model of the electrode considered two different plastic deformation mechanisms, i.e., consolidation and sliding.The former refers to the change of shape of the particle, while the latter represents the relative motion or re-location of the particles.The good correlation in Fig. 13 suggested that our hypothesis about the effect of plastic deformation on the thickness change is promising.This was also supported by the non-linear prediction in Fig. 14, which showed that the thickness changes during the charge and discharge processes deviate from each other, even within one cycle.For a commercialized battery cell, the level of degradation in one early cycle is believed to be small enough to not influence the thickness change.
Although good correlations between the predictions and the experimental data are demonstrated, it is still necessary to perform more experimental and computational studies to validate our hypothesis about the irreversible deformation casued by particle consolidation and rearrangement as illustrated in Fig. 15.This paper is expected to serve as an invitation of this important topic.
Various mechanical boundary conditions.-Inthis study, we only considered one type of boundary condition (BC), which is a given external stack pressure on the electrode.In reality, this type of BC only exist in pouch batteries, which use a laminated architecture of the jellyroll in a bag.During manufacturing, the cell is vacuumized, leading to a boundary condition of p c = p atm .During charge-discharge, this pressure varies as the pouch material deforms, but the BC can still be safety treated as a pressure boundary.However, this type of BC is not valid for cylindrical cells or prismatic cells.Cylindrical cells has a wound structure of the jellyroll.This structure can withstand high internal pressures without significant deforming.Therefore, it is more reasonable to use a zerodisplacement BC, i.e., u 3 = 0. Prismatic cells have more complicated BCs.Their jellyrolls are encased in aluminum or steel with small initial gaps and allowed room for expansion.Therefore, a tractionfree BC is suggested before significant swelling, i.e., σ 33 | t=0 = 0, but when contact happens, the BC can only be determined by considering the deformation of the casing.
Our code can be easily modified for zero-displacement or traction-free BCs.For the more complicated contact BC, the properties of the casing must be obtained for developing a useful model.The assumption of 1D deformation.-Despitethe initial successful demonstrations, the 1D deformation assumption does not allow our model to capture the 3D deformation of the battery cell.According to Li et al., 49 the deformation of a pouch battery cell showed complex 3D features through DIC measurement.In order to characterize these features, the 1D deformation assumption must be reconsidered.First, the effect of the geometric edges must be included in the theory.The simplification of the stress and strain tensors in Fig. 2 is only valid for the materials far away from the edges and, therefore, cannot be used for 3D prediction.Second, the 1D deformation treatment relies on the assumption that the metal current collector can largely constrain the lateral deformation of the powders.According to Li et al., 49 the lateral deformation of a freestanding pouch battery cell is about 10% of the thickness change.In some cases, this lateral deformation should not be neglected.Lastly, it is clear that the stress and strain fields are significantly influenced by the material heterogeneity and the manufacturing process.Considering all these aspects, traditional finite element methods (FEM) is probably inevitable to solve the mechanical problem with complex geometries.As comparison, the clear advantage of the 1D deformation assumption is its simplicity to be incorporated into the PET theory.
Computational time.-Theelectro-chemo-mechanical theory developed in this study is general and, therefore, can be implemented into any computational platform to solve the PET equations.Here, we used the PETLION platform developed by Berliner et al. 19 for its advantage of high computational efficiency.The demonstrative cases in Fig. 9 took 100 ms.This short responding time enables the application of the model in advanced BMSs.To our best knowledge, existing BMSs are not designed to consider the cycling-induced deformation or swelling.Incorporating our PETLION-based electrochemo-mechanical model is expected to bridge this gap.In this study, we showed that the computational code can well predict the reversible and irreversible deformation during cycling.Particularly, it can capture the highly non-linear thickness change of a battery cell that is free from external loads.In real-world battery modules or packs, the cells are usually under constrained conditions, i.e., fixed space in enclosure structures or adjusted pressure by foams or gels.Given the design of the structure and boundary conditions, it is convenient to modify the model to determine the internal stress.Hence, our model will enable a BMS to perform pressure management of a battery module, which is of great importance for battery lifetime. 51

Conclusions
Battery cells are typical complex engineered systems.Understanding and predicting their behavior require a system-level characterization, including the electro-chemo-mechanical interaction.Characterizing battery swelling has become a challenge for the EV and electronics industry as it is directly related to the lifetime of the battery-reliant products.This study presents an early attempt at incorporating powder plasticity into P2D modeling of porous electrodes.Plastic deformation was described using an appropriate yield function and flow rule that combines the sliding behavior of powders and consolidation in foams.Parameters for the yield surfaces for NMC and graphite identified using indentation experiments when combined with the coupled electrochemical-mechanical model revealed that diffusion-induced stress in the negative electrode has a significant impact on cell thickness change.Comparisons were made against two experimental datasets wherein the trends in thickness change were found to be in good agreement with the experiments.
The main contribution of this model is that it for the first time revealed that plastic deformation of a porous electrode can lead to significant irreversible thickness change, even though no degradation mechanism was considered.This interesting finding sheds light on characterizing the dramatic thickness increase in the early cycles after formation, before degradation mechanisms become dominant.We also provided potential solutions to overcome the limitations of this model, including the lack of saturation mechanism in the accumulative plastic strain, the 1D deformation assumption, the reduced computational efficiency compared with the original PETLION platform, as well as the lack of degradation mechanisms and complicated boundary conditions.This model could be foundational in motivating the development of multi-scale electro-chemo-mechanical porous electrode models.At the same time, the model requires validation which can be made possible through dilatational experimental measurements to quantify cell-level displacement.Further analysis of spatio-temporal evolution of porosity using X-Ray CT could be conducted to validate the trends predicted by the model.

Figure 1 .
Figure 1.Overall flowchart of the computational framework of this study combining PET with powder plasticity.

Figure 2 .
Figure 2. Illustration of the stress and deformation states in the porous electrode.

Figure 3 .
Figure 3. Material-level strain as a function of lithiation ratio fitted to Shishvan et al. 43 .

Figure 4 .
Figure 4. Drucker-Prager cap yield surface and its parameters.
state of positive and negative electrode.-Thesteps in a single cycle of the CC-CV protocol are • Charge with I = 1C until V = 4.1 V • Hold V = 4.1V until I ⩽ 1/20 C • Rest for 30 min • Discharge with I = 0.5C until V = 2.7 V • Rest for 30 min

Figure 6 .
Figure 6.Simulation results of CC-CV cycles in the modified PETLION platform.

Figure 7 .
Figure 7. Yield surface and stress states at the center through the thickness of the (a) negative and (b) positive electrodes.

Figure 8 .
Figure 8. Vertical strains in the fifth negative electrode volume element.Shaded regions indicate plastic strain accumulation.Figure 9. Behavior at different levels of p c (material property p a = 2.75 MPa): relative thickness change over time normalized by period of one CC-CV cycle.

Figure 9 .
Figure 8. Vertical strains in the fifth negative electrode volume element.Shaded regions indicate plastic strain accumulation.Figure 9. Behavior at different levels of p c (material property p a = 2.75 MPa): relative thickness change over time normalized by period of one CC-CV cycle.
Comparison against realistic measurements.-Experimentsperformed on 20 Ah NMC batteries (around 100 layers) by Li et al. 49 using 3D Digital Image Correlation (DIC) are shown in Fig. 13.The increasing trend of relative thickness change for the 8 mm pouch cell is in remarkable agreement with the model predictions for a single cell.Again, note that the cycling tests performed by Li et al.only focused on the beginning of the life cycle before significant degradation happens.It therefore justifies the application of our model, which has not considered degradation effects.

Figure 11 .
Figure 11.Voltage curve under the standard CC-CV cycle with (dashed) and without (solid) plastic deformation enabled.

Figure 12 .
Figure 12.Relative thickness change for an exponential hardening law over 300 cycles.

Figure 13 .
Figure 13.Relative thickness change measurements by Li et al. 49 compared to unit cell model prediction plotted against time normalized by period of CC-CV cycle.

Figure 14 .
Figure 14.Thickness change over SOC for (a) dilatational Panasonic NCR experiment 50 and (b) normal simulation.

Figure 15 .
Figure 15.(a) shows the original position and shape of particles.After charging, in fig.(b), the particles will swell and rearrange with the change of shape.Later on, the reversible thickness change will vanish during discharging while the rearrangement and consolidation of particles will induce the irreversible thickness change.

Table I .
Electrode-level mechanical parameters.

Table II .
PDAE system corresponding to the Porous Electrode Model adapted from Ref. 21.Cathode/Anode (i ∈ {p, n})