An Electro-chemo-thermo-mechanical Coupled Three-dimensional Computational Framework for Lithium-ion Batteries

The finite strain theory in continuum mechanics is used to characterize the homogenized mechanical response of a battery. Interested readers are directed to Ref. 40 for details. Using the finite strain theory allows us to accurately compute the SLTD in the deformed configuration and thus the σⁿⁿ used in the 1D DualFoil model.

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 ${\boldsymbol{F}}={\partial }_{{\boldsymbol{X}}}{\boldsymbol{\varphi }}={\bf{1}}+\tfrac{\partial {\boldsymbol{u}}}{\partial {\boldsymbol{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

$$\begin{eqnarray}&&{\boldsymbol{F}}={{\boldsymbol{F}}}_{\mathrm{elastic}}{{\boldsymbol{F}}}_{\mathrm{swelling}}{{\boldsymbol{F}}}_{\theta }\quad \mathrm{with}\quad {{\boldsymbol{F}}}_{\mathrm{swelling}}={\bf{1}}+\bar{\omega }{\boldsymbol{n}}\otimes {\boldsymbol{N}}\end{eqnarray} \tag{ 1 }$$

where F _elastic, F _swelling, and F _θ describe the elastic deformation, intercalation-induced volume change, and thermal expansion, respectively. The term $\bar{\omega }$ 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 Ref. 41. The thermal expansion coefficients for the electrodes and the separator are on the order 10⁻⁵ K⁻¹. ⁴² 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%. ²⁷ Based on the above comparison, we neglect thermal-induced expansion by setting F _θ = 1. This assumption is supported by Ref. 37. The deformation and the temperature fields φ and θ are obtained by solving the balance of linear momentum and the heat equation

$$\begin{eqnarray}&&\mathrm{Div}[{\boldsymbol{P}}]={\bf{0}}\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&\rho {C}_{p}\displaystyle \frac{\partial \theta }{\partial t}+{\rm{\nabla }}\cdot \left(-{\kappa }_{\theta }{\rm{\nabla }}\theta \right)=\dot{q}\end{eqnarray} \tag{ 3 }$$

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

$$\begin{eqnarray}&&W({\boldsymbol{F}})=\displaystyle \frac{\mu }{2}\left(\mathrm{tr}({{\boldsymbol{C}}}_{\mathrm{elastic}})-3-2\mathrm{ln}{{ \mathcal J }}_{\mathrm{elastic}}\right)+\displaystyle \frac{\lambda }{2}{\left({{ \mathcal J }}_{\mathrm{elastic}}-1\right)}^{2},\end{eqnarray} \tag{ 4 }$$

where μ and λ are elastic material parameters, C _elastic is the right elastic Cauchy-Green tensor with ${{\boldsymbol{C}}}_{\mathrm{elastic}}={{\boldsymbol{F}}}_{\mathrm{elastic}}^{T}{{\boldsymbol{F}}}_{\mathrm{elastic}}$, and ${{ \mathcal J }}_{\mathrm{elastic}}$ is the determinant of F _elastic, respectively. With 4, the Cauchy stress σ is computed as

$$\begin{eqnarray}&&{\boldsymbol{\sigma }}=\displaystyle \frac{{\boldsymbol{P}}{{\boldsymbol{F}}}^{T}}{{ \mathcal J }}=\displaystyle \frac{1}{1+\bar{\omega }}\left(\lambda ({{ \mathcal J }}_{\mathrm{elastic}}-1){\bf{1}}+\displaystyle \frac{\mu }{{{ \mathcal J }}_{\mathrm{elastic}}}\left({{\boldsymbol{F}}}_{\mathrm{elastic}}{{\boldsymbol{F}}}_{\mathrm{elastic}}^{T}-{\bf{1}}\right)\right).\end{eqnarray} \tag{ 5 }$$

In Eq. 3, $\dot{q}$ is calculated based on Eq. 19 from the electrochemical model. The Eqs. 2 and 3 are solved with an in-house finite element code based on the deal.II library. ⁴³ Interested readers are directed to Ref. 44 for details of the finite element method.

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–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 Ref. 7. This work employs a mechanically coupled DualFoil model, ¹² 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 I.

Table I. Summary of the equations to describe the electrochemistry in the mechanically coupled DualFoil model ¹² with the averaged volume change ratio of an active particle calculated as ${\bar{{ \mathcal J }}}_{\mathrm{active}}=1+{\varepsilon }_{\mathrm{active}}^{{nn}}/{\epsilon }_{\mathrm{active}}^{0}$.

Equation description	Equations
Electrolyte material balance	$$\begin{eqnarray}&&\displaystyle \frac{\partial ({\epsilon }_{e}^{0}{C}_{e})}{\partial t}={{\rm{\nabla }}}_{X}\cdot \left(\displaystyle \frac{{D}_{e,\mathrm{eff}}}{1+{\varepsilon }^{{nn}}}{{\rm{\nabla }}}_{X}\left({C}_{e}\right)\right)+{{\rm{\nabla }}}_{X}\cdot \left(\displaystyle \frac{{{\boldsymbol{I}}}_{e}(1-{t}_{+})}{F}\right)\end{eqnarray} \tag{ 10 }$$
Electrolyte-phase Ohm's law	$$\begin{eqnarray}&&{{\boldsymbol{I}}}_{e}=\displaystyle \frac{-{\kappa }_{e,\mathrm{eff}}{{\rm{\nabla }}}_{X}{\phi }_{e}+2\tfrac{{\kappa }_{e,\mathrm{eff}}R\theta }{F}(1-{t}_{+})\left(1+\tfrac{\partial \mathrm{ln}{f}_{\pm }}{\partial \mathrm{ln}{C}_{e}}\right){{\rm{\nabla }}}_{X}\mathrm{ln}{C}_{e}}{1+{\varepsilon }^{{nn}}}\end{eqnarray} \tag{ 11 }$$
Solid-phase Ohm's law	$$\begin{eqnarray}&&{{\boldsymbol{I}}}_{s}=\displaystyle \frac{1}{1+{\varepsilon }^{{nn}}}\left(-{\kappa }_{s,\mathrm{eff}}{{\rm{\nabla }}}_{X}{\phi }_{s}\right)\end{eqnarray} \tag{ 12 }$$
Charge conservation	$$\begin{eqnarray}&&{{\rm{\nabla }}}_{X}\cdot ({{\boldsymbol{I}}}_{s}+{{\boldsymbol{I}}}_{e})=0\,\mathrm{with}\,{{\rm{\nabla }}}_{X}\cdot {{\boldsymbol{I}}}_{e}=-{a}_{0}{{FJ}}_{n}\,\mathrm{where}\,{a}_{0}=\displaystyle \frac{3{\epsilon }_{s}^{0}}{{R}_{p}^{0}}\end{eqnarray} \tag{ 13 }$$
Butler- Volmer insertion kinetics	$$\begin{eqnarray}&&{J}_{n}=\displaystyle \frac{{k}_{0}{C}_{s}^{{\alpha }_{c}}{C}_{e}^{{\alpha }_{a}}{\left({C}_{s}^{\max }-{C}_{s}\right)}^{{\alpha }_{a}}\left[\exp \left(\tfrac{{\alpha }_{a}F}{R\theta }{\eta }_{s}\right)-\exp \left(-\tfrac{{\alpha }_{c}F}{R\theta }{\eta }_{s}\right)\right]}{{\left({\bar{{ \mathcal J }}}_{\mathrm{active}}\right)}^{1/3}}\end{eqnarray} \tag{ 14 }$$
Intercalate material balance	$$\begin{eqnarray}&&\displaystyle \frac{d}{{dt}}\bar{{C}_{s}}+3\displaystyle \frac{{J}_{n}}{{R}_{p}^{0}}=0\end{eqnarray} \tag{ 15 }$$
Volume-averaged flux relation	$$\begin{eqnarray}&&\displaystyle \frac{d}{{dt}}\bar{Q}+30\displaystyle \frac{{D}_{s}}{{\bar{{ \mathcal J }}}_{\mathrm{active}}^{2/3}{\left({R}_{p}^{0}\right)}^{2}}\bar{Q}+\displaystyle \frac{45}{2}\displaystyle \frac{{J}_{n}}{{\left({R}_{p}^{0}\right)}^{2}}=0\end{eqnarray} \tag{ 16 }$$
Intercalate boundary condition	$$\begin{eqnarray}&&35\displaystyle \frac{{D}_{s}}{{\bar{{ \mathcal J }}}_{\mathrm{active}}^{2/3}{R}_{p}}({C}_{s}-\bar{C})-8\displaystyle \frac{{D}_{s}}{{\bar{{ \mathcal J }}}_{\mathrm{active}}^{2/3}}\bar{Q}=-{J}_{n}\end{eqnarray} \tag{ 17 }$$

Mechanical coupling. The effective electrochemical properties in Table I, such as diffusivity and conductivity, depend on the porosity of the electrodes. For this, we assume the common Bruggeman's relationship ³

$$\begin{eqnarray}&&{D}_{e,\mathrm{eff}}={D}_{e}{\epsilon }_{e}^{\alpha },\quad {\kappa }_{e,\mathrm{eff}}={\kappa }_{e}{\epsilon }_{e}^{\alpha },\quad {\kappa }_{s,\mathrm{eff}}={\kappa }_{s}{\epsilon }_{s}^{\alpha },\end{eqnarray} \tag{ 6 }$$

where α is the Bruggeman exponent. As discussed in Ref. 12, the electrolyte volume fraction (or porosity) _e in the two electrodes is updated as

$$\begin{eqnarray}&&{\epsilon }_{e}=\displaystyle \frac{1-{\epsilon }_{\mathrm{add}}^{0}-{\epsilon }_{\mathrm{active}}^{0}}{1+{\varepsilon }^{{nn}}}=\displaystyle \frac{{\epsilon }_{{\rm{e}}}^{0}}{1+{\varepsilon }^{{nn}}},\end{eqnarray} \tag{ 7 }$$

where εⁿⁿ is the total mechanical strain ^b in the SLTD, ${\epsilon }_{\mathrm{add}}^{0}$ is the initial volume fraction of all components other than active material, e.g., binder or carbon black, ${\epsilon }_{\mathrm{active}}^{0}$ is the initial volume fraction of the active material, and ${\epsilon }_{{\rm{e}}}^{0}$ is the initial volume fraction of the electrolyte. The porosity of the separator is updated as

$$\begin{eqnarray}&&{\epsilon }_{e}=\displaystyle \frac{1-{\epsilon }_{\mathrm{solid}}^{0}}{1+{\varepsilon }^{{nn}}}=\displaystyle \frac{{\epsilon }_{{\rm{e}}}^{0}}{1+{\varepsilon }^{{nn}}}\end{eqnarray} \tag{ 8 }$$

with ${\epsilon }_{\mathrm{solid}}^{0}$ 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 εⁿⁿ is computed as

$$\begin{eqnarray}&&{\varepsilon }^{{nn}}={\varepsilon }_{\mathrm{elastic}}^{{nn}}+{\varepsilon }_{\mathrm{active}}^{{nn}}\quad \mathrm{with}\quad {{\boldsymbol{\varepsilon }}}_{\mathrm{elastic}}^{{nn}}=\displaystyle \frac{1}{E}\,{\sigma }^{{nn}}\end{eqnarray} \tag{ 9 }$$

where ${\varepsilon }_{\mathrm{elastic}}^{{nn}}$ and ${\varepsilon }_{\mathrm{active}}^{{nn}}$ 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. σⁿⁿ 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 ${{\boldsymbol{\varepsilon }}}_{\mathrm{elastic}}^{{nn}}$ is different for each porous region. For the electrodes, the Li insertion/extraction induced volume change ratios are different. ${\varepsilon }_{\mathrm{active}}^{{nn}}$ is inhomogeneous in both electrodes along the SLTD due to the nonuniform Li distribution during charge/discharge. For the separator, ${\varepsilon }_{\mathrm{active}}^{{nn}}$ 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 εⁿⁿ 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 I.

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 I, follow an Arrhenius-type temperature dependence in the form ¹⁰

$$\begin{eqnarray}&&(\bullet )={\left(\bullet \right)}_{\mathrm{ref}}\exp \left[\displaystyle \frac{{E}_{a}}{R}\left(\displaystyle \frac{1}{{\theta }_{0}}-\displaystyle \frac{1}{\theta }\right)\right]\end{eqnarray} \tag{ 18 }$$

where (•) represents D_s , D_e , κ_s , and k₀, and (•)_ref represents their corresponding reference values. Some of the electrolyte transport properties have a different dependence on temperature (see Table IV). Second, the local volumetric heat generation rate is computed by the DualFoil model based on Ref. 6

$$\begin{eqnarray}\begin{array}{rcl}\dot{q} & = & -\displaystyle \frac{1}{{L}_{0}}\left({\displaystyle \int }_{L}\displaystyle \sum _{k}{a}_{k}{i}_{n,k}{U}_{H,k}\,{dl}-{IV}\right),\quad \mathrm{with}\\ {U}_{H,k} & = & -T\displaystyle \frac{\partial {U}_{k}}{\partial T}+{U}_{k}\end{array}\end{eqnarray} \tag{ 19 }$$

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₀ of the sub-cell to compute the volume-averaged heat generation rate $\dot{q}$.

Following Ref. 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. ⁴⁵ 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 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.

Averaged quantities for the mechanical field

The mechanically coupled DualFoil model takes the averaged stress σⁿⁿ of the jellyroll in the SLTD to update the porosity in Eq. 9. σⁿⁿ represents the stress component of the homogenized Cauchy stress tensor σ computed from the thermomechanical model, 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 ₁, e ₂, e ₃) to denote the basis vectors of the Cartesian coordinate system. For an arbitrary direction n , to obtain the value of σⁿⁿ , we need to first rotate σ along e ₃ with an angle of δ to an intermediate coordinate system $({{\boldsymbol{e}}}_{1}^{{\prime} }$, ${{\boldsymbol{e}}}_{2}^{{\prime} }$, ${{\boldsymbol{e}}}_{3}^{{\prime} })$ with a rotation tensor Q ₃ defined as

$$\begin{eqnarray}{{\boldsymbol{Q}}}_{3}=\left[\begin{array}{ccc}\cos (\delta ) & -\sin (\delta ) & 0\\ \sin (\delta ) & \cos (\delta ) & 0\\ 0 & 0 & 1\end{array}\right].\end{eqnarray} \tag{ 20 }$$

The corresponding stress tensor ${\boldsymbol{\sigma }}^{\prime}$ expressed in the coordinate system $({{\boldsymbol{e}}}_{1}^{{\prime} }$, ${{\boldsymbol{e}}}_{2}^{{\prime} }$, ${{\boldsymbol{e}}}_{3}^{{\prime} })$ is computed as ${\boldsymbol{\sigma }}^{\prime} ={{\boldsymbol{Q}}}_{3}^{T}{\boldsymbol{\sigma }}{{\boldsymbol{Q}}}_{3}$. We then rotate ${\boldsymbol{\sigma }}^{\prime}$ along ${{\boldsymbol{e}}}_{2}^{{\prime} }$ with an angle of ϕ to the local coordinate system ( n , s , t ) with a rotation tensor Q ₂ defined as

$$\begin{eqnarray}{{\boldsymbol{Q}}}_{2}=\left[\begin{array}{ccc}\cos (\phi ) & 0 & \sin (\phi )\\ 0 & 1 & 0\\ -\sin (\phi ) & 0 & \cos (\phi )\end{array}\right].\end{eqnarray} \tag{ 21 }$$

The corresponding stress tensor $\tilde{{\boldsymbol{\sigma }}}$ expressed in the coordinate system ( n , s , t ) is computed as $\tilde{{\boldsymbol{\sigma }}}={{\boldsymbol{Q}}}_{2}^{T}{\boldsymbol{\sigma }}^{\prime} {{\boldsymbol{Q}}}_{2}={{\boldsymbol{Q}}}_{2}^{T}{{\boldsymbol{Q}}}_{3}^{T}{\boldsymbol{\sigma }}{{\boldsymbol{Q}}}_{3}{{\boldsymbol{Q}}}_{2}$. 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 σⁿⁿ for a specific electrochemical element is computed as

$$\begin{eqnarray}&&{\sigma }^{{nn}}=\displaystyle \frac{1}{{\displaystyle \sum }_{l=1}^{k}{\int }_{{{ \mathcal B }}_{l}}{dV}}\displaystyle \sum _{l=1}^{k}{\int }_{{{ \mathcal B }}_{l}}{\tilde{\sigma }}_{11}{dV},\end{eqnarray} \tag{ 22 }$$

where the integration is performed via the Gaussian quadrature rule with two Gauss points ^c being used in each direction. The same σⁿⁿ is assigned to each node of the 1D electrochemical mesh to compute the ${{\boldsymbol{\varepsilon }}}_{\mathrm{elastic}}^{{nn}}$ defined in 9. The thermomechanical model requires $\bar{\omega }$ to compute the F _swelling in 1, which is computed from the mechanically coupled DualFoil model as

$$\begin{eqnarray}&&\bar{\omega }=\displaystyle \frac{{\int }_{L}\overline{{\rm{\Omega }}(C)}{\epsilon }_{\mathrm{active}}^{0}{dl}}{{L}_{0}},\end{eqnarray} \tag{ 23 }$$

where $\overline{{\rm{\Omega }}(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. ^d One can refer to Ref. 12 for a more detailed description. The same $\bar{\omega }$ 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.

Zoom In Zoom Out Reset image size

Averaged quantities for the thermal field

The DualFoil model takes a volume-averaged temperature $\bar{\theta }$ 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 $\bar{\theta }$ for a specific sub-cell is computed via

$$\begin{eqnarray}&&\bar{\theta }=\displaystyle \frac{1}{{\displaystyle \sum }_{l=1}^{k}{\int }_{{{ \mathcal B }}_{l}}{dV}}\displaystyle \sum _{l=1}^{k}{\int }_{{{ \mathcal B }}_{l}}\theta {dV}.\end{eqnarray} \tag{ 24 }$$

where the integration is performed via the Gaussian quadrature rule with two Gauss points being used in each direction. The same $\bar{\theta }$ 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 rate $\dot{q}$ is calculated based on 19. The same $\dot{q}$ 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.

The electrochemical model is coupled with the thermomechanical model via volume-averaged quantities, such as temperature $\bar{\theta }$, heat generation rate $\dot{q}$, intercalation-induced volume change ratio $\bar{\omega }$, and mechanical stress σⁿⁿ . This information is exchanged between the 3D finite element code based on deal.II and the 1D DualFoil model code through 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 σⁿⁿ and volume-averaged temperature $\bar{\theta }$ from the 3D thermomechanical model as inputs and passes the heat generation rate $\dot{q}$ and intercalation-induced volume change ratio $\bar{\omega }$ 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 II.

Zoom In Zoom Out Reset image size

Table II. Summary of the simulation steps in the proposed multiphysics battery simulation framework.

A. Initialize a simulation with the proper initial state variables, boundary conditions, and physical parameters.

B. Compute quantities of interest at a new time step t

(a) Update the thermal loading and the mechanical loading. The former consists of an external convective boundary condition and an internal heat generation rate $\dot{q}$. The latter consists of an externally applied loading and an internal intercalation-induced swelling $\bar{\omega }$.

(b) Solve the temperature and displacement fields with the finite element method.

(c) Pass the volume-averaged temperature $\bar{\theta }$ and the projected stress value σnn to the DualFoil model via the API.

(d) Apply the electrical loading through the API. This is done indirectly by applying a guessed, uniform value of the voltage to each electrochemical sub-cell (e.g., converged voltage of the previous time step).

i. Independently solve the electrochemical properties and states of each sub-cell with finite difference method based on the input potential V, temperature $\bar{\theta }$, and stress σnn .

ii. Compute the sum of the output volumetric current densities of all the electrochemical sub-cells and check if it matches the applied volumetric current density.

iii. If the relative error err > tol with tol = 10−8, as shown in Fig. 3, adjust the applied sub-cell voltage V and repeat steps B(d)i and B(d)ii.

(e) Once converged, pass the volumetric quantities $\dot{q}$ and $\bar{\omega }$ from each sub-cell to its associated thermomechanical element(s) via the API.

(f) Output state variables of interest.

C. Move forward to the next time(loading) step.

D. Stop the simulation if certain criteria are met (e.g., a specified cutoff potential is reached).

First, we consider a brick-shape jellyroll with a size of 2.5 mm × 0.167 mm × 1.0 mm. The SLTD is parallel to the Y-axis. The jellyroll is meshed with 250 elements for the thermomechanical model, as shown in Fig. 4. Each thermomechanical element corresponds 1:1 to an electrochemical sub-cell, i.e., the fine mesh and the coarse mesh illustrated in Fig. 1 are identical with the ratio between the number of electrochemical elements and the number of thermomechanical elements being one. Each sub-cell is discretized into 80 intervals in the SLTD to solve for the electrochemical state variables. The material parameters listed in Table III–V are used. The jellyroll is gavanostatically discharged at I = 2 A with a fixed homogeneous temperature at θ = 20 °C under three different scenarios considering material parameter variations and mechanical loading: (i) uniform material distribution without external mechanical loading, (ii) non-uniform material distribution without external mechanical loading, where the initial anode porosity of electrochemical elements has a random distribution in the range of [0.226, 0.25], equivalent to a ±5% initial anode porosity variation, and (iii) uniform material distribution with non-uniform external mechanical loading in the SLTD that varies in a linear way along the X-axis with σⁿⁿ = −5 MPa at X = 0 and σⁿⁿ = 5 MPa at X = 2.5 mm. In all three cases, the porosity may change due to swelling and shrinkage of the active material.

Zoom In Zoom Out Reset image size

The current density and SOC ^e 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 Figs. 5a, 5b. 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 Figs. 5c–5f. These results confirm the capability of the newly proposed framework to capture spatially varying quantities in a Li-ion battery.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Next, an 18 650 cell is discharged at 4A rate under both natural (h = 12.5 W/(m² · K)) and forced (h = 49.5 W/(m² · K)) convection conditions, with h as the external heat transfer coefficient. The original experimental setup can be found in Ref. 14. In the simulation, we use a fine mesh with 61,700 elements for the thermomechanical model, as shown in Fig. 6, 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).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

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 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 18 650 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. 7b, 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 non-uniform 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.

Zoom In Zoom Out Reset image size

Lastly, a prismatic cell, which has been previously studied and parameterized in Refs. 12,59–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/8 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 Ref. 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₂ active particles as inputs and outputs the averaged thickness change of the whole jellyroll. As pointed out in Ref. 12, the observed difference in the thickness ratio change arises because the swelling ratio vs SOC for both graphite and LiCoO₂ used in the simulation is measured from a different cell.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

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. 12b, 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.

Zoom In Zoom Out Reset image size