An electro-chemo-thermo-mechanical coupled three-dimensional computational framework for lithium-ion batteries

Thermal and mechanical effects play a vital role in determining the electrochemical behavior of lithium-ion batteries (LIBs). Non-uniform temperature distribution and mechanical deformation can result in uneven electrochemical states, leading to spatially varying aging rates that signiﬁcantly shorten cell lifetime. In order to improve simulation accuracy and thus the quality of computational battery design optimization, it is therefore essential to capture these coupled phenomena in a simulation model of a full battery cell. In this work, an electro-chemo-thermo-mechanical coupled framework is proposed to simulate LIBs in the three-dimensional space. In this new framework, a recently proposed one-dimensional electrochemical model, which includes the impact of mechanical deformation and local lithiation state on the effective transport properties of the charged species, is coupled with a three-dimensional thermomechanical model. A unique coupling scheme is proposed to handle information exchange between these two models. This framework allows us to accurately and efﬁciently study the behavior of three-dimensional cells with realistic geometry and resolve the spatial variation of interested ﬁelds. Two commercial cells are studied to show the performance of the newly proposed battery simulation framework. a uniform SOC and current density distribution in the jellyroll. (c-f) show that the newly proposed simulation framework is capable of capturing the spatially varying cell state quantities. The inhomogeneity of these quantities arises in (c,d) because a random initial porosity for the anode is used in the simulation. In (e,f), the inhomogeneity comes from the externally applied inhomogeneous mechanical loading.

A c c e p t e d M a n u s c r i p t  Figure 1: Illustration of two solution domains for the new 3D ECTM model, adapted from ref. [14]. (a) A fine finite element mesh of the jellyroll of an 18650 cell with n TM elements, on which the thermomechanical model is solved for the spatially varying temperature and stress fields. (b) A coarse mesh with each element spatially covering multiple thermomechanical elements. Each coarse element is associated with a 1D electrochemical sub-cell. The electrochemical behavior of the whole jellyroll is described by connecting all sub-cells, with a total number of n EC , in parallel with identical voltage. The sum of the current output from sub-cells is enforced to match the applied total current I with I = nEC l=1 i l . Volume-averaged quantities are used to exchange information between the thermomechanical model and the electrochemical model. in a single cell sandwich [6,7]. For larger cells or packs, sufficiently high applied current gives rise to nonuniform electrochemical states, temperature, and mechanical deformation, which in turn can result in spatial variations in the local current distribution, state of charge (SOC), or aging rate [13]. To investigate such spatial inhomogeneities, three-dimensional (3D) models are needed. Existing models to study thermal effects in the 3D setting [14][15][16][17][18] use a 1D DualFoil model to describe species transport and material balance, and a 3D model to calculate heat conduction. These two models are generally solved in a staggered manner based on the fact that heat conduction is several orders of magnitude faster than mass transport in the battery systems.
In lithium-ion batteries (LIBs), mechanical deformation can arise both from a volume change within the cell components, such as intercalation induced swelling of the active material, thermal expansion, film growth due to side reactions, and externally applied mechanical loadings [19][20][21][22]. Particularly, in battery cells that contain high energy density electrode materials, such as Silicon or Germanium, the large volume change of the active material during lithium (Li) intercalation can cause fracture or pulverization of the electrodes [23][24][25][26]. For these material systems, mechanical effects predominantly define the electrochemical performance and lifetime of the battery cells. For battery cells made with conventional electrode materials, such as graphite and LiCoO 2 , the volume change of the active material during Li intercalation is comparatively small [27][28][29]. However, it may still impact the electrochemical performance of a battery cell due to a change of the local porosity and because of an increase (decrease) of the distance that the charged species need to traverse as the electrodes expand (contract) [12,20,30]. These effects are typically cyclic and mostly reversible. However, irreversible processes, e.g., due to side-reactions such as Li-plating, may also occur and need to be considered in the definition of the optimal mechanical pressure applied to a battery cell [31].
In this work, we focus on cells made with conventional electrode materials. Compared to the number of 3D battery models that include thermal effects [14][15][16][17][32][33][34][35][36], battery models that consider the mechanical impact on the electrochemical performance of LIBs in the 3D dimensional setting are scarce. In [37], a one-way coupled electrochemo-thermo-mechanical (ECTM) model is proposed to study the stress level due to Li intercalation and thermal expansion for one sandwich layer, but the authors did not investigate how mechanical deformation affects the electrochemical behavior of LIBs. A 3D implementation of an extended DualFoil model including mechanical deformation, porosity change, and thermal expansion is presented in [38] and used to study coupling effects in a sandwich layer geometry. In a different work [39], the same authors investigate these effects by explicitly resolving the geometry of the pores and particles inside of the electrodes and separator. However, both models are computationally expensive and make it challenging to study the electrochemical behavior on the full-cell level.
To overcome these limitations, in this work, an efficient 3D ECTM coupled framework is proposed that allows studying the interplay of different physical processes occurring at multiple scales up to the cell level. For this purpose, we adapt the mesh mapping technique described in [14]. As illustrated in Fig. 1, a fine mesh is used to compute the homogenized thermomechanical behavior of the cell. This allows us to obtain a high spatial resolution for the temperature and displacement fields. To simulate the cell's electrochemical response, we generate a coarse mesh with each element representing a 1D electrochemical sub-cell. The sub-cell is discretized in the SLTD, which distinguishes the two 2  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 F o r R e v i e w O n l y electrodes and the separator. In our simulation, we do not explicitly create the coarse 3D mesh as shown in Fig. 1(b). Instead, a lookup table is generated to group multiple thermomechanical elements to represent the coarse mesh 1 . The overall electrochemical behavior of the jellyroll is obtained by connecting all the sub-cells in parallel, as shown in Fig.  1. In this framework, the sub-cells' parameters and states can be specified independently and updated in parallel on a high-performance computer. It allows us to assign different parameters (e.g., porosity and thickness) to each sub-cell to account for a non-uniform material distribution in the jellyroll. A two-way coupling scheme is used to exchange information between the thermomechanical and the electrochemical part of the simulation through volume-averaged quantities and projected stress.

Model Description
In this section, we describe the proposed coupled framework in detail. First, a 3D thermomechanical model formulated in the large deformation setting is presented. Next, a recently proposed mechanically coupled 1D DualFoil model [12] is summarized, which accounts for the change of the local porosity and the distance that charged species cross from one electrode to the other. We then describe the two-way information exchange scheme between the electrochemical model and the thermomechanical model.

Thermomechanical model
The finite strain theory in continuum mechanics is used to characterize the homogenized mechanical response of a battery. Interested readers are directed to [40] for details. Using the finite strain theory allows us to accurately compute the SLTD in the deformed configuration (see Section 2.3) and thus the σ nn used in the 1D DualFoil presented in Section 2.2.
The deformation of a battery is denoted by ϕ with ϕ = X + u, where X is the position vector of a point inside the battery in the reference configuration and u is the displacement vector of that point. The deformation gradient F = ∂ X ϕ = 1 + ∂u ∂X is a quantity to characterize the shape change of a solid body with 1 as a second order identity tensor. We consider a multiplicative decomposition of F with where F elastic , F swelling , and F θ describe the elastic deformation, intercalation-induced volume change, and thermal expansion, respectively. The termω in (1) is a volume change ratio due to Li-intercalation. F swelling describes the anisotropic volume change occurring in the SLTD, as the strain is found to be negligible in the two in-plane directions compared to the SLTD with three-dimensional digital image correlation technique in [41]. The thermal expansion coefficients for the electrodes and the separator are on the order 10 −5 K -1 [42]. Considering that the operation temperature range of a battery is typically less than 60 K, thermal-induced volume change is smaller than 0.06%. This value is much smaller compared to the reversible intercalation-induced volume change of around 2.0% [27]. Based on the above comparison, we neglect thermal-induced expansion by setting F θ = 1 . This assumption is supported by [37]. The deformation and the temperature fields ϕ and θ are obtained by solving the balance of linear momentum and the heat equation together with the properly imposed Dirichlet and Neumann boundary conditions. In Eq. (2), P is the first Piola-Kirchhof stress with P = ∂ F W . For the strain energy density W , we assume a Neo-Hookean relation where µ and λ are elastic material parameters, C elastic is the right elastic Cauchy-Green tensor with C elastic = F T elastic F elastic , and J elastic is the determinant of F elastic , respectively. With (4), the Cauchy stress σ is computed as 1 The lookup table can be created either based on the spatial coordinates, where thermomechanical elements in close proximity are grouped together, or based on the state variables, where thermomechanical elements with similar temperature and projected stress are grouped together. The latter case allows us to dynamically increase or decrease the total number of electrochemical subcells based on the evolution of state variables during charge/discharge to improve the efficiency of the electrochemical simulation.  (3) are solved with an in-house finite element code based on the deal.II library [43]. Interested readers are directed to Ref. [44] for details of the finite element method.

Mechanically coupled 1D DualFoil model
The DualFoil model and its modified versions describe the charge-transfer kinetics at the electrode-electrolyte interfaces, mass transfer and potential variation in both the active solid electrode materials and the liquid electrolyte based on the porous electrode and concentrated solution theory. The governing equations of this model are well described in the literature [2][3][4]. In order to capture thermal effects, carefully determined temperature dependent material parameters and a model for the various sources of heat generation are required. A comprehensive description of the individual heat generation rates in an electrochemical cell is given in [7]. This work employs a mechanically coupled DualFoil model [12], which accounts for the change of porosity and distance of charged species transport due to Li intercalation and externally applied mechanical loading. In the following, we give a short outline on how these effects are modelled. A summary of the resulting equations is given in Table 1.
Mechanical coupling. The effective electrochemical properties in Table 1, such as diffusivity and conductivity, depend on the porosity of the electrodes. For this, we assume the common Bruggeman's relationship [3] D e,eff = D e ǫ α e , κ e,eff = κ e ǫ α e , κ s,eff = κ s ǫ α s , where α is the Bruggeman exponent. As discussed in [12], the electrolyte volume fraction (or porosity) ǫ e in the two electrodes is updated as where ε nn is the total mechanical strain 2 in the SLTD, ǫ 0 add is the initial volume fraction of all components other than active material, e.g., binder or carbon black, ǫ 0 active is the initial volume fraction of the active material, and ǫ 0 e is the initial volume fraction of the electrolyte. The porosity of the separator is updated as with ǫ 0 solid as the initial volume fraction of the solid phase separator material. In this work, we do not consider film growth due to side reactions. However, the presented framework is capable of accounting for it. In (7) and (8), the strain field ε nn is computed as where ε nn elastic and ε nn active describe the strain due to elastic deformation and Li-intercalation, respectively. In (9), E is the effective Young's modulus, whose value is different for each porous region. σ nn is the projected stress in the SLTD computed from the thermomechanical model, which is the same for all the three porous regions. Thus, the value of ε nn elastic is different for each porous region. For the electrodes, the Li insertion/extraction induced volume change ratios are different. ε nn active is inhomogeneous in both electrodes along the SLTD due to the nonuniform Li distribution during charge/discharge. For the separator, ε nn active is zero. Thus, ǫ e in the electrodes is location dependent and varying along the SLTD, whereas ǫ e in the separator is homogeneous. A non-zero ε nn leads to an elongation or contraction of the porous composite materials and thus to a change of the distance that the charged species need to travel. This effect is reflected in the equations summarized in in Table 1.
Thermal coupling. The DualFoil model and the temperature field are coupled in two ways. First, the electrochemical properties, such as conductivity, diffusivity, and electrode kinetics in the equations in Table 1, follow an Arrhenius-type temperature dependence in the form [10] where (•) represents D s , D e , κ s , and k 0 , and (•) ref represents their corresponding reference values. Some of the electrolyte transport properties have a different dependence on temperature (see Table 4). Second, the local volumetric 2 The DualFoil model [12] is developed for battery cells made of conventional active materials with small volume change during charge/discharge. In small deformation theory, the mechanical strain tensor ε is defined as the symmetric part of the displacement gradient with ε = 1 2 ∇u + (∇u) T . Because the DualFoil model is a 1D model, the component of ε in the SLTD, which is a scalar quantity, is used instead of the full three-dimensional tensor ε.  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
Electrolyte material balance Electrolyte-phase Ohm's law Solid-phase Ohm's law Charge conservation Butler-Volmer insertion kinetics Intercalate boundary condition Table 1: Summary of the equations to describe the electrochemistry in the mechanically coupled DualFoil model [12] with the averaged volume change ratio of an active particle calculated asJ active = 1 + ε nn active /ǫ 0 active .
heat generation rate is computed by the DualFoil model based on [6] which accounts for irreversible heat generation due to electrochemical processes, Joule heating due to ohmic losses, and reversible entropic heat. As the heat generation rate is inhomogeneous along the SLTD, an integration over the thickness of the sandwich layer is performed in (19) to compute the heat generation rate per unit area, which is then divided by the initial thickness L 0 of the sub-cell to compute the volume-averaged heat generation rateq.

Coupling between electrochemical and thermomechanical models
Following [14], a fine mesh with n TM elements is used for the thermomechanical model, which allows us to obtain a high spatial resolution for the temperature and displacement fields. A coarse mesh with n EC elements is created. Each coarse element represents a 1D electrochemical sub-cell, which spatially covers multiple thermomechanical elements, as shown in Fig. 1. The 1D electrochemical sub-cell is discretized in the SLTD, along which the electrochemical model is solved. Note that the coarse element may be comprised of multiple and/or partial cell sandwich layers, and the 1D discretization and direction corresponds to the average of the corresponding layers. Due to the high conductivity of the current collectors, the electrical potential between two electrodes at any location of the jellyroll is approximated as identical, although the approximation introduces some error at high C rates in cells without multiple or continuous tabs [45]. Thus, the overall electrochemical behavior of the jellyroll can be modeled by connecting n EC sub-cells in parallel, where all the sub-cells have an identical electrical potential, but in principle differing temperature, current density, stress level, and porosity. The governing equations are solved with a Crank-Nicolson method in the 1D setting for each sub-cell. In this framework, the n EC sub-cells can be specified independently from one other. It allows us to assign different initial parameters (e.g., thickness, porosity) to each sub-cell to account for a non-uniform material distribution in the jellyroll.
A two-way coupling scheme is used to exchange information between the thermomechanical and the electrochemical submodels through volume-averaged quantities and projected stresses, as shown in Fig. 1. Even though the electrochemical model is solved in the 1D setting, due to the complex physics and the large number of total sub-cells, the electrochemical simulation can still be expensive. In order to obtain an optimal compromise between simulation 5  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 accuracy and runtime, we can vary the ratio between the number of thermomechanical elements and the number of electrochemical sub-cells. To obtain high resolution of the electrochemical state variables, we can define an electrochemical sub-cell mesh that aligns with the layers of the jellyroll. But often this is not necessary, as the adjacent cell sandwich layers can have similar electrochemical properties, which allows us to use one coarse electrochemical element to cover multiple cell sandwich layers to improve the computational efficiency. In addition, the presented approach also allows us to assign more sub-cells in those regions that require a higher spatial resolution of the electrochemical state variables and fewer sub-cells in those regions that tend to have more uniform states.  Fig. 1, where a local coordinate system with basis vectors (n, s, t) is defined with n indicating the SLTD for each element, along which the mechanically coupled 1D DualFoil model is solved. (b) Coordinate system transformation between Cartesian system with basis vector of (e 1 , e 2 , e 3 ) and the local coordinate system (n, s, t). As indicated by the red arrows, n is changing with the jellyroll geometry. The input stress value to the electrochemical model is σ nn , which is the component value of σ in the n direction. To obtain σ nn , we first calculate σ ′ in an intermediate is rotated along e ′ 2 with an angle of φ to obtain σ in the n direction.
The mechanically coupled DualFoil model takes the averaged stress σ nn of the jellyroll in the SLTD to update the porosity in Eq (9). σ nn represents the stress component of the homogenized Cauchy stress tensor σ computed from the thermomechanical model in Section 2.1, in the direction of n = F −T N . N is the SLTD defined in the reference configuration and n describes the SLTD obtained in the actual configuration. As illustrated in Fig. 2, we use (e 1 , e 2 , e 3 ) to denote the basis vectors of the Cartesian coordinate system. For an arbitrary direction n, to obtain the value of σ nn , we need to first rotate σ along e 3 with an angle of δ to an intermediate coordinate system (e ′ 1 , e ′ 2 , e ′ 3 ) with a rotation tensor Q 3 defined as The corresponding stress tensor σ ′ expressed in the coordinate system We then rotate σ ′ along e ′ 2 with an angle of φ to the local coordinate system (n, s, t) with a rotation tensor Q 2 defined as The corresponding stress tensorσ expressed in the coordinate system (n, s, t) is computed asσ Since each coarse electrochemical element is associated with multiple (assuming a ratio of k) fine thermomechanical elements, as shown in Fig. 1, the volume-averaged stress σ nn for a specific electrochemical element is computed as where the integration is performed via the Gaussian quadrature rule with two Gauss points 3 being used in each direction. The same σ nn is assigned to each node of the 1D electrochemical mesh to compute the ε nn elastic defined in (9). The 3 Interested readers are directed to Ref. [44] for details of numerical integration on Gaussian quadrature.
Output: thermomechanical model requiresω to compute the F swelling in (1), which is computed from the mechanically coupled DualFoil model asω where Ω(C) is the average relative volume change of the active material. The integral in (23) calculates the total volume change of the sub-cell due to Li insertion and extraction. 4 One can refer to [12] for a more detailed description. The sameω computed from one specific sub-cell is assigned as a Gauss point value to those thermomechanical elements that are associated with it to update the 3D stress and the displacement fields.

Averaged quantities for the thermal field
The DualFoil model takes a volume-averaged temperatureθ that is computed by the 3D thermomechanical model as an input parameter. Since each coarse electrochemical element is associated with multiple (assuming a ratio of k) fine thermomechanical elements, as shown in Fig. 1, the volume-averaged temperatureθ for a specific sub-cell is computed viaθ where the integration is performed via the Gaussian quadrature rule with two Gauss points being used in each direction. The sameθ is assigned to each node of the 1D electrochemical sub-cell. Thus, the temperature in each sub-cell is homogeneous, whereas it can be different across the electrochemical sub-cells. For each sub-cell, a volumetric heat generation rateq is calculated based on (19). The sameq is assigned as Gauss point values to those thermomechanical elements that are associated with the sub-cell to update the temperature field. This is transferred to all associated thermomechanical elements and used in the next iteration step to update the temperature field.

Simulation flowchart
The electrochemical model is coupled with the thermomechanical model via volume-averaged quantities, such as temperatureθ, heat generation rateq, intercalation-induced volume change ratioω, and mechanical stress σ nn . This information is exchanged between the 3D finite element code based on deal.II and the 1D DualFoil model code through 4 In a 1D electrochemical sub-cell, the area of the separator is assumed to be constant. Thus, the volume change of the sub-cell is equivalent to its thickness change in the SLTD.ω is a hypothetical thickness change that would occur assuming no external mechanical pressure on the sub-cell with incompressible electrolyte.  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 a coupling and electrical balance application programming interface (API). The API ensures that the sum of the current output from each individual sub-cell matches the applied current to the full jellyroll. In the simulation, this is achieved by iteratively updating the electrical potential applied to the sub-cells and checking the relative current difference. This API takes the normal stress σ nn and volume-averaged temperatureθ from the 3D thermomechanical model as inputs and passes the heat generation rateq and intercalation-induced volume change ratioω to the 3D thermomechanical elements. This exchange is performed in an alternating fashion. Hence, electrochemical model and thermomechanical model are solved in a staggered manner. The information flow of the framework is illustrated in Fig. 3 and summarized in Table 2.

Model Parameterization
We apply the model to investigate how mechanical deformation affects the electrochemical behavior of an 18650 cell and a prismatic cell. Both cells consist of a graphite anode, a LiCoO 2 cathode, a Celgard 2400 separator, and a 1M LiPF 6 salt in EC/DMC electrolyte. The parameterization of the simulation model corresponding to these cells is presented in detail in previous works [12,14] and is taken over in this study. For the cylindrical cell, the previous work [14] was carried out without accounting for mechanical deformation. Thus, we re-parameterize the cell for the mechanically coupled electrochemical model by combining reported values from the literature, measured values from experiments, and optimized values from simulations. Several of the parameters, e.g., diffusion coefficient, reaction rate, activation energy, etc., were re-tuned by numerical optimization, i.e., by minimizing the L2-norm of the difference of the original experimental voltage profile reported in [14] and the simulated voltage profile at given C-rates. The full set of material parameters used in the simulations are summarized in Tables 3 through 5. Following [46], the homogenized Young's modulus of the jellyroll in the SLTD in Table 5 is computed based on Voigt averaging via where i refers to the individual components in the jellyroll, e.g., the current collectors, electrodes, and the separator.  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 4: List of parameters and properties of the electrolyte in the 18650 cell [14].
The current density and SOC 6 distribution at the end of discharge of the three cases are shown in Fig. 5. Based on the specified discharge protocol, about 45% of the cell capacity is obtained for each case, resulting in the final value of the SOC around 0.55. For case (i), because of the uniformity of both internal and external states, a uniform SOC and current density distribution in the jellyroll is observed, as shown in Fig. 5(a,b). For case (ii) and (iii), due to either the non-uniform initial cell state or non-uniform external mechanical loading, the jellyroll has unbalanced internal electrochemical states. A spatially varying SOC and current density distribution is observed for both cases, as shown in Fig. 5(c-f). These results confirm the capability of the newly proposed framework to capture spatially varying quantities in a Li-ion battery.

Simulation of a cylindrical cell
Next, an 18650 cell is discharged at 4A rate under both natural (h = 12.5 W/(m 2 ·K)) and forced (h = 49.5 W/(m 2 ·K)) convection conditions, with h as the external heat transfer coefficient. The original experimental setup can be found in [14]. In the simulation, we use a fine mesh with 61,700 elements for the thermomechanical model, as shown in Fig.  7, with the jellyroll being divided into 90 sub-cells (3 in the radical direction, 6 in the azimuthal direction, and 5 in the height direction) for the electrochemical simulation (see Fig. 8).
In the simulation, the proper mechanical constraint at the interface between different parts of a battery cell is not easy to measure. In addition, the cell initial internal pressure might vary among vendors. Future experimental work is needed to determine if there is any relative movement between the jellyroll and other parts during charge/discharge, and what is the appropriate range of initial internal pressure level inside a battery cell. The presented model can 6 In this work, SOC is defined as the ratio between the remaining cell capacity and the total cell capacity between equilibrium voltages of 4.11 V and 2.8 V for the 18650 cell, and between 4.36 V and 2.8 V for the prismatic cell. During discharge, the SOC has an initial value of 1.0.  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 8: Illustration of the SOC and current density distribution in the coarse electrochemical sub-cells at the end of discharge under the forced convection condition. The non-uniformity of the SOC and current density is more correlated to the non-uniform temperature than the projected stress field in Fig. 7. be used to address these uncertainties and evaluate how different cell configurations will impact the electrochemical response of battery cells. In this simulation, we neglect the initial internal pressure applied to the jellyroll. We assume that all the parts are fully connected with no relative movement between each other during charge/discharge.
The simulated cell voltage and core temperature of the 18650 cell is compared with experimental measurements for both natural and forced convection conditions, as shown in Fig. 9. We observe that the simulated cell voltage profiles match the experimental discharge curves very well. The predicted core temperature differs from the measured values by less than 3 • C. High resolution non-uniform temperature and projected stress distributions from the simulation are shown in Fig. 7. Because of the aforementioned assumption, a high stress is observed at the bottom of the jellyroll, as shown in Fig. 7(b), which differs from other regions primarily due to the mechanical interaction between the jellyroll and other parts. Under the forced convection condition, about 35% of the cell capacity is discharged (see Fig. 9), resulting in the final value of the SOC around 0.65. The SOC and current density distribution in the coarse electrochemical sub-cells are shown in Fig. 8, from which we can observe that, in this particular simulation, the nonuniform SOC and current density is more correlated to the non-uniform temperature than the projected stress, where the SOC is lower in the region with higher temperature due to its lower electrochemical resistance during 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  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

Simulation of a prismatic cell
Lastly, a prismatic cell, which has been previously studied and parameterized in [12,[59][60][61], is investigated. For the sake of simplicity, we neglect the can of the battery cell and simulate only the jellyroll. In addition, because the jellyroll is symmetric with respect to x − y plane, y − z plane, and x − z plane, only 1/8th of it is simulated, as shown in Fig. 10. 511 elements for both the thermomechanical model and the electrochemical model are used. Furthermore, as pointed out in [59], this prismatic cell has a fairly uniform temperature distribution due to a very small Biot number. To better understand the impact on the electrochemical behavior of the jellyroll from the mechanical deformation, an isothermal condition is considered. The jellyroll is discharged at the C/2 rate at a temperature of 25 • C with an initial capacity of ∼ 2.7 Ah. Both the voltage profile and the thickness ratio change of the prismatic cell are shown in Fig. 11, where a very good match of the voltage profile between simulation and experimental data is observed. The simulation takes the experimentally measured swelling ratio vs SOC for both graphite and LiCoO 2 active particles as inputs and outputs the averaged thickness change of the whole jellyroll. As pointed out in [12], the observed difference in the thickness ratio change arises because the swelling ratio vs SOC for both graphite and LiCoO 2 used in the simulation is measured from a different cell.
To understand how mechanical deformation affects the electrochemical behavior of this prismatic cell, the distributions of the projected stress, SOC, and current density are shown in Fig. 12. From Fig. 12(b), one can observe that a high stress concentration rises in the simulation due to intercalation-induced volume change and the curved geometry of the jellyroll. This high stress concentration results in a non-uniform SOC and current density distribution. At the end of the discharge, a ∼ 8.5% difference is observed between the maximum and the minimum SOC values. The current density shows a ∼ 5.8% difference between the maximum and the minimum values. The unbalanced electrochemical states resulting from inhomogeneous mechanical deformation due to the jellyroll geometry may lead to non-uniform Li-plating during charging, as observed at the edge of prismatic cells [62,63]. This effect is the subject of a manuscript in preparation, in which the same model is applied to cell charging.

Conclusion
In this work, a coupled framework of an electro-chemo-thermo-mechanical model is proposed that is capable of efficiently and accurately simulating LIBs at the cell/pack level. A cylindrical cell and a prismatic cell are studied to demonstrate the performance of the presented framework. The simulation results obtained from this new framework reveal that localized stress concentration arises due to intercalation-induced volume change and jellyroll geometry. Such localized mechanical deformation can result in unbalanced electrochemical states that could potentially be related to non-uniform Li-plating. These results demonstrate the impact of including mechanical effects in a battery simulation. However, it is worth to point it out that the quantitative results due to mechanical deformation presented in this work should only be interpreted qualitatively. Future research concerning the accurate measurement of effective mechanical properties of each component in a Li-ion battery, the understanding of how to properly impose mechanical constraints inside the cell, and the careful experimental design to thoroughly measure the electrochemical, thermal, and mechanical parameters is needed to achieve a comprehensive understanding of the coupling between the thermomechanical effects and the electrochemical behavior of LIBs. In addition, for cells operating under high external mechanical loading and/or for cells designs with larger active material volume change (e.g. Si-based), a more sophisticated coupling scheme between mechanics and electrochemistry will be needed to account for the bulk electrolyte movement.   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 Poisson's ratio [-] ω