1 Introduction

Lithium-ion batteries (LIBs) are a popular type of rechargeable battery which have wide applications in portable devices such as cell phones, cameras, laptops, and even for vehicles and smart grids [1, 2]. However, fire safety issues remain a severe challenge for their further development. Safe storage and transport of these batteries with high energy density is becoming a concern to LIB manufacturers, transport industries as well as recycling plants, because of many fire incidents having been reported in the recent years. In the period spanning 1991 to 2018, 238 airport incidents related to Lithium-ion battery transportation as cargo or baggage have been reported [3]. Figure 1 shows two large fires that took place at LIB manufacturer and warehouses, causing serious damage. It is crucial to understand the fundamental mechanisms of the cause of these large-scale LIB storage fires to help develop strategies for preventing further accidents. Unfortunately, the causes of these fires are still not understood, with few published studies discussing LIB storage fires [4]. The commonly known causes of LIB fires such as electrical abuse, mechanical abuse and thermal abuse [5] are not applicable for these open-circuit stored LIBs, and many of these fires did not show signs of battery abuse. One possible explanation is the internal short circuits due to manufacturing defects [5, 6]. The problem is internal short circuits are almost impossible to verify because fire destroys any evidence.

Figure 1
figure 1

Two large LIB storage fires. (Left) a battery recycle plant in the UK caught fire causing the burning of 4 tonnes of waste LIBs [7], and (right) a LIB factory warehouse in China caught fire, the area burnt was approximately 1000 m2 [8]

While the electrochemical community seeks explanations mostly from the battery’s chemistry, one important factor that has been omitted is that heat transfer also plays a significant role in ignition [9, 10]. Ignition event occurs because the heat generation of the system exceeds the system’s ability to dissipate it. Fires could be triggered by a small heat generation with poor heat dissipation. Self-heating ignition [9, 10] is the spontaneous ignition of a material in an environment due to the exothermic reactions within the material [11]. It is a known cause of fires in many reactive materials such as coal [12], shale [13], and biomass [14] when stacked in large piles. The weak heat transfer within the large system enables the low heat generation from the spontaneous reactions to raise temperature. When the temperature reaches to a critical level, the system goes into a special state, Thermal runaway (TR), which the reaction rate of an exothermic reaction increases due to an increase of temperature leading to a further temperature increase and a further increase of reaction rate [6]. This positive feedback raises the temperature rapidly and triggers more violent reactions and finally flaming combustion happens. Driven by the low heat transfer, large stacked fuels could self-ignite at natural ambient temperatures.

Unfortunately, this kind of self-heating ignition has rarely been discussed for LIBs because self-heating of LIBs is often confused with another concept, TR propagation [15,16,17,18]. TR propagation studies when there is a failure on one single cell, how the fire propagates to the nearby cells. While self-heating ignition discusses when a large number of LIB cells are stacked, the weak heat transfer within the LIB ensemble might trigger a self-ignition. Most studies on thermal hazards of LIBs focusing on modelling chemistry [19,20,21,22,23,24,25,26,27,28] or experimental testing on TR fire behaviour [29, 30], failed to recognize that the low heat transfer for a large ensemble of LIBs may also contribute to triggering a fire. Our recent experimental study [31] showed that critical ambient temperature for self-ignition of LIBs decreases with the increase in the number of cells. While the classical self-heating theories assume continuum and homogeneous distribution of reactants, the real storage and transport of LIBs requires packaging materials to provide electrical insulation and mechanical cushion [32, 33]. These packaging materials might significantly change the self-heating behaviour of the system.

In this study, we build a numerical model to analyse the self-heating behaviour of a box of cells (100 cylindrical cells) in storage and try to investigate the effects of insulations on its self-heating behaviour. LiCoO2 (LCO) type of battery is chosen as a case study, because abundant kinetic data is available for this kind of batteries. We first compare the self-heating behaviour of a box of cells with one single cell and analyse the effects of the kinetics. Afterwards, the impacts of different insulations and packing configurations on self-heating behavior of the battery box are discussed.

2 Method

Semenov’s theory and Frank-Kamenetskii’s theory [9, 10, 34] are two classical theories which are mostly used for self-heating analyses of different materials. These two steady-state theories provide approximate analytical solutions to the relations between the geometrical size of the sample and the critical ambient temperature triggering self-heating ignition, and can be applied for a large-scale estimation. These theories are based on the assumptions of a one-step global reaction and an infinite amount of reactants. However, each LIB is a closed system with limited reactants, and multiple reactions have been reported to take place at elevated temperature [26, 27, 35,36,37,38]. Studies in the literature [10, 39] have already shown that for reactants with multiple reactions, extrapolation using Frank-Kamenetskii’s theory of high-temperature reactions to estimate self-heating behaviour at low temperature could result in wrong predictions. Therefore, the classical theories cannot necessarily be used to investigate the complex geometry introduced by insulation and its effects.

In this study, we adopt multi-step reaction kinetics [26, 27] and model heat transfer across the complex 3D geometry. The packing configuration of the box is based on a 18650 cell storage, as shown in Fig. 2. Variety types of LIB with different kinetics and thermal response are available in the market, it is impossible to investigate all of them. In this study, we choose the battery from Hatchard et al. [26]., E-One/Moli Energy ICR18650 (18 mm diameter, 65 mm length) 1.65 Ah cobalt cells in 100% state of charge (SOC) as the prototype and adopt their kinetics and parameters. In the box scenario, the distance between two cells is 20 mm, and there is a 5 mm gap between the front layer of cells and the cover of the box, as shown in Fig. 2a. The total dimensions of the box are 0.208 m × 0.208 m × 0.075 m. We assume all gaps are filled with air. In this packing configuration, the volume ratio of battery cells with respect to the whole box χ = Vb/Vbox is 0.51. Due to the symmetry of the geometry, one-eighth of the box is chosen as a computational domain to save computational costs, as shown in Fig. 1b. The numerical calculations are performed using the finite element method (FEM) software COMSOL Multiphysics 5.4.

Figure 2
figure 2

Domains of the battery boxes studied in this work. a The whole battery box containing 100 cylindrical cells, and b one-eighth of the box is chosen as a computational domain using geometrical symmetry

2.1 Heat Transfer

We simulate an oven test on a box containing 100 cylindrical LIB cells. The initial temperature of the box is T0, placed into a hot environment with temperature Ta. Inside the box, two materials are considered, the battery, where exothermic reactions take place, and the inert air. Heat conduction is considered as the main heat transfer. The 3D governing equation is the heat diffusion equation:

$$\rho_{i} c_{p, i} \frac{\partial T}{\partial t} = k_{i} \nabla^{2} T + q_{{{\text{tot}}, i}}^{'''}$$
(1)

where ρi is the density, cp, i is the heat capacity, ki is the thermal conductivity of i material, and \(q_{{{\text{tot}},i}}^{'''}\) is the volumetric heat source in i. Here, heat transfer properties ρ, cp, and k for i material are considered as constant effective values and do not change with temperature. In this box scenario, the air gaps are thin. The flow of air is constrained, so only heat conduction is considered and natural convection within the constrained narrow gaps is negligible. The heat transfer at interfaces between battery and air domain is calculated using the following equation

$$\left. {\rho_{b} c_{p, b} \frac{\partial T}{\partial t}} \right|_{{{\text{int}},{\text{b}}}} = \left. {\rho_{air} c_{p, air} \frac{\partial T}{\partial t}} \right|_{{{\text{int}}, {\text{air}}}}$$
(2)

At the free surface of the box, radiative and convective heat transfer are considered. Equation (3) is used for the convective heat flux component

$$q_{\text{conv}}^{''} = h\left( {T_{s} - T_{a} } \right)$$
(3)

where h is the convection heat transfer coefficient, Ts is the temperature of the free surface boundary, and Ta is the ambient temperature. Equation (4) is the radiative heat flux component at the boundary

$$q_{\text{rad}}^{''} = \varepsilon \sigma \left( {T_{s}^{4} - T_{a}^{4} } \right)$$
(4)

where ε is the surface emissivity, and σ is the Stefan-Boltzmann constant. No heat passes through the symmetric boundaries.

$$\left. {\rho_{i} c_{p, i} \frac{\partial T}{\partial t}} \right|_{\text{sym}} = 0$$
(5)

2.2 Chemical Kinetics

We adopt kinetics proposed by Hatchard et al. [26] and Kim et al. [27]. The four exothermic processes considered are: the solid electrolyte interphase (SEI) decomposition, negative-electrolyte reaction, positive-electrolyte reaction, and electrolyte decomposition. It is commonly believed that the thermal abuse reaction starts at the SEI layer, which is a thin passivating layer formed around the negative electrode. At around 100°C, the meta-stable component inside the SEI layer starts to decompose. While the SEI is decomposing, the intercalated Lithium inside the negative electrode is exposed to electrolyte and can react to form a new SEI layer. This reaction is called negative-electrolyte reaction. For the positive electrode, MacNeil et al. [37] found an autocatalytic mode of reaction for LiCoO2 and electrolyte. Meanwhile, the organic solvent of the electrolyte can also decompose at a relatively high temperature. All of these four reactions are exothermic and take place at different temperature ranges. The total volumetric heat source consists of four terms in Eq. (6).

$$q_{{{\text{tot}}, {\text{b}}}}^{'''} = q_{\text{sei}}^{'''} + q_{\text{n}}^{'''} + q_{\text{p}}^{'''} + q_{\text{e}}^{'''}$$
(6)

The detailed kinetics of each reaction are listed in Table 1. For a single cell, the internal structure of components is quite complex and has a heterogeneous structure requiring high computational resources to model local heating inside the cell. In this study, we are interested in the self-heating behaviour of the whole box. In this case, we assume that a battery cell is a continuous homogeneous material which is chemically and geometrically uniform. This assumption considers the whole heating effects of the cell as well as the heat transfer between cells. This simplification slightly decreases accuracy, but it greatly reduces the computational complexity. Kim et al. [27] compared their 3D model results with the lumped model result for one single cylindrical cell, and the results show small differences, which means this simplification is reasonable.

Table 1 The kinetic model for LiCoO2 Lithium-ion battery from the literature [26, 27]

2.3 Experimental Validation for a Single Cell

To validate our model and assumptions, simulations are performed first to study the thermal response of a cylindrical cell exposed to a preheated oven. The results are then compared with experimental results in the literature [26], as shown in Fig. 3. The values of the parameters used in the model are listed in Table 2 and Table 3. The numerical model gives a good prediction of the critical ambient temperature at which TR started to take place, but underestimates the temperature after ignition. It is because we only focus on the ignition stage. After the ignition, more complex combustion reactions occur, which is out of the scope of this study.

Figure 3
figure 3

Comparison of temperature history between our predictions (solid lines) and experiments (dot points) by Hatchard et al. [26]

Table 2 Physical and kinetic parameters used in the model
Table 3 Thermophysical properties of battery cells and insulating materials used in the model

2.4 Mesh Independence Analysis

The mesh used for the single-cell contains 384 elements. It is relatively easy to grid such a small homogeneous domain. However, for a complex heterogeneous multi-cells system, mesh quality might play a crucial role in the accuracy of calculations. It is necessary to conduct a mesh independence analysis on such a system. Three sets of meshes are used with 3872, 5150, and 7480 elements respectively. The temperature of the central cell and total heating power of the whole battery box are two crucial parameters and are chosen as the criteria. The result of the mesh independence analysis is listed in Table 4. It shows that all three meshes give almost the same results. Using the results of 7480 elements as a baseline, the meshes with 3872 elements and 7480 elements provide good accuracy. To ensure the accuracy for more complex simulation conditions, we chose the mesh with 5150 elements.

Table 4 Mesh independence analysis. The calculated results of the temperature of the central cell and total heating power of the whole battery box at the time of 2 h, 5 h, and 10 h based on three sets of meshes

3 Results and Discussion

A series of numerical calculations under different ambient temperatures are conducted to determine the critical ambient temperature Ta,cr that triggers self-heating ignition. All battery cells in the box are fully charged and in open circuit condition. The effects of chemical defects, insulating material, and packing configurations on the critical ambient temperature Ta,cr are analysed.

3.1 Critical Ambient Temperature

Figure 4 shows the temperature history of both the single-cell and the battery box exposed to different ambient temperatures. The central cell in the box is the local hot spot due to geometric symmetry, and its temperature is used to assess the onset of ignition. The Ta,cr for the battery box is predicted to be 125°C, 30°C lower than that of the single cell. This is because the heat dissipation condition for the central cell in the box is less effective than the single cell, which requires less heat generation and, therefore lower temperature to undergo TR.

Figure 4
figure 4

The calculated temperature history of the single-cell (left) and the box (right) under different oven temperatures. The vertical red line represents the time when TR starts to take place (Color figure online)

3.2 Time to TR

The numerical results show that although the battery box can self-ignite at lower ambient temperature, the time it needs to runaway is much longer. As Fig. 4 shows, under the critical ambient temperature (155°C for the single cell, and 125°C for the box), the time scale for the single cell is less than 100 min, while that for the battery box is over 10 h. In this study, the onset state of TR is defined as the time where the second derivative of temperature over time \(\frac{{d^{2} T}}{{dt^{2} }}\) the first time turns from negative to positive. After this critical time, not only does the temperature increase with time, but also the temperature increase rate \(\frac{dT}{dt}\) starts to increase, which leads to uncontrollable temperature increase and runaway. This time is ton, and the cell temperature at this time is defined as the onset cell temperature Ton. To be noted, Ton is the temperature of the central point of the system, while Ta, cr is the critical temperature of the ambient. In this case, the single-cell starts to thermally runaway at ton = 38 min, while the ton of the battery box is 7.6 h. The ton for both the single-cell and the box is marked with a red vertical line shown in Fig. 4. The ton is the criterion to judge the onset of TR. After this point, without any external measures being applied to cool down the system, the system can heat itself up, which means it already goes into a hazardous state. The ton is ahead of the time when there is a sharp temperature increase, which can be called fully developed TR state.

3.3 The Temperature Distribution Inside the Battery Box

Figure 5 shows the calculated temperature distribution inside the box at 11 h (Ta = 125°C). At this point, the box is already undergoing TR. The hot spot is located at the central point of the box because of its lowest heat dissipation due to geometric symmetry. There is an obvious temperature difference between different cells, while the temperature gradient inside one cell can be omitted. This is caused by the difference of thermophysical properties between the batteries and air. The batteries have much larger volumetric heat capacity ρcp than that of the air, which means that battery cells have a stronger ability to store energy and maintain their temperature. The temperature distribution of the central line (the black dashed line in Fig. 5) as a function of time is monitored to analyse the development of TR, as shown in Fig. 6. Due to the geometry symmetry, the temperature distribution in the y-direction is the same as that in the x-direction. The x-axis of Fig. 5 is normalised to the length of the box. In the early stage, the temperature of the outer layer of cells increases sharply. After 5 h, the battery cells in the central part of the box reach the environmental temperature. The heat generated by side reactions heats up the battery cells, leading the temperature of cells to be higher than the environment. After accumulating heat for another 2 h in the box, the temperature of the central cell reaches TR conditions.

Figure 5
figure 5

The calculated temperature distribution along middle cross-section planes inside the whole box (t = 11 h, Ta = 125°C). The dark dashed line is the line we used to analyse the temperature distribution history afterwards

Figure 6
figure 6

Predicted temperature distribution at the centre line in the x-direction at different times. The x-axis is normalised by the total length of the box (Ta = 125°C)

3.4 Chemical Reactions

To analyse the mechanisms that initiate TR, the dimensionless concentration of the reactants is studied, as shown in Fig. 7. All reactants show very similar trends for both cell and box. The meta-stable SEI component for both scenarios is already consumed by the time TR initiated, while the electrolyte remains unreacted because the temperature we focusing is lower than the temperature to initiate electrolyte decomposition. The intercalated Lithium in the negative electrode decreases slowly, while the positive electrode decomposes rapidly just after TR. To analyse the most important reactions triggering TR, the heat generation of each reaction is calculated. Heating power Qi is defined as the volume integration of the heat generation rate \(q_{i}^{'''}\) of i reaction.

Figure 7
figure 7

Predicted dimensionless concentration of reactants for the single-cell (left) and the box (right) at their critical ambient temperature. The vertical red line represents the time when TR is initiated

$$Q_{i} = \mathop \smallint \limits_{V}^{{}} q_{i}^{'''}$$
(16)

The predictions are shown in Fig. 8. The left plot shows the heating power of one single cell, while the right graph shows the heating power of the whole box with 100 cells. Qtot peaks at an early stage, around 0.2 h for the single-cell and 3 h for the battery box due to the SEI decomposition. TR is initiated at 38 min for the cell and 7.6 h for the box. After this critical time, the positive-electrolyte reaction is the dominant reaction to trigger TR for the single cell, while for the battery box, the positive-electrolyte reaction and the negative-electrolyte reaction contribute equally on triggering TR. It is because the battery box requires a lower temperature to initiate TR, which the reactions at low temperature contribute more heat to trigger TR. While the single-cell TR at a relatively high temperature that only the reactions at high temperature could initiate TR.

Figure 8
figure 8

Predicted heating power of reactions for the single-cell (left) and the box (right) at their critical ambient temperature. The vertical red line represents the time when TR starts to take place (Color figure online)

3.5 Effects of Activation Energy

We conduct a sensitivity analysis on 24 input parameters of the single-cell case. The results show that the activation energy of the positive reaction Ep is the most sensitive parameter for the single cell. In this case, we also evaluate the influence of Epe on Ta,cr for the whole box case. The Epe varies from 90 to 110% of the baseline from Table 2. The calculated results are shown in Fig. 9. It shows that Epe also has a large influence on the Ta,cr of the box of cells. With ± 10% variation of Epe, the Ta,cr could either decrease to 65°C or increase to 165°C for the box.

Figure 9
figure 9

Predicted Ta,cr of the battery box using different Epe (Basic value is taken from Hatchard et al. [26])

3.6 Effects of Packaging Materials

For the storage and transport of LIBs, electrical and mechanical insulation between cells is required to avoid short-circuits and mechanical collisions. Several kinds of materials are available with different thermophysical properties, which may change the self-heating behaviour of the box. In this study, we analyse three commonly used packaging materials: bubble wrap, polystyrene, and polyurethane. The thermophysical properties of these three materials are listed in Table 3. The predictions are shown in Fig. 10. Ta,cr of the box case with different packaging materials are estimated to be the same, 125°C, indicating that the type of insulating materials has a negligible influence on Ta,cr. This is mainly because the thermal conductivities k for all three packaging materials are similar compared to the thermal conductivity of the batteries, which is 2 orders of magnitude higher. This means the packaging materials introduced the same order of magnitude of heat resistance inside the box system, while the heat resistance by battery cells is much smaller and can be neglected. This makes the self-heating ignition of the box with these packaging materials almost identical.

Figure 10
figure 10

Predicted temperature history of the central cell with different insulating materials. The boxes with different packaging materials all TR at 125°C

3.7 Effects of Packing Configuration

There is no standard or guideline for the distance of gaps between adjacent cells for the package of LIBs. However, the industry tends to use smaller gaps to allow more cells to be stacked as long as the gaps provide reasonable protections for electrical insulation and mechanical cushion. The gap between cells could influence the critical condition for self-ignition. We have studied the effects of packing configuration with different spacing of cells. The same box with dimensions of 0.208 m × 0.208 m × 0.075 m is chosen for all packing configurations, with different numbers of cells inside (100, 49, 25 and 9 respectively), and the volume ratio of batteries over the total box χ = Vb/Vtot are 0.510, 0.250, 0.127, and 0.046 respectively, as shown in Fig. 11. While the former analysis has shown that the box with all three types of packaging materials has almost the same self-heating behaviour. We assume gaps between cells are filled with bubble wrap. The theoretic case, χ = 1, is considered as the extreme case for comparison. The predictions are shown in Fig. 12.

Figure 11
figure 11

The Sketch of domains for different packing configurations using the same box (0.208 m × 0.208 m × 0.075 m) with 100, 49, 25, and 9 cells. The volume ratio of battery cells χ = Vb/Vbox is 0.510, 0.250, 0.127, and 0.046 respectively

Figure 12
figure 12

Predicted Ta,cr for different packing configurations. The same box containing 9, 25, 49, 100, and 196 (theoretical) battery cells, with χ = 0.046, 0.127, 0.250, 0.510, and 1 respectively

Interestingly, Fig. 12 shows that the box with χ = 0.510 has lowest Ta,cr, which means it has the worst thermal stability and can self-ignite at the lowest ambient temperature. Compared with the theoretical condition, χ = 1, the Ta,cr for the scenario with χ = 0.510 is around 10°C lower, which means the presence of packaging material actually could significantly promote self-heating ignition. This is because the presence of thermal insulation significantly changes the effective thermophysical properties of the system.

For the theoretical case, χ = 1, it has the maximum number of battery cells as well as the highest heat generation. However, its ability to dissipate heat is also strongest because the heat conductivity of batteries is over 100 times higher than the packaging material. In this ideal scenario, heat resistance of the whole box is much lower than other scenarios, leading to a higher tolerance of the system against ambient temperature. On the other hand, the box with χ = 0.046 has the minimum number of battery cells as well as the lowest heat generation. However, it has the same order of magnitude of heat resistance as boxes with other χ. The lower heat generation with the same heat dissipation leads the box with χ = 0.046 a larger thermal tolerance and self-igniting at a higher temperature.

4 Conclusions

In this study, we build a model to analyse the self-heating ignition behaviour of a box of LiCoO2 (LCO) type of cylindrical cells under storage. We analyse the impacts of packaging materials and packing configurations during LIB package on the self-heating behaviour, and to provide insights for predicting the self-heating behaviour of large-scale LIB storage system. It is found that the critical ambient temperature triggering self-ignition of the battery box with 100 cells is 30°C lower than that of a single cell, indicating that the safety guidance based on tests on single cells does not guarantee the safety of a large pack of battery cells. The model predicts that the presence of packaging materials could accelerate self-heating ignition, while the type of materials used has a negligible influence on the self-heating ignition of the box as long as thermal conductivities of the insulation are the same order of magnitude of that of the air. It is predicted that the box with the volume ratio of battery cells in the box of around 0.51 has the lowest tendency to self-ignite. This study shows that the insulations during LIB package have crucial impacts on the self-heating behaviour of the LIB storage system. To make a more reasonable prediction on the self-heating behaviour of a large-scale LIB ensemble, the effects of insulations should not be omitted.

The results presented in this paper focus on LiCoO2 batteries with no defects: however, other types of LIB and batteries with defects would have different kinetics and could increase the tendency to self-ignite.