## Abstract

The behavior of lithium-ion batteries (LIBs) under mechanical loading is a complex multiphysics process including mechanical deformation, internal short circuit, and thermal runaway. To deeply understand the mechanism of battery failure and accurately predict the onset of internal short circuit and thermal runaway, a multiphysics-based computation framework of LIBs is in pressing need. In this article, a multiphysics model that couples five submodels (mechanical model, internal short-circuit model, battery model, heat transfer model, and thermal runaway model) is established to predict the evolution of force, voltage, and temperature under steel ball compression. The suitable agreement between simulation results and experimental data of batteries with different state of charges demonstrates that the proposed model is capable of predicting the multiphysical behavior of the battery. Further, a systematic parametric study is conducted to investigate the short-circuit triggering and temperature rise of batteries under different conditions, and the workflow of battery safety optimal design is proposed by applying the multiphysics model.

## 1 Introduction

Due to the high energy density and long lifetime, lithium-ion batteries (LIBs) have been regarded as one of the most promising alternative energy sources [1], with the extensive application from consumer electronic devices to electric-drive vehicles [2,3]. Internal short circuit and thermal runaway can result in catastrophic consequences such as smoke, fire, and explosion [4]. Generally speaking, the failure of LIBs can be caused by both external abuse conditions (impact [5,6], overcharge [7], overheating, [8]) and internal abuse conditions (internal short circuit [9,10] and aging [11,12]), of which the internal short circuit caused by external mechanical abuse [3,13,14] is a very practical problem.

Mechanical loading such as bending [15], compression [16], and indentation [17] can lead to the deformation of LIBs, which causes the failure of components and further the internal short circuit of LIBs. The Joule heat generated from the internal short circuit can result in a drastic temperature increase. Then, thermal runaway is triggered, and the related chemical reactions [18,19] are gradually activated once the temperature reaches the threshold. Finally, the extremely high pressure inside the battery shell due to chemical reactions can cause a possible explosion. Above all, the behavior of LIBs under mechanical abuse is the typical multiphysics progress involving mechanical deformation, electrochemical reaction, chemical reaction, and thermal reaction.

In recent years, multiphysics models for mechanical abuse conditions have been gradually developed to unveil the failure mechanism of the battery and predict the internal short circuit and thermal runaway behavior at the cell level [20–24], component level [25,26], and particle level [27,28]. Jia et al. [20] investigated the relationship between electrical behavior and force response of LIBs under impact loading based on a validated mechanical-electrochemical model. To obtain a comprehensive and accurate multiphysics model, a coupled mechanical-electrical-thermal model for cylindrical LIBs, which contains five submodels, is established by Yuan et al. [21] and validated based on indentation, compression, and three-point bending experiments. Liu et al. [22] combined the multiphysics models and physical tests to investigate the underlying mechanism of the internal short circuit. They proved that mechanical deformation, stress distribution, and failure processes act together to determine the particular mode of the internal short circuit (ISC). Moreover, Jia et al. [24] explored the safety characteristics of defective LIBs and proved that thermal runaway was more prone to be triggered in defective LIBs. Gao et al. [25] developed a multiphysics-multiscale model of Si–C composite anode to study the effect of Si percentage, mechanical constraint, and charging rate. Further, they established an electro-chemo-mechanical model to investigate the behavior of five Si/C composite nanostructures [28].

However, current researches mainly focus on the accuracy of the multiphysics model and further the mechanism analysis based on the model. For better guidance of battery safety design, the safety boundary of battery need to be obtained based on the comprehensive parametric study. In this article, a multiphysics model is established based on the experimental results to describe the typical behavior of pouch LIB cell under steel ball compression, including mechanical failure, internal short circuit, and thermal runaway. Based on the verified multiphysics model, the parametric study is conducted to figure out the safety boundary of battery, which is vital for the guidance of battery safety design.

## 2 Methods

### 2.1 Coupling Strategies.

The Li-ion battery with the size of 97 mm × 58 mm × 4.5 mm supplied by CATL is chosen as the target sample in this study. The nominal capacity of this battery is 3.577 Ah, and its charge/discharge cutoff voltages are 4.2 and 2.5 V, respectively. A multiphysics model of the LIB is established in this section, which consists of five submodels, i.e., mechanical model, internal short-circuit model, battery model, heat transfer model, and thermal runaway model (Fig. 1). The five submodels are coupled with one another through key parameters. With the input boundary conditions, the mechanical model can compute the mechanical deformation and transmit the strain field to the internal short-circuit model. The internal short-circuit model calculates the internal short-circuit resistance according to the strain field and transmits it to the battery model. In exchange, the battery model returns the calculated voltage and internal short-circuit current back to the internal short-circuit model. The heat transfer model, internal short-circuit model, battery model, and thermal runaway model are coupled through heat and temperature.

To simplify the computation of the multiphysics model, different geometries are adopted in submodels considering that only relevant variables are exchanged between submodels. According to the symmetry of loading condition, the mechanical model is simplified to a 2D axisymmetric solid mechanical model, which only describes the stress field around the loading region. The 1D lithium-ion battery model is adopted since the battery model only needs to provide voltage and current for the internal short-circuit model and thermal model. Due to the need to compute the temperature distribution in the whole region of the battery, a homogenized 3D solid heat transfer model was adopted to reach a rapid computational convergence and lower the computation cost.

Applying the multiphysics model, we can simulate the whole process of battery under mechanical abuse conditions from mechanical failure to internal short circuit and thermal runaway. The force, voltage, and temperature distributions are predicted for different input conditions such as loading, state of charge (SOC), and initial temperature.

### 2.2 Battery Model.

*ɛ*

^{l}is the volume fraction,

*c*

^{l}represents the Li

^{+}concentration, $Deffl$ indicates the modified diffusion coefficient,

*i*

^{l}is the current density, $\kappa effl$ represents the modified electrical conductivity,

*φ*

_{l}is the potential of the electrolyte,

*t*

_{+}is the Li

^{+}transfer number, $1+dlnfdlncl$ is the molar activity coefficient,

*R*

_{g}is the gas constant, and

*F*is Faraday’s constant.

*q*

_{r}and

*q*

_{j}are produced reaction heat and joule heat, respectively,

*L*is the length of the 1D battery model, and

*ϕ*

^{s}represents the potential of the electrode. The key parameters of the battery model are listed in Table 1.

Symbol | Parameter | Unit | Anode | Cathode | Separator |
---|---|---|---|---|---|

κ_{s} | Solid-phase conductivity | S/m | 270 | 100 | N/A |

D | Solid-phase diffusivity | m^{2}/s | 4e-14 | 1.8e-14 | N/A |

R | Particle radius | m | 1.6e-5 | 9e-6 | N/A |

C | Capacity | Ah | 3.577 |

Symbol | Parameter | Unit | Anode | Cathode | Separator |
---|---|---|---|---|---|

κ_{s} | Solid-phase conductivity | S/m | 270 | 100 | N/A |

D | Solid-phase diffusivity | m^{2}/s | 4e-14 | 1.8e-14 | N/A |

R | Particle radius | m | 1.6e-5 | 9e-6 | N/A |

C | Capacity | Ah | 3.577 |

Charge and discharge experiments with rates of 0.1C, 0.2C, 0.5C, and 1C are performed on the target battery cells. The simulated voltage corresponds well with the experimental results for all C rates (Fig. 2).

### 2.3 Mechanical Model.

A detailed mechanical model, which includes a separator, positive and negative active material, aluminum, and copper collector, was built in the present study. To improve the computational efficiency, the 2D axisymmetric geometry was adopted for the mechanical modeling (Fig. 3(a)). The *z*-direction of the battery lower surface and the *r*-direction of the battery central axis are fixed, while the steel ball is given an *z*-direction velocity of 1 mm/min. A penalty-based contact is set for the contact parts assuming a friction coefficient of 0.3.

*σ*

_{p}and plastic strain

*ɛ*

_{p}satisfy the following equation:

*a*,

*b*, and

*n*are fitting parameters, which are validated by comparing the simulated and experiment results in terms of tension and compression tests of the components (Figs. 4(a)–4(e)). Furthermore, the comparison of the experiment and simulation results of steel ball compression indicates that the mechanical model is capable to accurately predict the mechanical behavior of the battery.

### 2.4 Internal Short-Circuit Model.

*A*is defined as the area where the strain of the separator is greater than the initial internal short-circuit strain $\epsilon initaleq$, which can be observed in Fig. 3(b) and expressed as follows:

*t*is the separator thickness and

*V*is the volume of which the strain is larger than $\epsilon initaleq$. The area resistance of internal short-circuit contact $RArSt$ is obtained from experiments, and there are different area resistances for different contact types. In this study, the area resistance range can be expressed as follows considering the complex contact situation of the real internal short-circuit process often mixed with several contact types:

*σ*is the failure stress of the internal short circuit.

*α*is 0.1. $RArSt$ can be obtained that

*f*(SOC) ∈ (0, 1] is defined as the correction coefficient. The correction coefficients at different SOCs and corresponding internal short-circuit resistances are obtained (Table 2) based on the experimental results.

### 2.5 Thermal Transfer Model.

*ρ*is the density of the material,

*C*

_{p}is the heat capacity of the LIB,

*T*is the temperature,

*Q*is the total heat source, and

**is the heat flux written as follows:**

*q**Q*

_{cell}is the heat generated during the battery discharge process, provided by the battery model. The thermal convection power with the air is expressed as follows:

*A*

_{c}is the contact area between the battery and air,

*T*

_{amb}is the ambient temperature, and

*h*

_{c}is the convection coefficient.

*A*

_{t}includes the lower surface area of the battery and the upper surface area in contact with the indenter and

*h*

_{t}is the equivalent heat transfer coefficient. The heat sources and thermal boundary conditions are presented in Fig. 3(c).

### 2.6 Thermal Runaway Model.

The thermal runaway model describes the decomposition reaction rates and calculates the heat generated from reactions through Arrhenius-type equations. The main reactions and equations are as follows:

- Solid electrolyte interphase (SEI) decomposition:The decomposition rate is:(22)$dcseidt=dcsei,ddt+dcsei,gdt$The regeneration rate is:(23)$dcsei,ddt=\u2212Aseiexp[\u2212Ea,seiRgT]csei,d$where(24)$dcsei,gdt=Ksei,gcsei,g$
*c*_{sei},*c*_{sei,d}, and*c*_{sei,g}represent the dimensionless amount of lithium-containing metastable species in SEI, decomposited SEI, and regenerated SEI, respectively.*A*_{sei}and*E*_{a,sei}are the Arrhenius constant and experimental activation energy of SEI decomposition, respectively.*R*_{g}is the gas constant. - The cathode–electrolyte reaction:where(25)$d\alpha c\u2212edt=Ac\u2212e\alpha c\u2212e(1\u2212\alpha c\u2212e)exp[\u2212Eac\u2212eRgT]$
*α*^{c−e},*A*^{c−e}, and $Eac\u2212e$ are the conversion degree, Arrhenius constant, and experimental activation energy of cathode–electrolyte reaction, respectively. - Electrolyte decomposition:where(26)$dcedt=\u2212Aeexp[\u2212EaeRgT]ce$
*c*_{e}represents the dimensionless amount of electrolytes, and*A*^{e}and $Eae$ are the Arrhenius constant and experimental activation energy of electrolyte decomposition, respectively. - Anode decomposition:where(27)$dcadt=\u2212Aaexp[\u2212cseicsei0]exp[\u2212EaaRgT]ca$
*c*_{a}represents the dimensionless amount of lithium intercalated within carbon.*c*_{sei0}represents the initial value of the dimensionless amount of lithium-containing metastable species in the SEI.*A*^{a}and $Eaa$ are the Arrhenius constant and experimental activation energy of anode decomposition, respectively. - Separator decomposition:where(28)$dcsepdt=\u2212Asepexp[\u2212EasepRgT]csep$
*c*_{sep}represents the dimensionless amount of seperator.*A*^{sep}and $Easep$ are the Arrhenius constant and experimental activation energy of separator decomposition, respectively. - The internal short circuit caused by direct contact between cathode and anode due to separator decomposition, which is related to SOC:where(29)$dSOCdt=\u2212Aecexp[\u2212EaecRgT]SOC\u2212iC$
*A*^{ec}and $Eaec$ are the Arrhenius constant and experimental activation energy of internal short circuit, respectively.The total heat production rate in the thermal runaway model is the summation of all reactions:where(30)$dQxdt=Hxmxdcxdt(Tx>Tcx)dQdt=\u2211dQxdt$*H*^{x}is the reaction enthalpy change,*m*_{x}is the mass density of the reactants, and $Tcx$ is the trigger temperature of the reaction. The enthalpy change of reaction (6) can be expressed as follows:where(31)$Hec=3600Ucell2Ccell\eta $*U*_{cell}is the cell voltage,*C*_{cell}is the capacity, and*η*is the heat efficiency (the value of*η*is adjusted according to the temperature peak).

## 3 Results

The force, voltage, and temperature behavior of batteries with different SOCs under steel ball compression are recorded and analyzed (Figs. 5 and 6). The results indicate that an internal short circuit occurs in all SOC conditions, while thermal runaway only occurs in cases of high SOCs (0.8 and 1).

### 3.1 Internal Short Circuit Without Thermal Runaway.

For cases of low SOCs (0.1, 0.3, and 0.6), the battery exhibits typical internal short-circuit behavior with no thermal runaway (Fig. 5). The force increases continuously at first and then drops drastically accompanied by voltage drop and temperature increase caused by the internal short circuit. The comparison between the simulation and experiment results under different SOC conditions (0.1, 0.3, 0.6) indicates that the multiphysics model can well predict the force peak, voltage drop, and temperature rise in cases of internal short circuit without thermal runaway.

### 3.2 Thermal Runaway Cases.

For cases of high SOCs (0.8, 1), both internal short circuit and thermal runaway happen in the steel ball compression tests (Fig. 6). The multiphysics model can accurately predict the evolution process of force, voltage, and temperature, especially the temperature peak and cooling history of the battery.

## 4 Discussion

Applying the validated multiphysics model, a series of parametric studies are conducted to investigate short-circuit triggering and temperature rise of batteries under different conditions. It mainly includes three parts:

The effect of steel ball size, steel ball mass, loading speed, and the number of battery layers on internal short-circuit triggering under quasi-static and impact loading conditions.

The effect of internal short-circuit position and battery size on temperature rise of the whole battery after short-circuit triggering.

Critical internal short-circuit resistance of thermal runaway triggering.

Based on the aforementioned three parts, a total of about 6500 groups of parametric calculations were performed.

### 4.1 Internal Short-Circuit Triggering Under Quasi-Static Loading.

For cases under quasi-static loading, the influences of steel ball size (*a*, *b*) and battery thickness (number of layers *N*_{layer}) on internal short-circuit triggering force *F* and displacement *S* are studied in this section (Fig. 7). The steel ball size *a* and *b* change with the range from 1 mm to 4 mm, while the number of battery layers *N*_{layer} is set from 6 to 18. It can be observed that the ISC displacement *S* and ISC force *F* tend to decrease with increasing *a*, decreasing *b*, and decreasing *N*_{layer}. The internal short circuit is triggered in advance due to the earlier failure of the separator under the circumstances of a sharper steel ball and thinner battery. Besides, it can be found that *S* is more sensitive to *b* and *N*_{layer} and *F* is more sensitive to *b*.

### 4.2 Internal Short-Circuit Triggering Under Impact Loading.

For cases under impact loading, five critical parameters (Fig. 8(a)) including steel ball size *a* (1–100 mm) and *b* (1–100 mm), steel ball mass *m* (1–10 kg), loading speed *v* (1–10 m/s), and the number of battery layers *N*_{layer} (6–48) are investigated to figure out the relationship between these parameters and internal short-circuit triggering.

First, the coupling effect of *a*, *b*, *m*, and *v* is presented in Figs. 8(b)–8(e), where the blue and red marks represent cases with no internal short circuit (safe cases) and cases with internal short circuit (dangerous cases), respectively. The results show that the internal short circuit is more prone to happen while increasing *a*, *m*, *v,* and decreasing *b* since the separators of the battery are easier to fail under the impact of a sharper steel ball with larger momentum. In addition, among all these four parameters, *b* is the primary factor that influences the triggering of the internal short circuit. The boundary of the internal short circuit is clarified in Fig. 8(d).

The coupling effect of *a* and *N*_{layer} is further analyzed (Fig. 9). The larger occupying space of ISC cases in the parameter domain indicate higher risk of internal short circuit as *N*_{layer} increases. Internal short circuit is more prone to be triggered in thinner batteries under impact loading. The internal short circuit will be triggered when the strain of the separator reaches the internal short circuit strain. Under the same deformation, the thinner battery will have a larger strain.

### 4.3 Effect of Internal Short-Circuit Positions.

The maximum (time) average (battery) temperature rise of the battery after an internal short circuit is greatly affected by the internal short-circuit position. Different internal short-circuit positions *x* (0.25–0.75), *y* (0.25–0.75), and *z* (0–0.5) have been set to investigate their influence on the maximum (time) average (battery) temperature increase under several SOCs (0.1–1). First, the coupling effect of *x* and *y* is studied with *z* = 0 (Fig. 10(b)). For batteries of low SOCs (0.1–0.6), compared to the internal short circuit occurring at the battery center, the internal short circuit occurring at four battery corners leads to higher temperature rise. However, for high SOC batteries (0.8–1), a higher temperature rise occurs when the internal short circuit occurs at the center of the battery. For low SOC cells, the TR will not be triggered. The heat generation rate is relatively low and only affects the ISC area. The heat transfer in the cell domain will influence the temperature rise more. The internal short-circuit heat applied at the center of the battery will be distributed quickly. For high SOC cells, the TR will be triggered. The TR heat power is very large and applied to the whole domain. Thus, the influence of heat transfer within the cell domain is very small. The internal short circuit at the center of the cell will lead to a slightly large temperature rise due to the boundary conditions (heat transfer at the boundary). The influence of longitudinal coordinate *z* is presented in Figs. 10(c)–10(e). It is obvious that when the short-circuit position is in the center of the battery, changing the coordinate *z* of internal short-circuit position hardly influences the temperature rise of the battery. Once the internal short circuit occurs in the corner of the battery, the maximum temperature rise is sensitive to *z*. The battery temperature rise reaches the highest when *z* = 0.

### 4.4 Effect of Battery Dimensions.

Battery dimensions significantly influence the maximum (time) average (battery) temperature rise after internal short circuit triggers. Here, *w* (0.5–2), *l* (0.5–2), and *t* (0.5–2) are adopted to represent the scaling coefficient of the three dimensions of battery (Fig. 11(a)). Keeping the internal short-circuit position at battery center, the coupling effect of *w*, *l*, and *t* for different SOCs is depicted in Figs. 11(b)–11(e). For cases with low SOCs (0.1–0.6), the temperature rise first decreases and then increases when increasing *w* and *l*. For cases with high SOCs (0.8–1), the changing trend of temperature rise behaves differently according to different *t*.

The coupling effect of internal short circuit positions and battery dimensions is further investigated. Figures 12(b)–12(e) presents the maximum (time) average (battery) temperature rise for four short-circuit positions P1 (surface center), P2 (body center), P3 (surface corner), and P4 (surface edge), respectively (Fig. 12(a)). It can be found that changing internal short-circuit position on the upper surface of battery (P1, P3, and P4) will not influence the relationship between temperature rise and battery dimensions obviously. However, when the internal short-circuit position changes from the upper surface (P1, P3, and P4) to body center (P2), the overall temperature rise decreases for cases of 0.6 SOC and 0.8 SOC due to the changing of heat dissipation condition.

### 4.5 Thermal Runaway Triggering Boundary.

*R*

_{c}, thermal runaway will not be triggered.

### 4.6 Optimal Design Methodologies.

The battery safety optimal design can be realized by applying the multiphysics model (Fig. 14(a)). With the input parameters including target condition and battery parameters, the output results can be obtained from the multiphysics model such as internal short-circuit triggering force/displacement and temperature distribution. Based on the output results, the battery parameters can be optimized according to the optimized design including minimalizing temperature rise and increasing tolerable load. Then, the optimized parameters become the input parameters of the multiphysics model to achieve a better safety performance. Taking the battery dimensions optimization as an example (Fig. 14(b)), switching the battery dimensions from line 1 (current design) to line 2 (decreasing *w* and increasing *l*) can decrease the temperature rise at high SOCs while keeping the good performance at low SOCs.

## 5 Conclusions

In summary, a multiphysics model that contains mechanical model, internal short-circuit model, battery model, heat transfer model, and thermal runaway model is established and validated based on the experimental results. The model can well predict the evolutions of force, voltage, and temperature during the whole process of mechanical loading from mechanical deformation to internal short circuit and thermal runaway. Parametric studies in terms of governing factors, such as steel ball size, loading speed, and number of battery layers are conducted and discussed. According to the results, several conclusions can be summarized as follows:

The ISC displacement

*S*and ISC force*F*tend to decrease with increasing*a*, decreasing*b*, and decreasing*N*_{layer}under quasi-static loading. For cases under impact loading, the internal short circuit is more prone to happen while increasing*a*,*m*,*v,*and decreasing*b*.For batteries of low SOCs (0.1–0.6), the internal short circuit occurring at four battery corners leads to a higher temperature rise. However, for high SOC batteries (0.8–1), a higher temperature rise occurs when the internal short circuit occurs at the center of the battery.

For cases with low SOCs (0.1–0.6), the temperature rise first decreases and then increases when increasing

*w*and*l*. For cases with high SOCs (0.8–1), the changing trend of temperature rise behaves differently at different*t*.

This model can lay a solid foundation to understand the mechanical integrity of LIBs and guide the safety design of the next-generation battery.

## Conflict of Interest

There are no conflicts of interest.

## Data Availability Statement

The datasets generated and supporting the findings of this article are obtainable from the corresponding author upon reasonable request. The authors attest that all data for this study are included in the paper.