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 [2024], 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.

Fig. 1
Multiphysics modeling framework
Fig. 1
Multiphysics modeling framework
Close modal

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.

The battery model was developed based on the classical 1D battery model, which is first established by Newman and coworkers [29]. This model is used to predict the voltage drop and heat produced under discharging conditions by solving the following governing equations [30]:
εlclt=Defflcl+ilt+F
(1)
il=κefflφl+2κefflRgTF(1+dlnfdlncl)(1t+)lncl
(2)
where ɛl is the volume fraction, cl represents the Li+ concentration, Deffl indicates the modified diffusion coefficient, il is the current density, κ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, Rg is the gas constant, and F is Faraday’s constant.
The total heat generation of the battery system includes produced reaction, joule, the resistance of current collector, polarization, and irreversible heat [31]. In this study, only the produced reaction and joule heat are considered:
qtotal=qr+qj
(3)
qr=0Li(ϕsϕlEeq)dx
(4)
qj=0Lκeffs(δϕsδx)2+κeffl(δϕlδx)2+κeffD(δlnclδx)(δϕlδx)dx
(5)
where qr and qj 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.
Table 1

Key parameters of the battery model

SymbolParameterUnitAnodeCathodeSeparator
κsSolid-phase conductivityS/m270100N/A
DSolid-phase diffusivitym2/s4e-141.8e-14N/A
RParticle radiusm1.6e-59e-6N/A
CCapacityAh3.577
SymbolParameterUnitAnodeCathodeSeparator
κsSolid-phase conductivityS/m270100N/A
DSolid-phase diffusivitym2/s4e-141.8e-14N/A
RParticle radiusm1.6e-59e-6N/A
CCapacityAh3.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).

Fig. 2
Validation of a battery model. Comparison of cycling curves between experimental results and simulation results at different C-rate: (a) 0.1C, (b) 0.2C, (c) 0.5C, and (d) 1C.
Fig. 2
Validation of a battery model. Comparison of cycling curves between experimental results and simulation results at different C-rate: (a) 0.1C, (b) 0.2C, (c) 0.5C, and (d) 1C.
Close modal

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.

Fig. 3
Setups of the submodels: (a) setup of the mechanical model, (b) boundary conditions of the thermal transfer model, and (c) schematic diagram of the internal short-circuit model
Fig. 3
Setups of the submodels: (a) setup of the mechanical model, (b) boundary conditions of the thermal transfer model, and (c) schematic diagram of the internal short-circuit model
Close modal
The main equation given by Newton’s second law is expressed as follows:
ρ2ut=(FS)T+FV
(6)
where u is the displacement field, S is the Piola–Kirchhoff stress tensor, FV is the body force, and F is the deformation gradient, which are denoted as follows:
F=I+u
(7)
The linear elastic model and Johnson–Cook plastic model are selected to characterize the mechanical properties of the components. The plastic stress σp and plastic strain ɛp satisfy the following equation:
σp=a+bεpn
(8)
where 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.
Fig. 4
Validation of the mechanical model. Comparison of force–displacement curve between experimental results and simulation results: tension test results of (a) anode, (b) cathode, and (c) separator; compression test results of (d) anode and (e) cathode; (f) cell test results of 8 mm sphere indentation.
Fig. 4
Validation of the mechanical model. Comparison of force–displacement curve between experimental results and simulation results: tension test results of (a) anode, (b) cathode, and (c) separator; compression test results of (d) anode and (e) cathode; (f) cell test results of 8 mm sphere indentation.
Close modal

2.4 Internal Short-Circuit Model.

The internal short-circuit model is built to calibrate the internal short-circuit resistance, internal short-circuit area, and heat produced. The internal short-circuit area is computed based on the strain field obtained from the mechanical model and equivalent to the regional heat source in the thermal model. Specifically, the internal short-circuit area A is defined as the area where the strain of the separator is greater than the initial internal short-circuit strain εinitaleq, which can be observed in Fig. 3(b) and expressed as follows:
A=1tcdV,c={0εeq<εinitaleq1εeqεinitaleq
(9)
where t is the separator thickness and V is the volume of which the strain is larger than ε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:
RArSt(RAAnCaSt(σ),RAAlCuSt(σ))=(1.52×107,1.27×105)(Ωm2)
(10)
where σ is the failure stress of the internal short circuit.
It is assumed that the contact area of the positive and negative active material is one order of magnitude larger than the contact area of the positive and negative collector, which means:
ArSt=AAnCaSt+AAlCuSt
(11)
AAnCaSt/AAlCuSt=(1α)/α
(12)
where the value of α is 0.1. RArSt can be obtained that
RArSt=RAAlCuSt(σ)RAAnCaSt(σ)(1α)RAAlCuSt(σ)+αRAAnCaSt(σ)
(13)
The internal short-circuit resistance can be expressed as follows:
Rst=RArSt/ArSt
(14)
Considering the effect of SOC, the internal short-circuit resistance can be modified as follows:
Rst=Rst0f(SOC)
(15)
where Rst0 is the internal short-circuit resistance at 0.6 SOC, 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.
Table 2

Correction coefficients of the internal short-circuit resistances at different SOCs

SOC0.10.30.60.81
Rst(soc) (Ω)0.611.20.230.2
f(soc)0.50.8310.190.17
SOC0.10.30.60.81
Rst(soc) (Ω)0.611.20.230.2
f(soc)0.50.8310.190.17

2.5 Thermal Transfer Model.

The governing equation of the thermal transfer model is as follows:
ρCpTt+q=Q
(16)
where ρ is the density of the material, Cp is the heat capacity of the LIB, T is the temperature, Q is the total heat source, and q is the heat flux written as follows:
q=kT
(17)
The total heat source is equal to the sum of all heat sources:
Q=Qst+Qcell+Qconv+Qtrans
(18)
where the heat source from the internal short-circuit model is expressed as follows:
Qst=I2Rst
(19)
Qcell 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:
Qconv=Achc(TTamb)
(20)
where Ac is the contact area between the battery and air, Tamb is the ambient temperature, and hc is the convection coefficient.
The equivalent thermal transfer power between supporter and battery is expressed as follows:
Qtrans=Atht(TTamb)
(21)
where At includes the lower surface area of the battery and the upper surface area in contact with the indenter and ht 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:

  1. Solid electrolyte interphase (SEI) decomposition:
    dcseidt=dcsei,ddt+dcsei,gdt
    (22)
    The decomposition rate is:
    dcsei,ddt=Aseiexp[Ea,seiRgT]csei,d
    (23)
    The regeneration rate is:
    dcsei,gdt=Ksei,gcsei,g
    (24)
    where csei, csei,d, and csei,g represent the dimensionless amount of lithium-containing metastable species in SEI, decomposited SEI, and regenerated SEI, respectively. Asei and Ea,sei are the Arrhenius constant and experimental activation energy of SEI decomposition, respectively. Rg is the gas constant.
  2. The cathode–electrolyte reaction:
    dαcedt=Aceαce(1αce)exp[EaceRgT]
    (25)
    where αce, Ace, and Eace are the conversion degree, Arrhenius constant, and experimental activation energy of cathode–electrolyte reaction, respectively.
  3. Electrolyte decomposition:
    dcedt=Aeexp[EaeRgT]ce
    (26)
    where ce represents the dimensionless amount of electrolytes, and Ae and Eae are the Arrhenius constant and experimental activation energy of electrolyte decomposition, respectively.
  4. Anode decomposition:
    dcadt=Aaexp[cseicsei0]exp[EaaRgT]ca
    (27)
    where ca represents the dimensionless amount of lithium intercalated within carbon. csei0 represents the initial value of the dimensionless amount of lithium-containing metastable species in the SEI. Aa and Eaa are the Arrhenius constant and experimental activation energy of anode decomposition, respectively.
  5. Separator decomposition:
    dcsepdt=Asepexp[EasepRgT]csep
    (28)
    where csep represents the dimensionless amount of seperator. Asep and Easep are the Arrhenius constant and experimental activation energy of separator decomposition, respectively.
  6. The internal short circuit caused by direct contact between cathode and anode due to separator decomposition, which is related to SOC:
    dSOCdt=Aecexp[EaecRgT]SOCiC
    (29)
    where Aec 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:
    dQxdt=Hxmxdcxdt(Tx>Tcx)dQdt=dQxdt
    (30)
    where Hx is the reaction enthalpy change, mx 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:
    Hec=3600Ucell2Ccellη
    (31)
    where Ucell is the cell voltage, Ccell 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).

Fig. 5
Results of the nonthermal runaway cases: (a) cell tests experimental setup and force/voltage/temperature time curves at different SOCs: (b) 0.1, (c) 0.3, and (d) 0.6
Fig. 5
Results of the nonthermal runaway cases: (a) cell tests experimental setup and force/voltage/temperature time curves at different SOCs: (b) 0.1, (c) 0.3, and (d) 0.6
Close modal
Fig. 6
Typical results of thermal runaway cases. Force/voltage/temperature time curves at different SOCs: (a) 0. 8 and (b) 1.
Fig. 6
Typical results of thermal runaway cases. Force/voltage/temperature time curves at different SOCs: (a) 0. 8 and (b) 1.
Close modal

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:

  1. 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.

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

  3. 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 Nlayer) 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 Nlayer 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 Nlayer. 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 Nlayer and F is more sensitive to b.

Fig. 7
Parametric analysis of ISC triggering under quasi-static loading: (a) selected parameters, (b) ISC triggering displacement, and (c) ISC triggering force
Fig. 7
Parametric analysis of ISC triggering under quasi-static loading: (a) selected parameters, (b) ISC triggering displacement, and (c) ISC triggering force
Close modal

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 Nlayer (6–48) are investigated to figure out the relationship between these parameters and internal short-circuit triggering.

Fig. 8
Parametric analysis of ISC triggering under impact loading: (a) selected parameters, and ISC triggering behaviors at (b) a = 5, (c) a = 20, (d) a = 50, and (e) a = 100
Fig. 8
Parametric analysis of ISC triggering under impact loading: (a) selected parameters, and ISC triggering behaviors at (b) a = 5, (c) a = 20, (d) a = 50, and (e) a = 100
Close modal

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 Nlayer is further analyzed (Fig. 9). The larger occupying space of ISC cases in the parameter domain indicate higher risk of internal short circuit as Nlayer 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.

Fig. 9
Effects of indenter width and battery layer numbers on ISC triggering under impact loading. ISC triggering behaviors when a = 5 at (a) N = 6, (b) N = 12, (c) N = 24, and when a = 20 at (d) N = 6, (e) N = 12, (f) N = 24.
Fig. 9
Effects of indenter width and battery layer numbers on ISC triggering under impact loading. ISC triggering behaviors when a = 5 at (a) N = 6, (b) N = 12, (c) N = 24, and when a = 20 at (d) N = 6, (e) N = 12, (f) N = 24.
Close modal

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.

Fig. 10
Parametric analysis of temperature increase for different ISC positions. (a) Selected parameters, and (b) maximum (time) average (battery) temperature increase for different ISC positions with z = 0. Temperature increase behaviors with different z values: (c) z = 0, (d) z = 0.25, and (e) z = 0.5.
Fig. 10
Parametric analysis of temperature increase for different ISC positions. (a) Selected parameters, and (b) maximum (time) average (battery) temperature increase for different ISC positions with z = 0. Temperature increase behaviors with different z values: (c) z = 0, (d) z = 0.25, and (e) z = 0.5.
Close modal

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.

Fig. 11
Parametric analysis of temperature increase for different battery dimensions. (a) Selected parameters, and average temperature increase for different battery dimensions: (b) t = 0.5, (c) t = 1, (d) t = 1.25, and (e) t = 2.
Fig. 11
Parametric analysis of temperature increase for different battery dimensions. (a) Selected parameters, and average temperature increase for different battery dimensions: (b) t = 0.5, (c) t = 1, (d) t = 1.25, and (e) t = 2.
Close modal

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.

Fig. 12
Coupling effects of ISC positions and battery dimensions on temperature increase. (a) Selected ISC positions, and average temperature increase at different battery dimensions with different ISC positions: (b) P1, (c) P2, (d) P3, and (e) P4
Fig. 12
Coupling effects of ISC positions and battery dimensions on temperature increase. (a) Selected ISC positions, and average temperature increase at different battery dimensions with different ISC positions: (b) P1, (c) P2, (d) P3, and (e) P4
Close modal

4.5 Thermal Runaway Triggering Boundary.

The triggering boundary of thermal runaway is investigated and clarified in this section. The critical internal short circuit resistance of thermal runaway is highlighted as curves under different ambient temperatures (Fig. 13(a)), SOCs (Fig. 13(b)), and short-circuit contact types (Fig. 13(c)). Based on the results, the formula is derived to represent the critical resistance as follows:
lnRc=0.15e0.03T0+1.5lnSOC0.3ln(RA)4.7
(32)
When the short-circuit resistance is greater than the critical resistance Rc, thermal runaway will not be triggered.
Fig. 13
Thermal runaway triggering boundary of internal short-circuit resistance under different (a) ambient temperatures, (b) SOCs, and (c) short-circuit contact types
Fig. 13
Thermal runaway triggering boundary of internal short-circuit resistance under different (a) ambient temperatures, (b) SOCs, and (c) short-circuit contact types
Close modal

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.

Fig. 14
(a) Workflow and (b) example for the battery safety optima l design
Fig. 14
(a) Workflow and (b) example for the battery safety optima l design
Close modal

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:

  1. The ISC displacement S and ISC force F tend to decrease with increasing a, decreasing b, and decreasing Nlayer 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.

  2. 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.

  3. 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.

References

1.
Wang
,
C.-Y.
,
Zhang
,
G.
,
Ge
,
S.
,
Xu
,
T.
,
Ji
,
Y.
,
Yang
,
X.-G.
, and
Leng
,
Y.
,
2016
, “
Lithium-Ion Battery Structure That Self-Heats at Low Temperatures
,”
Nature
,
529
(
7587
), pp.
515
518
.
2.
Jia
,
Y.
,
Li
,
J.
,
Yuan
,
C.
,
Gao
,
X.
,
Yao
,
W.
,
Lee
,
M.
, and
Xu
,
J.
,
2021
, “
Data-Driven Safety Risk Prediction of Lithium-Ion Battery
,”
Adv. Energy Mater.
,
11
(
18
), p.
2003868
.
3.
Liu
,
B.
,
Jia
,
Y.
,
Yuan
,
C.
,
Wang
,
L.
,
Gao
,
X.
,
Yin
,
S.
, and
Xu
,
J.
,
2020
, “
Safety Issues and Mechanisms of Lithium-Ion Battery Cell Upon Mechanical Abusive Loading: A Review
,”
Energy Storage Mater.
,
24
, pp.
85
112
.
4.
Fu
,
Y.
,
Lu
,
S.
,
Shi
,
L.
,
Cheng
,
X.
, and
Zhang
,
H.
,
2018
, “
Ignition and Combustion Characteristics of Lithium Ion Batteries Under Low Atmospheric Pressure
,”
Energy
,
161
(
15
), pp.
38
45
.
5.
Yu
,
D.
,
Ren
,
D.
,
Dai
,
K.
,
Zhang
,
H.
,
Zhang
,
J.
,
Yang
,
B.
,
Ma
,
S.
,
Wang
,
X.
, and
You
,
Z.
,
2021
, “
Failure Mechanism and Predictive Model of Lithium-Ion Batteries Under Extremely High Transient Impact
,”
J. Energy Storage
,
43
, p.
103191
.
6.
Xu
,
J.
,
Liu
,
B.
,
Wang
,
L.
, and
Shang
,
S.
,
2015
, “
Dynamic Mechanical Integrity of Cylindrical Lithium-Ion Battery Cell Upon Crushing
,”
Eng. Fail. Anal.
,
53
, pp.
97
110
.
7.
Baker
,
D. R.
, and
Verbrugge
,
M. W.
,
2019
, “
Modeling Overcharge at Graphite Electrodes: Plating and Dissolution of Lithium
,”
J. Electrochem. Soc.
,
167
(
1
), p.
013504
.
8.
Lopez
,
C. F.
,
Jeevarajan
,
J. A.
, and
Mukherjee
,
P. P.
,
2015
, “
Characterization of Lithium-Ion Battery Thermal Abuse Behavior Using Experimental and Computational Analysis
,”
J. Electrochem. Soc.
,
162
(
10
), pp.
A2163
A2173
.
9.
Liu
,
B.
,
Duan
,
X.
,
Yuan
,
C.
,
Wang
,
L.
,
Li
,
J.
,
Finegan
,
D. P.
,
Feng
,
B.
, and
Xu
,
J.
,
2021
, “
Quantifying and Modeling of Stress-Driven Short-Circuits in Lithium-Ion Batteries in Electrified Vehicles
,”
J. Mater. Chem. A
,
9
(
10
), pp.
7102
7113
.
10.
Liu
,
L.
,
Feng
,
X.
,
Zhang
,
M.
,
Lu
,
L.
,
Han
,
X.
,
He
,
X.
, and
Ouyang
,
M.
,
2020
, “
Comparative Study on Substitute Triggering Approaches for Internal Short Circuit in Lithium-Ion Batteries
,”
Appl. Energy
,
259
(
1
), p.
114143
.
11.
Paul
,
N.
,
Wandt
,
J.
,
Seidlmayer
,
S.
,
Schebesta
,
S.
,
Mühlbauer
,
M. J.
,
Dolotko
,
O.
,
Gasteiger
,
H. A.
, and
Gilles
,
R.
,
2017
, “
Aging Behavior of Lithium Iron Phosphate Based 18650-Type Cells Studied by In Situ Neutron Diffraction
,”
J. Power Sources
,
345
(
31
), pp.
85
96
.
12.
Yang
,
X.-G.
,
Leng
,
Y.
,
Zhang
,
G.
,
Ge
,
S.
, and
Wang
,
C.-Y.
,
2017
, “
Modeling of Lithium Plating Induced Aging of Lithium-Ion Batteries: Transition From Linear to Nonlinear Aging
,”
J. Power Sources
,
360
(
31
), pp.
28
40
.
13.
Jia
,
Y.
,
Uddin
,
M.
,
Li
,
Y.
, and
Xu
,
J.
,
2020
, “
Thermal Runaway Propagation Behavior Within 18,650 Lithium-Ion Battery Packs: A Modeling Study
,”
J. Energy Storage
,
31
, p.
101668
.
14.
Liu
,
B.
,
Zhao
,
H.
,
Yu
,
H.
,
Li
,
J.
, and
Xu
,
J.
,
2017
, “
Multiphysics Computational Framework for Cylindrical Lithium-Ion Batteries Under Mechanical Abusive Loading
,”
Electrochim. Acta
,
256
(
1
), pp.
172
184
.
15.
Wang
,
L.
,
Yin
,
S.
, and
Xu
,
J.
,
2019
, “
A Detailed Computational Model for Cylindrical Lithium-Ion Batteries Under Mechanical Loading: From Cell Deformation to Short-Circuit Onset
,”
J. Power Sources
,
413
(
15
), pp.
284
292
.
16.
Xu
,
J.
,
Liu
,
B.
,
Wang
,
X.
, and
Hu
,
D.
,
2016
, “
Computational Model of 18650 Lithium-Ion Battery With Coupled Strain Rate and SOC Dependencies
,”
Appl. Energy
,
172
(
15
), pp.
180
189
.
17.
Jia
,
Y.
,
Gao
,
X.
,
Mouillet
,
J.-B.
,
Terrier
,
J.-M.
,
Lombard
,
P.
, and
Xu
,
J.
,
2021
, “
Effective Thermo-Electro-Mechanical Modeling Framework of Lithium-Ion Batteries Based on a Representative Volume Element Approach
,”
J. Energy Storage
,
33
, p.
102090
.
18.
Huang
,
Z.
,
Liu
,
J.
,
Zhai
,
H.
, and
Wang
,
Q.
,
2021
, “
Experimental Investigation on the Characteristics of Thermal Runaway and Its Propagation of Large-Format Lithium Ion Batteries Under Overcharging and Overheating Conditions
,”
Energy
,
233
(
15
), p.
121103
.
19.
Feng
,
X.
,
Ren
,
D.
,
He
,
X.
, and
Ouyang
,
M.
,
2020
, “
Mitigating Thermal Runaway of Lithium-Ion Batteries
,”
Joule
,
4
(
4
), pp.
743
770
.
20.
Jia
,
Y.
,
Yin
,
S.
,
Liu
,
B.
,
Zhao
,
H.
,
Yu
,
H.
,
Li
,
J.
, and
Xu
,
J.
,
2019
, “
Unlocking the Coupling Mechanical-Electrochemical Behavior of Lithium-Ion Battery Upon Dynamic Mechanical Loading
,”
Energy
,
166
(
1
), pp.
951
960
.
21.
Yuan
,
C.
,
Gao
,
X.
,
Wong
,
H. K.
,
Feng
,
B.
, and
Xu
,
J.
,
2019
, “
A Multiphysics Computational Framework for Cylindrical Battery Behavior Upon Mechanical Loading Based on LS-DYNA
,”
J. Electrochem. Soc.
,
166
(
6
), pp.
A1160
A1169
.
22.
Liu
,
B.
,
Jia
,
Y.
,
Li
,
J.
,
Yin
,
S.
,
Yuan
,
C.
,
Hu
,
Z.
,
Wang
,
L.
,
Li
,
Y.
, and
Xu
,
J.
,
2018
, “
Safety Issues Caused by Internal Short Circuits in Lithium-Ion Batteries
,”
J. Mater. Chem. A
,
6
(
43
), pp.
21475
21484
.
23.
Yin
,
S.
,
Hong
,
Z.
,
Hu
,
Z.
,
Liu
,
B.
,
Gao
,
X.
,
Li
,
Y.
, and
Xu
,
J.
,
2020
, “
Fabrication and Multiphysics Modeling of Modified Carbon Fiber as Structural Anodes for Lithium-Ion Batteries
,”
J. Power Sources
,
476
(
15
), p.
228532
.
24.
Jia
,
Y.
,
Liu
,
B.
,
Hong
,
Z.
,
Yin
,
S.
,
Finegan
,
D. P.
, and
Xu
,
J.
,
2020
, “
Safety Issues of Defective Lithium-Ion Batteries: Identification and Risk Evaluation
,”
J. Mater. Chem. A
,
8
(
25
), pp.
12472
12484
.
25.
Gao
,
X.
,
Lu
,
W.
, and
Xu
,
J.
,
2020
, “
Modeling Framework for Multiphysics-Multiscale Behavior of Si–C Composite Anode
,”
J. Power Sources
,
449
(
15
), p.
227501
.
26.
Wang
,
L.
,
Jia
,
Y.
, and
Xu
,
J.
,
2021
, “
Mechanistic Understanding of the Electrochemo-Dependent Mechanical Behaviors of Battery Anodes
,”
J. Power Sources
,
510
(
31
), p.
230428
.
27.
Gao
,
X.
,
He
,
P.
,
Ren
,
J.
, and
Xu
,
J.
,
2019
, “
Modeling of Contact Stress Among Compound Particles in High Energy Lithium-Ion Battery
,”
Energy Storage Mater.
,
18
, pp.
23
33
.
28.
Gao
,
X.
,
Lu
,
W.
, and
Xu
,
J.
,
2021
, “
Unlocking Multiphysics Design Guidelines on Si/C Composite Nanostructures for High-Energy-Density and Robust Lithium-Ion Battery Anode
,”
Nano Energy
,
81
, p.
105591
.
29.
Doyle
,
M.
,
Fuller
,
T. F.
, and
Newman
,
J.
,
1993
, “
Modeling of Galvanostatic Charge and Discharge of the Lithium/Polymer/Insertion Cell
,”
J. Electrochem. Soc.
,
140
(
6
), pp.
1526
1533
.
30.
Miranda
,
D.
,
Costa
,
C. M.
,
Almeida
,
A. M.
, and
Lanceros-Méndez
,
S.
,
2016
, “
Computer Simulations of the Influence of Geometry in the Performance of Conventional and Unconventional Lithium-Ion Batteries
,”
Appl. Energy
,
165
(
1
), pp.
318
328
.
31.
Nyman
,
A.
,
Zavalis
,
T. G.
,
Elger
,
R.
,
Behm
,
M.
, and
Lindbergh
,
G.
,
2010
, “
Analysis of the Polarization in a Li-Ion Battery Cell by Numerical Simulations
,”
J. Electrochem. Soc.
,
157
(
11
), p.
A1236
.