A Modified Electrochemical Model to Account for Mechanical Effects Due to Lithium Intercalation and External Pressure

For a battery cell, both the porosity of the electrodes/separator and the transport distance of charged species can evolve due to mechanical deformation arising from either lithium intercalation-induced swelling and contraction of the active particles or externally applied mechanical loading. To describe accurately the coupling between mechanical deformation and the cell’s electrochemical response, we extend Newman’s DualFoil model to allow variable, non-uniform porosities in both electrodes and the separator, which are dynamically updated based on the electrochemical and mechanical states of the battery cell. In addition, the ﬁnite deformation theory from continuum mechanics is used to modify the electrochemical transport equations to account for the change of the charged species transport distance. The proposed coupled electrochemomechanical model is tested with a parameterized commercial cell. Our simulation results conﬁrm that mass conservation is satisﬁed with the new formulation. We further show that mechanical effects have a signiﬁcant impact on the cell’s electrochemical response at high

Lithium-ion batteries (LIBs) are complex systems where electrochemistry, mechanics, and thermal phenomena are closely coupled [1][2][3][4]. Mechanical deformation can impact a battery cell's electrochemical response in several ways. For example, the energy density of a battery cell can be significantly improved by selecting silicon or germanium as the anode material [5,6]. However, the diffusion-induced deformation is large enough to cause fracture of the active materials [7], leading to poor battery cyclic performance. For LIBs with conventional electrode materials, Li-intercalation causes the active materials to swell and contract, leading to a change of both the porosity of the electrodes and the species transport distance, which in turn affects the cell's electrochemical behavior [4,8]. In addition, mechanical deformation can in some circumstances increase the driving force for Li-plating [9]. And mechanical abuse can trigger short circuit, causing catastrophic failure of battery cells [10]. Thus, understanding the impact of mechanical deformation on the cell's electrochemical response can provide a basis for innovations that ensure battery safety, increase cycle life, and improve performance (higher energy density, faster charge/discharge rate).
Newman's physics-based porous electrode model (the so-called DualFoil model) [11][12][13][14] has laid the foundation for modeling battery systems. Various extensions have been made in the past to account for thermal effects [15,16] and other physical phenomena, such as side-reactions, film formation, and dendrite growth [17][18][19][20][21][22]. Though stress generation of active particles inside battery cells have been studied in [3,4,23,24], the development of complete, coupled electrochemomechanical models for full cells is still scarce. Many existing mechanics-related modeling works are carried out only on the particle level [4,[23][24][25][26][27][28][29]. In a recent work [30], the authors solve a full set of electro-chemothermo-mechanical coupled equations in the three-dimensional (3D) space with the finite element method, which accounts for the change of both porosity and species transport distance. However, their model is computationally expensive to solve and would not be fast enough for optimization or exploring a large parameter space .  1  2  3  4  5  6  7  8  9  10  11  12  13  14  15  16  17  18  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58 59 60 A c c e p t e d M a n u s c r i p t In this work, a coupled electrochemomechanical model is proposed for battery cells made with conventional electrode materials in the one-dimensional (1D) setting. During charge/discharge, both the porosity of the electrodes/separator and the species transport distance evolves because of Li intercalation and the externally applied pressure, as illustrated in Fig. 1. In the new model, the porosity of each porous region is dynamically updated based on the electrochemical and mechanical state variables. The finite deformation theory from continuum mechanics is used to describe accurately the change of charged species transport distance. In addition, a modified version of the 4th order polynomial approximation for computing the concentration profile in the solid active material [31] is provided to account for the diffusion distance change due to Li-intercalation. In the proposed model, we do not account for the electrolyte movement, thus the total volume of electrolyte is assumed to be constant at each mesh location. It is worth noting that such assumption is valid only when the level of the externally applied mechanical loading is low.

Model description
In the original DualFoil model [11][12][13][14], a constant porosity is assumed for each porous region in a battery cell. However, such assumption is not accurate and does not account for volume changes caused by Li intercalation or externally applied mechanical loading, which can result in non-uniform porosity along the sandwich layer through direction (SLTD), as illustrated in Fig. 1. To account for mechanical effects, we need to distinguish the undeformed configuration and deformed configuration, which are also called the reference configuration and the current configuration in the continuum mechanics theory. We use upper case letters to denote terms defined in the reference configuration, and lower case for these defined in the current configuration. In the following, we first summarize the existing DualFoil model with a constant porosity. Following that, we discuss how to account for the change of the porosity and the species transport distance with a new set of equations.  Figure 1: Illustration of a jellyroll sandwich layer in the 3D space (left) and the 1D modeling space for a discharge process (right). The letter Z is used to represent the sandwich layer through direction (SLTD), which points from one electrode to the other. In the right figure, the solid line and the dashed line represent the undeformed (reference) configuration along Z and the deformed (current) configuration along z, respectively. As the overall deformation in the X-Y plane is confined by the metal current collector and the cell case, we assume the averaged deformation only occurs in the Z direction with a fixed cross-section area A for both configurations 1 . During a discharge process, the volume contraction of graphite is larger than the volume swelling of LiCoO 2 , resulting in a total volume contraction as illustrated in the figure. The blue arrows indicate the modeling domain change due to external pressure, and red arrows for domain change due to Li insertion/removal. The actual length of the blue/red arrows for each porous region will be determined by the multiphysics model. Due to the volume change of the active particles, the negative electrode has an increased porosity, whereas the positive electrode has a decreased porosity. As shown in the deformed configuration, due to the spatial variation of Li concentration, the porosity across the electrode is non-uniform.

Existing electrochemical model
In conventional LIBs, two porous electrodes, which contain active energy storage materials, binders, and conductive additives, are separated by a porous separator, with the pores of all three regions being filled with liquid electrolyte 1 This assumption is supported by the strain study with three-dimensional digital image correlation in [32], where the strain is found to be negligible in the X and Y directions compared to the Z direction (the SLTD).  1  2  3  4  5  6  7  8  9  10  11  12  13  14  15  16  17  18  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59  60 A c c e p t e d M a n u s c r i p t

Equation description
Equation Table 1: Summary of equations to describe the electrochemistry in the DualFoil model with constant porosity [33]. [33]. During charge/discharge cycles, Li ions travel from one electrode to the other through the electrolyte, whereas electrons are conducted between the electrodes through an external circuit. In this process, electrochemical reactions occur at the electrode-electrolyte interfaces. The applied current as well as the distribution of electronic versus ionic current causes the electric potential to vary in both the liquid electrolyte phase and the solid phase. Heat is generated due to irreversible resistive heating, reversible entropic heating, and chemical reactions [15]. The DualFoil model and its extended versions for LIBs are capable of describing these physical processes and computing the voltage response, potential profiles, concentration profiles, and the heat generate rate by utilizing the porous electrode and concentrated solution theories. The related equations are summarized in Table 1, which are solved with a Crank-Nicolson method in the 1D setting.

Porosity change
Porosity ǫ e (or the volume fraction of the electrolyte) has a significant impact on the solution of equations listed in Table 1, as it is fully integrated into them through either its explicit appearance or effective quantities. The latter are computed based on the Bruggeman's relationship [34] (1) with α as the so-called Bruggeman's exponent.
Definition of porosity. To properly account for the porosity change and obtain an algebraic relation between the initial porosity ǫ 0 e and the current porosity ǫ e , we consider the volume change of different components in each porous region (i.e. the electrodes or the separator). In the reference configuration, ǫ 0 e is computed as for each region, where the total volume V tot equals to the sum of the partial volumes of the electrolyte V electrolyte , active materials V active , additional components V add , and the film V film produced by side reactions. Under the assumption that the liquid electrolyte is incompressible, similar to (2), the porosity ǫ e in the current configuration is computed as 3 Here, V denotes an initial volume defined in the reference configuration and v is for a volume defined in the current configuration. In the separator, we have V add = V film = v add = v film = 0.
Description of deformation. Since the reversible intercalation-induced volume change of a cell made with conventional electrode materials is relatively small (ca. 2.0% [35] and ca. 3.3% from our experiments), and only mild external mechanical loading is considered, the linear elastic theory is sufficient to describe the deformation and the volume change of each porous component. The displacement vector of each component of a cell is denoted as u(x, t), which is a function of the coordinate x and time t. The total strain tensor is defined as the symmetric part of the displacement gradient ε = 1 2 ∇u + (∇u) T (4) with the superscript T denoting a transpose operation. The ratio of the total volume change is expressed through the trace of the total strain tensor ε as which results in the following relation between the initial total volume V tot and the changed total volume v tot . The entire volume change of each region is assumed to arise either due to elastic deformation from an externally applied mechanical loading or due to an internal volume change induced by Li-intercalation in active materials or the film growth, which has an additive relation ∆v tot = ∆v elastic + ∆v active + ∆v film .
Thermal expansion induced volume change is neglected, as it is orders of magnitude smaller than these factors accounted in this work [36][37][38]. Within the scope of linear elastic theory, the strain tensor is expressed as where ε elastic , ε active , and ε film represent the elastic strain from external mechanical loading, the intercalation induced strain from active materials, and the film growth induced strain, respectively. The insertion of (7) into (5) leads to As the strain in the X and Y directions is found to be negligible compared to it in the SLTD [32], the strain tensor only has one non-zero diagonal term. Eqs. (5) and (10) can be simplified as where the superscript z denotes the SLTD, as illustrated in Fig. 1. In (11), ∆v film is obtained by accounting for different side-reactions [19][20][21] in the electrochemical model. Though the proposed formulations account for porosity change due to film growth, we neglect the side reactions in this work for simplicity.
Intercalation induced deformation of active materials. For the active materials, the intercalation induced volume change is much larger than the elastic deformation caused by externally applied mechanical loading. For example, the volume change by Li insertion and extraction is ca. 10% for graphite [35,39] and ca. 2% for LiCoO 2 [40], whereas for an externally applied high stack pressure of 5 MPa [9], the resulting volumetric strain is ca. 0.01% for graphite (with a bulk modulus of 35.8 GPa [41]) and ca. 0.002% for LiCoO 2 (with a bulk modulus of 205.6 GPa [42]). Thus, we can neglect the elastic deformation of the active material that results from the homogenized stress in the composite materials and express v active in the current configuration purely due to Li-interaction as where Ω(C) is the average relative volume change ratio of the active material. In (12), Li concentration C(X) is a location dependent quantity defined in the reference configuration, Ω(C) is the local relative volume change,V is the total volume of an active particle, and B is the domain describing the particle. In general, the intercalationinduced strain and Li concentration C has a non-linear relation [40,43], as shown in Fig 1  2  3  4  5  6  7  8  9  10  11  12  13  14  15  16  17  18  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59  60 A c c e p t e d M a n u s c r i p t SLTD undeformed shrink swell Figure 2: Illustration of the swelling/shrinkage of active particles in the 1D homogenized electrochemical model along the SLTD. Considering a representative 1D element with the given detailed microstructure in the undeformed configuration (only the active particle and the electrolyte are shown), the active particles could have different shapes and deform in any direction. Because the lateral deformation is confined by the case and metal current collector foils, the final deformation of the 1D element can only occur in the SLTD. As we do not account for the electrolyte movement across different elements, the total volume of the electrolyte is the same in the undeformed/shrink/swell configurations.
where ω is a swelling coefficient andC 0 is the averaged initial Li concentration in the active particle. By inserting (12) into (11), we have Since the value of Ω(C) depends on C, and the distribution of C along the SLTD, in general, is non-uniform, ε zz active is inhomogeneous along the SLTD.
Remark: In a battery cell, the active particles have different shapes (e.g. flake, spherical, etc.) that are not aligned along a specific direction, as shown in Fig. 2. The swelling/shrinkage of active particles could occur in any direction in the 3D space. When active particles deform in any direction, the particles will interact with the electrolyte. Since the lateral deformation is confined by the case and the metal current collector foils, the electrolyte can only flow in the SLTD. Thus, all the volume change of the active particles will result in a volume change of a porous region only in the SLTD due to the movement of the incompressible electrolyte. As we do not explicitly model electrolyte movement, the implication of the electrolyte movement in the neighborhood of a mesh point is thus reflected in Eq. (13) by stating that the volume change of the active particles causes a total volume change of a porous region only in the SLTD. Eq. (13) should not be interpreted as active particles are only allowed to deform in the SLTD.
Homogenized elastic deformation. The elastic strain tensor ε elastic describes the homogenized deformation of an entire composite layer (e.g., electrode or separator), rather than its individual components, because it is experimentally easier to measure the average deformation and mechanical properties of a porous region 2 . The non-zero component of the elastic strain tensor ε elastic is computed as where E and σ are the effective Young's modulus and the homogenized Cauchy stress for the composites in the SLTD, respectively. Each porous region has a different E, which needs to be measured by experiments. Thus, the resulting ε zz elastic is the same inside each region, but is different among regions. Updated porosity. Now, the porosity of electrode regions in the current configuration given in (3) can be rewritten as with ǫ 0 active and ǫ 0 add as the initial volume fraction for the active materials and the additive materials, respectively, which are computed as For the porous separator, which does not contain the active materials, Eq. (15) can be simplified as 2 For each porous region, as the electrolyte is incompressible, and the elastic deformation of the active materials caused by externally applied stress is negligible, the homogenized elastic strain εelastic is mainly attributed to the compressible additive materials in the electrodes and the solid phase in the separator, and potentially any unfilled pores in these regions.  1  2  3  4  5  6  7  8  9  10  11  12  13  14  15  16  17  18  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59  60 A c c e p t e d M a n u s c r i p t with ǫ 0 solid as the volume fraction of the solid phase materials. The strain tensor ε for the separator only contains the elastic component ε elastic with tr[ε film ] = 0 and tr[ε active ] = 0.

Modified electrochemical formulations
In addition to the porosity change, mechanical deformation impacts the solution of the equations listed in Table 1 via changing the domain over which these equations are solved, as illustrated in Fig. 1. To ensure the equivalence of equations that describe physics in both the reference configuration and the deformed configuration, we follow the continuum mechanics theory, which ensures consistency of physical properties of a continuum defined in arbitrary coordinate systems. Terms to account for the effect of deformation will appear in the relevant equations. We start with a finite strain formulation and define the mapping of physical properties to different configurations. To make the result more accessible to readers without a mechanics background, we apply the approximation of linear elastic theory to derive formulations that resemble those in Table 1.
Coordinate transformation. The finite strain formulation requires the differentiation between upper-and lower-case symbols, with the former describing quantities defined in the reference configuration and the latter those in the current configuration. The displacement u is defined as u = x − X, with x and X describing the current coordinate and the reference (initial) coordinate of the same part of a body, respectively. The deformation gradient F is defined as where 1 is the second-order identity tensor. For the problem considered in this work, where the in-plane strain is negligible and the cell behavior is described with the 1D model in the SLTD, the deformation gradient F has the form with e i and E J representing the basis vectors in the current and reference configurations, respectively. The term ∂u3 ∂X3 is equivalent to ε zz in (11). A volume dv defined in the current configuration is related to its equivalent representation dV in the reference configuration with where det[F ] is called Jacobian (equivalent to 1 + tr[ε] in (5) for the linear elastic theory) to describe the volume change ratio.
For a quantity defined per unit volume, such as concentration, we ensure the equivalent description by enforcing which leads to the relation c = C/J with c and C as the concentration in the current and reference configuration, respectively. The reference gradient ∇ X of a quantity (•) is transformed to its spatial gradient ∇ x by An area element da in the current configuration is related to its reference representation dA through To ensure that the rates of particles flowing through an area in two different configurations are identical, we have k · da = K · dA, resulting in k = J −1 F K (24) with k and K as the flux defined in current and reference configuration, respectively. Note that any flux variable (e.g., current density, which is a flux of charge) obeys the same relationship.
Electrolyte material balance. The electrolyte material balance in the reference configuration has the form where Q is the Li-ion flux due to the concentration gradient and I e is the current density in the reference configuration. The spatial flux q is computed as  1  2  3  4  5  6  7  8  9  10  11  12  13  14  15  16  17  18  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59  60 A c c e p t e d M a n u s c r i p t where the effective electrolyte diffusivity D e,eff is computed based on (1) with the updated porosity. Assuming the liquid electrolyte is incompressible, we have J e = 1, resulting in c e = C e . The readers should notice the difference between J e and J , where the former describes the volume change of the liquid electrolyte, and the latter is the total volume change of a porous region. Based on (22), (24), and (26), Q is computed as Inserting (27) into (25) yields the final form for the electrolyte mass balance equation as Electrolyte-phase Ohm's law. The electrolyte current density i e in the current configuration is computed as where the effective electrolyte conductivity κ e,eff is computed based on (1) with the updated porosity. Based on (22) and (24), the reference current density I e is computed as Solid-phase Ohm's law. The solid phase current density i s in the current configuration is computed as Based on (22) and (24), we obtain the reference current density Butler-Volmer insertion kinetics. For a spherical particle, we use j and J to denote the Li-ion surface fluxes in the current and reference configuration, n and N to denote the surface normal associated to j and J , respectively. The normal flux of Li ions on a particle's surface in the current configuration is expressed as Under the assumption of isotropic deformation, we havē for a spherical particle, whereF active is the deformation gradient to represent the overall deformation of the whole particle rather than any specific location, andJ active is the averaged volume change ratio of the particle withJ active = 1 + Ω(C). Based on (24), the reference flux J is computed as The normal flux of Li ions on a particle's surface in the reference configuration is thus expressed as With c s = C s /J active , c e = C e , and α a = 1 − α c , inserting (33) into (36) yields 3 Charge conservation. The charge balance equation in the reference configuration is expressed as where a 0 = 3ǫ 0 s /R 0 p [44]. 3 Note that this relationship is more complex when we relax the assumption that αa + αc = 1.  1  2  3  4  5  6  7  8  9  10  11  12  13  14  15  16  17  18  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59  60 A c c e p t e d M a n u s c r i p t Mass balance in active particles. In the reference configuration, the mass diffusion in solids is governed by where Q is the molar flux, which in general is composition and pressure dependent [4]. As shown in [4], the composition distribution is not significantly influenced by pressure except at high dimensionless current, therefore we neglect the pressure dependence of Q in the present model. We only consider a composition dependent Q governed by (in the current configuration) Based on (22) and (24), Eq. (39) becomes where F active is the deformation gradient at a specific location inside the particle, and J active is the corresponding volume change ratio at that location, which are location dependent. We assume that the active particle is disordered with lithium diffusion and lattice expansion/contraction being isotropic [4,23]. Eq. (41) becomes with F active = J active 1/3 1 . Eq. (42) can be solved by one extra discretization in the particle radial direction in addition to the main discretization in the SLTD (so-called pseudo 2D approach). Alternatively, one can utilize the three-parameter polynomial approximation presented in [31], where no additional discretization is needed, and the computational cost can be reduced significantly. Note that the formulation in [31] must be modified further, as in [45], to handle concentration-dependent diffusivities. For simplicity, the present model is restricted to concentrationindependent solid-phase diffusivity. Even in this case, the formulation in [31] is not directly applicable to (42), because J active is a radial location dependent quantity. Thus, to facilitate the computationally efficient polynomial approach, we further approximate J active with its volume averaged valueJ active , which has no location dependence, to simplify or with the form written in the spherical coordinate system as Eq. (44) is solved with a slightly modified version (see (51)-(53) in Table 2) of the formulation presented in [31]. A comparison study of the surface concentration C s resulted from (42) and its simplified form (43) with finite element methods in a spherical particle is given in the results section to show the accuracy of (43). Our simulation results confirm that the total amount of Li in both electrode active materials and the electrolyte is conserved with the new formulations.
1D and small deformation simplification. In the 1D setting under the assumption of small deformation, we have which lead to the simplification shown in Eqs. (46)-(50) of Table 2.
The proposed coupled electrochemomechanical model is solved with a Crank-Nicolson method in the 1D setting. The same BCs as used in [12,13] are applied to Eqs. The homogenized stress σ can be obtained by solving the linear momentum balance equation with an appropriate numerical tool, such as finite element methods, which is out of the scope of this work. In the results section, the externally applied stress value is pre-defined as a time-invariant boundary condition. To fully utilize the capability of the proposed model, a coupling between the electrochemical model with a mechanical or thermomechanical solver is needed to dynamically update σ based on the cell electrochemical states, which will be the subject of a subsequent paper [46].

Experiments and model parameterization
We measured the electrochemical and mechanical behavior of several LG G5 cell phone batteries, purchased from the official website [71] with model number BL-42D1F. After removing the control circuit, the cells were mounted into a customized platform, as shown in Fig. 3. The structure of the platform was designed by following [72] with two highprecision displacement sensors (1 µm accuracy and 0.1 µm resolution, model number Keyence GT2-H12KL, Japan) mounted in opposite directions to measure the amount of swelling during charge/discharge. The entire platform with the specimen was placed into a Yamato convection oven (Model DKN400) under forced convective cooling at a temperature of 25 • C. Several cells were charged/discharged at various C-rates (C/10, C/5, C/3, C/2, C/1) for three times at each C-rate, while their voltage and thickness changes were measured. For each C-rate, the cells were discharged to 2.7 V. Then, the cells were rest for 10 mins before charge. At the end of each charge, to 4.4 V, the voltage was held at 4.4 V until the current dropped to a C-rate of C/20. The cells were rest for 10 mins before discharge.  1  2  3  4  5  6  7  8  9  10  11  12  13  14  15  16  17  18  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59  60 A c c e p t e d M a n u s c r i p t  Table 3: List of simulation parameters for electrodes and separator. The bounds of fitting parameters are chosen based either fitted or measured values given in the references.
to measure its design parameters, such as electrode thickness, separator thickness, and separator area. A number of Li/LiCoO 2 and Li/graphite coin cells were made with electrode materials harvested from fresh LGG5 cells to measure the electrode open-circuit potentials (OCPs), as shown in Fig. 4. Interested readers are directed to Ref [62] for a comparison of OCP measurements from the Galvanostatic intermittent titration technique (GITT) and the averaged voltage method for a LCO/graphite cell that is similar to the one studied in this work. The functional description of the electrolyte from [33] was used for our simulation, with their formulations/values summarized in Table 4. The mechanical properties for each porous region were chosen from [57]. For graphite, a measured non-linear volume change ratio for a different pouch cell, as shown in Fig. 5, is used for the simulation, where Ω(C) is computed based on (12). A linear swelling ratio of ω = 1.8% is assumed for LiCoO 2 [40] with Ω(C) being computed based on Ω(C) = ω(C −C 0 ). The remaining parameters that could not be found in either of the sources mentioned, e.g., porosity, diffusion coefficient, reaction rate, etc., were determined by numerical optimization, i.e., by minimizing the L2-norm of the difference of the experimental and simulated voltage profile at given C-rates. The parameters used in the results section are summarized in Table 3 1  2  3  4  5  6  7  8  9  10  11  12  13  14  15  16  17  18  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59  60 A c c e p t e d M a n u s c r i p t    1  2  3  4  5  6  7  8  9  10  11  12  13  14  15  16  17  18  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59  60 A c c e p t e d M a n u s c r i p t  Figure 6: Comparison of surface concentration between the case neglecting mechanical effects and the case considering mechanical coupling, but using different formulations, (42) and (43). The surface potential associated with the surface concentration is shown on the right, where both plots show negligible difference among these three cases. In these simulations, the dimensionless time τ is defined as

Diffusion in active solid particles
The mass balance equations for active particles summarized in Table 2 are obtained based on (43), which is a simplified form of (42). In (43), an averaged volume change ratioJ active is used to solve for Li concentration instead of a radial location dependent volume change ratio J active . To evaluate the consequence on the surface concentration from this assumption, (42) and (43) are solved with finite element methods (FEM) in the two dimensional setting for a spherical particle. The solid diffusion in the spherical particle without mechanical deformation is also solved with FEM for comparison. A dimensionless analysis as carried out in [31] is conducted. A surface flux linearly dependent on the dimensionless pore wall flux defined in [31] is applied to the particle surface. The surface concentration and its associated surface potential for the three cases are plotted in Fig. 6(a,b), which show negligible difference for the three equations. Such results imply that the swelling the active particles has little impact on the Li diffusion in them. Despite of this observation, Eq. (43) is used for the remaining mechanically coupled simulations.

Discharge study of a commercial cell
In this section, we simulate the discharge process of LG G5 cells at C/1, C/2, C/3, C/5, and C/10 with the new formulations. Considering the fact that the cell has a very small Biot number [73], the C-rates are small, and the cells were tested under forced convection, an isothermal condition is assumed for all the simulations. The simulated thickness ratio change of the cell is compared with the measured value in Fig. 7, where a reasonably good match is observed. Because the swelling relations use for graphite and LiCoO 2 are not measured specifically for the tested LG G5 cells, the observed difference in Fig. 7 is expected. Comparison of cell voltage at different C-rates shows a good match between simulation results and experimental measurement, as plotted in Fig. 8. In fact, for smaller C-rates, mechanical deformation has little impact on the electrochemical response of the cell. Thus, the simulated discharge curves for these small C-rates from the DualFoil model without mechanical coupling will perform equally well. The mechanical effect on the cell's electrochemical response is more significant at higher C-rates, which is discussed in Section 4.3.

Porosity change in the SLTD
The porosity change of the three porous regions is plotted in Fig. 9, which shows that intercalation induced swelling/shrinkage causes the anode porosity to increase somewhat and the cathode porosity to decrease slightly.  1  2  3  4  5  6  7  8  9  10  11  12  13  14  15  16  17  18  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59  60 A c c e p t e d M a n u s c r i p t Porosity change impacts the effective transport properties of the cell via the Bruggeman relationship (1). The simulation results show that both electrodes have non-uniform porosity distribution during the simulation, which is due to the non-uniform distribution of Li intercalation/deintercalation rates in the electrodes. Our results show a maximum ca. 7.1% relative porosity increase in the anode region, and a maximum ca. 1.2% relative porosity decrease in the cathode region for this particular simulation. This porosity change should not be neglected in the simulation with higher Crates. For example, the discharge voltage profiles at C/3, C/1, and 1.2C from simulations with and without accounting for the mechanical effect are compared in Fig. 10. As expected, the difference is more significant at higher C-rates. This difference suggests that including mechanical effects enables more accurate simulations of the electrochemical behavior when extrapolating beyond the fitted operating parameter space to higher C-rates. Measurements of the LG G5 cells were restricted to 1C, and hence the model parameters may not be accurate at higher C-rates. We therefore limited our simulations to 1.2C, although this is not a general limitation of the model itself. We expand investigation of mechanical effects up to 4C in a subsequent work [74].

External mechanical loading effects
In a real cell, the jellyroll normally is wound layer by layer and is confined in the can. The constraint from the neighbor sandwich layer and the relatively rigid can will cause significant stress level change in the jellyroll. (c) cathode Figure 9: Illustration of the porosity change of the three porous regions in a cell during a C/1 discharge process at different time steps without externally applied pressure. A visible non-uniform porosity distribution in both electrodes is observed, which is caused by non-uniform distribution of the Li deintercalation reaction rate. understand how the externally applied pressure impacts cell electrochemical behavior, different normalized external pressures are applied to the cell at C/3, C/1, and 1.2C rates. The voltage profiles are shown in Fig. 11, which show that mechanical deformation has evident impact on cell electrochemical behavior at higher C-rates. This is mainly because that, at higher C-rates, Li-ion diffusion in the electrolyte is the major limiting factor, whereas porosity change can significantly impact the diffusion process via the Bruggeman relationship (1).

Conclusion
In this work, a new coupled electrochemomechanical model is presented to account for the porosity change and the transport distance change due to Li intercalation and externally applied pressure. Though carefully designed experiments are needed to thoroughly validate the proposed model, the preliminary results from a parameterized commercial LG G5 smartphone battery have shown promising performance, where the new model can accurately describe cell behavior at different discharge C-rates and reveal the potential impact of mechanical deformation on the cell's electrochemical behavior at high C-rates. The capability to model cell behavior with varying porosity is particularly important for simulating fast charge/discharge behavior, where the Li concentration variation (equivalently, the porosity variations) along the SLTD will be much more severe than the examples studied in this work. Because of the non-negligible coupling between mechanical deformation and electrochemical behavior, we expect that in large-format cells the non-  1  2  3  4  5  6  7  8  9  10  11  12  13  14  15  16  17  18  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59  60 A c c e p t e d M a n u s c r i p t (c) 1.2C Figure 11: Comparison of voltage profiles for C/3, C/1, and 1.2C discharge under different externally applied pressure, with the legend ∼0% − ∼3% indicating the equivalent level of externally applied elastic strain. The simulation results show that mechanical deformation has evident impact on cell electrochemical behavior at higher C-rates.