and Experimental Characterization of High Energy Density 21700 Lithium-ion Battery Fires

High energy density lithium-ion batteries (LIBs) are well suited for electrical vehicle applications to facilitate extended driving range. However, the associated fire hazards are of concern. Insight is required to aid the development of protective and mitigation measures. The present study is focused on 4.8 Ah 21700 cylindrical LiNi x Co y Mn z O (NMC) LIBs at 100% state of charge (SOC) with the aim to develop a viable predictive tool for simulating LIB fires, quantifying the heat release rate and temperature evolution during LIB thermal runaway (TR). To aid the model development and provide input parameters, thermal abuse tests were conducted in extended volume accelerating rate calorimetry (EV-ARC) and cone calorimetry. Some cells were instrumented with inserted temperature probe to facilitate in-situ measurements of both cell internal and surface temperatures. The mean peak values of the heat release rate, cell surface and internal temperatures were experimentally found to be 3.6 kW, 753 °C and 1080 °C, respectively. An analytical model has been developed to predict cell LIB internal pressure evolution following vent opening. The model uses the measured cell internal temperature and EV-ARC canister pressure as input data. Its predictions serve as boundary condition in the three-dimensional computational fluid dynamics (CFD) simulation of TR induced fire using opensource code OpenFOAM. The predicted transient heat release rate compare favourably with the measurements in the cone calorimetry tests. Predictions have also been conducted for an open cluster to assess the likelihood of TR propagation in the absence of cell side rupture. The present modelling approach can serve as a useful tool to assess the thermal and environment hazards of TR induced fires and aid design optimization of mitigation measures in enclosed cell clusters/modules


A Vent opening area (m 2 ) Subscripts
Aab Pre-exponential factor in the Arrhenius

Introduction
Lithium-ion batteries (LIB) are increasingly used in electric vehicles, consumer electronics and stationary energy storage devices. However, the individual components of LIBs are flammable in nature. High energy density batteries offer extended operating range/duration or increased power output, but their contents are more densely packed. The potential fire hazards associated with these batteries are of major safety concern [1][2][3][4]. LIBs are mostly safe within their designed usage environment. However, accidents occur when LIBs are exposed to abuse conditions [5,6], which can be divided into three categories: (1) thermal abuse due to overheating and fire, (2) electrical abuse due to overcharge, overdischarge and external/internal short circuit, and (3) mechanical abuse incurred by crash, penetration, and bend. These abuse conditions can trigger rapid exothermic reactions inside the LIB, which generate excessive heat and gases, leading to increases in LIB internal temperature and pressure. Once the designed critical pressure for the LIB safety vent is reached, the vent opens to release the generated gases. This is quickly followed by thermal runaway (TR), ejecting more gases generated from the abuse reactions. The exact compositions of the vented gases vary with LIB chemistry but typically consist of carbon dioxide (CO2), carbon monoxide (CO), hydrogen (H2), and some hydrocarbons. The flammable gases can auto ignite or be ignited by the hot sparks ejected from the LIB [7], resulting in fire and even explosion. Numerous incidents associated with LIB fire and explosion have been reported involving consumer electronic devices, electric vehicles, and aircraft [4]. It is hence of importance to gain insight of the flaming behaviour of high energy density LIBs to aid the development of protective measures and inform fire rescue services.
As summarised in a recent review [6], a large body of literature exists for 18650 cylindrical LIBs. In comparison, relatively few investigators have addressed the safety related issues of 21700 LIBs [7][8][9][10]. The authors' team [7] experimentally investigated 5 Ah 21700 LIBs under uniform heating with a flexible heater and used established correlations between heat release rate (HRR) and the mean flame height of turbulent jet flames to estimate the thermal hazards of the resulting fire. They also conducted further tests for the same LIB with flexible and nichrome wire heaters [8,9]. Some others reported modelling of 21700 LIBs when subjected to nail penetration or thermal abuse [10][11][12]. A brief review will be given about the previous modelling efforts in the following.
While significant progress has been achieved in the application of Computational Fluid Dynamics (CFD) techniques to process safety [13] [27] summarised experimental analysis from 20 tests and found that the major species include H2, CO, CO2 and hydrocarbons. The fractions of H2 and CO were found to increase with SOC while the fraction of CO2 decreased with the increase of SOC for some LIB chemistries, resulting in increase in fire hazards due to less inert gas in the vent gas mixture. Similar gas analysis was also performed by Golubkov et al. [28] for the vent gases of LCO/NMC/LFP. While each LIB type showed a unique gas composition, the main components identified were H2 and CO2 with significant amount of CO plus smaller fractions of methane (CH4), ethylene (C2H4), and ethane (C2H6). Due to lack of measurements for time-resolved vent gas compositions, the overall measurements of the major vented species concentrations were adopted as uniform vent gas compositions in almost the previous CFD simulations [6,19]. It should be noted that no gas composition measurements have been reported for 21700 LIBs, and none of the published vent models and CFD modelling of TR addressed this type of LIB.
One of the important parameters for defining fire hazard is the heat release rate (HRR), which gives information about the intensity at which the fire releases energy in the form of heat. The quantification of the HRR is essential for predicting the propagation of TR to neighbouring cells. Different experimental methods such as thermogravimetric analysis (TGA) [20], differential scanning calorimetry (DSC) [21], accelerating rate calorimetry (ARC) [22,23] were employed to study the energy release during TR while cone/in-situ calorimetry [24][25] were used to quantify HRR from the resulting fire.

The basic CFD framework
The simulations were performed using open source CFD code OpenFOAM [29,30].
The code has been previously validated for a wide range of jet and pool fires by the authors' group for both single and multiple component fuels. The RANS approach is used to solve the compressible flow with the two-equation model k-ω shear stress transport turbulence model of Menter [31]. Combustion is treated using the Eddy Dissipation Concept (EDC) [30], which assumes that chemical reactions occur in the fine structures and these structures are in the order of the Kolmogorov scales. The formulations are listed in the Supplementary Material.

Analytical model for determining the LIB internal pressure
The internal pressure of LIB is an important parameter which affects the vent gas flow rate. It is also needed as boundary conditions for the fire simulation. To develop a simplified model for the cell internal pressure, the LIB and EV-ARC pressure canister are considered as the system and surrounding. Measurements were conducted for LIB surface temperature, canister temperature and its internal temperature. The ejected gases were collected in the pressure canister filled with inert atmosphere of argon or nitrogen. Applying the Dalton's law of partial pressure, the chamber pressure can be written as the sum of the partial pressures of the initial inert gas and the partial pressure generated due to the releases of the vent gases from the cell.
Differentiating equation (1) with respect to time dP c dt = dP ar dt + dP vent dt Using ideal gas law, the above equation can be expanded to The partial pressure of the vent gas is affected by the change in the chamber temperature and the rate of change of the moles of the vent gas into the chamber, whereas the partial pressure of the initial argon gas is only affected by the change in the temperature of the chamber. In equation (2), the rate of change of the canister pressure ( dP c dt ) and the canister temperature ( dT c dt ) are measured experimentally.
The venting is choked until the EV-ARC canister pressure is greater than the critical pressure at the nozzle, i.e. P amb > P int ( 2 γ+1 ) γ γ−1 , where 'γ' is ratio of specific heat capacities.
Hence the mass flow rate of the vent gases is given by the following expression, Rewriting equation (3) in terms of the of vent gas moles, The internal pressure can be written as, where, T int is LIB internal temperature, C d is coefficient of discharge, MW vent is vent gas mixture molecular weight and 'A' is the vent opening area of the LIB. The term ( dnḋ t ) int represents the net rate of change of the moles inside the LIB, applying mass conversation The first term in the right hand side of equation (6) requires modelling to estimate the amount of gas generated during self-heating phase of the LIB while the second term can be evaluated from the ARC canister pressure measurements through the ideal gas law. Using the analogy of an auto catalytic reaction rate, the rate of thermal abuse and heat generation can be estimated with a one-step global expression [7,9] for the self-heating phase of the cell as where, 'α' is the degree of reaction progress, 'm o ' is initial mass of the LIB ,'R' is universal gas constant, 'A ab ' & 'E ab ' are the Arrhenius model parameters and ′H′ is the total heat generation for the unit mass. The rate of conversion of LIB mass is therefore given by equation (7). The rate of conversion of the LIB mass into heat can be considered analogues to the rate of generation of the gases due to the thermal abuse and can be modelled as in where, ' ' represents the volatile moles fraction and is a modelling parameter. Substituting equation (8) into equation (6) and subsequently into equation (5), the cell internal pressure (P int ) during TR evolution can be obtained.

Experiments
Experiments were conducted to characterise the TR characteristics of the present 21700 NMC LIB at 100% SOC and provide essential measurements for the vent model and energy release. These include three repeated cone calorimetry tests, EV-ARC tests with pressure measurements in the canister and two LIBs with inserted internal temperature probes. Brief descriptions of individual experiments are presented in this section with the detailed specifications of the tested LIB listed in Table 1.

Cone calorimetry tests
Cone calorimetry tests were conducted to characterise the energy released during TR. Prior to the tests, fresh batteries were firstly discharged and then charged with the same setup used in this work is shown in Figure 1.  The measured heat release rates in the repeated tests are shown in Figure 2 (a) and temperature measurements in Figure 2 (b-d). The plotted HRR values were calculated e based on the oxygen consumption principle from the oxygen measurements. As oxygen was generated by the decomposition of the cathode at elevated temperatures. This would result in measured oxygen consumption rate to be smaller than the actual value. To avoid the problem, oxygen generated from the decomposition of the cathode material was calculated to correct the HRR. The total heat release (THR) 'Qt' can be expressed as: where 'Q' is the heat measured from the calorimeter, 'Qs' is the heat from the combustion with self-generated oxygen which is released from cathode materials. The energy of oxygen consumed per unit mass is 13.1 kJ according to the principle of oxygen consumption method. The peak HRR and the total energy release obtained by integrating the HRR curves are shown in Table 2  relatively big difference is also noted between the " QS" in Tests 1 and 3. Additionally, the temperature in Test 2 experienced a sudden drop. Although the LIB and thermocouples were initially securely fastened to the workbench and LIB sides in the cone calorimeter tests, the violent movement of the LIB during TR evolution rendered some thermocouples to become dislodged from the intended location, resulting in lower temperature measurements.
As the temperature measurements between Tests 1 and 2 have relatively better repeatability, their average are used to deduce the input conditions for the CFD simulations.

EV-ARC tests with pressure canister
EV-ARC tests were conducted to measure the battery surface temperature, pressure of the canister and the overall compositions of the vent gases. The LIB was tested in a canister arrangement to measure any volume of gas given off during its decomposition as shown in Figure 3. The "heat-wait-seek" approach was used to provide information about the temperature at which the battery becomes unsafe as well as the quantity of thermal energy released by the battery in decomposition and the time it takes for the battery to fully decompose. The tests simulate a "worst-case" adiabatic environment that is similar to the conditions in the middle of a large battery pack where heat from a LIB cannot easily escape to the surroundings. The measured surface temperature and the chamber pressure are shown in Figure 4. The peak temperature of 786 °C and peak chamber pressure of 7.12 bar were measured in one of the EV-ARC tests. The volume of the pressure canister is 4.3 Litres. This is equivalent to a vent gas volume at standard temperature and pressure conditions of 273.15 K and 10 5 Pascals to be around 9717 cm 3 . Figure 4. The EV-ARC measured chamber pressures and LIB surface temperature.

Vent gas composition
The vent gases were collected from the canister connected to the EV-ARC to analyse the compositions of the major gas species. In Tests 1 and 2, the vented gases were collected in the canister using nitrogen (N2) as an inert gas while Test 3 used argon (Ar) as inert. This approach gives the overall amounts of gases generated rather than the transient values which need can be obtained from in-situ detection. The collected gases were sent for further analysis to quantify the compositions.

Battery Internal temperature measurements
To facilitate measurements of LIB internal temperature during TR, a temperature probe was inserted into the battery central void space avoiding contacts with any of the battery components. The LIB was then re-sealed to ensure leak proof even under elevated temperatures in the thermal abuse tests. The internal temperatures measured in the two repeated tests are plotted in Figure 5.The peak internal temperatures were 1109.3 °C and 1051.5 °C. The corresponding peak surface temperatures were 712 °C and 793.5 °C, respectively. Clearly, the temperature inside the LIB where the abuse reactions take place is much higher than that on the surface, where heat is dissipated to the environment. The measured temperature difference for the 21700 LIB is approximately ~ 300 °C in the radial direction during their peak values.

Numerical setup for the fire simulation
The computational domain for the fire simulation is shown in Figure 6   The vent cap of the 21700 LIB considered in the present study has three openings of ~1.5 mm each. During the tests, the vent cap was found to become completely dislodged and the central opening was around 10 mm in diameter as shown in Figure 6 (d), indicating the potential of improving vent design to improve LIB safety [41]. In some of the tests, the jelly roll also ejected out, leaving the central vent opening of greater than 10 mm. The present simulation only considered the gases ejecta, neglecting the eruption of the jelly roll, the modelling of which would not only requires the track of the internal pressure but also finite element analysis to capture the jelly roll deformation and collapsing. The vent opening size was assumed to be 10 mm following the experimental observation.  Table 3, the mean overall composition of the vent gases calculated from the three sets of measurements is used to represent the ejected gases.  [27]. Strictly speaking, the vent gas compositions and temperature vary with time during the venting processes. Because of the lack of measurements for real time variations, adopting an experimental measured overall vent gas composition is a common practice adopted by most modellers [7-9, 11, 33].

Boundary and initial conditions
As discussed in Section 2.1, the single step infinitely fast chemistry is used to determine the reaction rate in the EDC combustion model. The single step was developed based on global irreversible chemical reaction as a surrogate fuel for the vent gas mixture.
The thermophysical properties of the surrogate fuel were evaluated from the sum of the mole fraction ratios of the considered individual vent gas properties.  Figure 8.
It should be clarified that the analytical model predicts the cell internal pressure from the moment of vent opening. This was also the time when the canister started collecting the vented gases (Figure 7a). This is unlike other published measurement or predictions which considered the pressure evolution in the canister or inside the cell for the whole duration of TR evolution [33][34][35][36][44][45]. Table 4. The Arrhenius model parameters for equation (7) Figure 8. The cell internal pressure evolution predicted by the analytical model.
The analytical model predicts the internal pressure from the moment of vent opening.
It assumes that the vent opens when the internal pressure reaches the critical vent opening pressure, which is set by the manufacturer as 27.9 bar as shown in  Table 4 are model constants for Equation (7), obtained by log curve fitting the EV-ARC LIB temperatures profiles following the same steps described in our previous work [8][9]. The predicted cell internal pressure evolution in Figure 8 shows two pressure peaks. The first peak is the critical vent opening pressure and the second peak is due to TR, which is predicted to be around 4.56 bars for the present 21700 NMC LIB.
Although no experimental measurements are available for direct comparison, the occurrence of the second peak during TR is in line with the experimental measurements for NCA 18650 LIB by Ziebert [34] and Lei et al. [35], which detected the second internal pressure peak to

Parameter Values
be around 5 bar. Similar value and trend were also predicted for 18650 LIB by the recent prediction of Kong et al. [19]. However, Kim et al. [11] and Mao et al. [33] predicted the second internal pressure peak to be around 2.5 bar and 1.056 bar, respectively for 18650 LIB. during the TR. The existence of the second internal pressure peak at TR was also experimentally found by Wang et al. [42] for 50 Ah prismatic cells and Zhao et al. [43] who measured the canister pressure for 18650 LIB. However, the inconsistency in the measured/predicted values of the second internal pressure peaks by different authors for 18650 LIBs suggests that the actual predicted and measured values of the internal pressure peak at TR should be considered with caution for both 18650 and 21700 LIBs in the lack of validation.
As the LIB internal and external ambient pressure ratio is higher than the critical value for choked flows, which is typically 0.58 for diatomic gases. To avoid simulating the shock-laden flow at the nozzle exit of an under expanded jet flow, the pseudo diameter approach is used. It is physics-based semi-empirical model derived from the conservation of mass and momentum to calculate the jet conditions across the pseudo diameter as sonic velocity at ambient pressure. In the present study, the pseudo diameter approach of Birch [32] is adopted to evaluate the vent gas velocity at the pseudo inlet which is the opened vent. The ARC canister temperature and pressure are considered as the ambient condition.
The gases were assumed to be venting at the experimentally measured cell internal temperature. The calculated pseudo diameter, which serves as the inlet boundary of the fire, changes with time due to the continuous variation of the LIB internal pressure. Once the internal pressure reaches equilibrium with the ambient pressure, the Pseudo diameter is the same as the vent opening diameter.

Results
The transient CFD predictions are presented for the temperature at various heights above the vent as well as the HRR and temperature at the top of the adjacent cell to evaluate the likelihood of TR propagation in LIB clusters.

Grid sensitivity study
The RANS CFD approach is known to be sensitive to mesh resolution. To establish that the numerical predictions are independent of the mesh size, preliminary simulations were performed with three different mesh sizes, i.e. 3.5, 4.5 and 5.0 million. The predicted flame temperature variations with time at 30 cm height are compared in Figure 9. It can be discerned that the differences in the predictions of the medium and fine meshes are very small. Therefore, the final simulations were performed with the medium resolution as a compromise between computational efficiency and accuracy. Figure 9. Mesh sensitivity study of flame temperature at centreline (30 cm ).

Battery fire dynamics
The underlying complex chemical and physical processes in TR involve decomposition of the electrode materials, burning of the electrolyte as well as ejection of sparks and flammable gases and the subsequent fire [7]. The whole process can be divided into soft venting, initial sparking and flaming.
(c) heater to trigger TR following the same procedure as described in our previous paper [7].
The ejected sparks can act as ignition sources for the vented flammable gases. But it might also be possible that these gases auto ignited as the measured internal temperatures, as shown in Figure 7(c), which are well above the auto-ignition temperature of the flammable gaseous constituents. The authors' group previously analysed [7] the high speed video files to establish these distinct phases mentioned afore and proposed a simplified approach to estimate the heat release rate from the resulting jets fires using established correlations for jet fire heights [36].   Figure 10. Some quantitative comparison between the predicted and measured HRR will be discussed in the next sub-section.

The heat release rate
The average values of the measured HRR and mid-height cell surface temperature from Tests 1 & 2 in the cone calorimeter are plotted in Figure 12. Two important observations can be made. Firstly, there is a slight time delay between the LIB surface temperature and HRR peaks. The duration of the stable LIB fire can be considered as having the same time span as the HRR curve. The slight time lag between the temperature peak and the HRR peak can be attributed to the ignition delay from the initial sparking stage.
Although the hot sparks carry energy and can transfer it to the surrounding cells, such energy release cannot be captured by the oxygen consumption method. This is indeed a limitation in using such method to quantify the energy release induced by LIB TR. In addition, some deviations can also occur due to oxygen generation in the decomposition of the cathode/anode/electrolyte materials itself, which is corrected in the current experimental analysis and found to be negligible by previous research. Huang et al. [38] compared the theoretical heat release rate from normalized electrolyte molecules with those obtained from the cone calorimeter tests and found that the HRR curves predicted by the two approaches had similar values and trends. This implies that the oxygen generated during the early stages of the decomposition reactions were ejected out before combustion even started. Figure 12. The experimentally measured heat release rate and cell surface temperature at mid-height in the cone calorimeter. Figure 13. Comparison between the predicted and measured heat release rate.
The predicted and measured HRR in the cone calorimeter are compared in Figure   13. The predictions match well with the measurements up to the peak HRR when the volatile moles fraction " " was set to 0.1125 in the analytical model for the internal pressure evolution. The predicted HRR deviate more from the measurements after the peak. This does illustrate the limitations of the present simulation but is of relatively less concern as the main purpose of the simulation was to capture the heat release rate during TR resulted fire at its peak. Quantitative comparison between predicted and measured LIB fires will be hard to achieve. Firstly, as discussed earlier no two LIB fires are the same. There is a great degree of stochasticity in the cell internal structure change due to material deformation and decomposition. Side ruptures can also occur and the actual size of the opening on the cell top is not limited to the vents and vary from cell to cell as shown in Figure 6  Radiative heat transfer plays a very important role in fires and its hazard for process safety. In many situations, radiation is the dominant mode of heat transfer in the vicinity of a fire source, especially for large-scale fires. For enclosed LIB clusters, radiative heat transfer is an important contributing factor for heating abuse of the neighbouring cells. To evaluate the effect of different radiative heat transfer treatment, predictions were performed using both the P1 and the more accurate finite volume based discrete ordinate method (fvDOM).
The P1 is a simplified approach but computationally more efficient than the more accurate fvDOM. In Figure 15,   The predicted LIB surface temperature is an indication whether the adjacent cell will enter TR based on the modelled critical self-heating temperature limit. Once the adjacent cell top temperature exceeds the critical temperature, the self-heating due to solid electrolyte interphase decomposition would be triggered. CFD simulations have also been conducted for an open cluster as shown in Figure 15. The central LIB was assumed to undergo TR and represented by the circular patch of 10 mm in diameter on the central LIB top as inlet. The predicted temperature on the top of the adjacent LIBs with the two different radiative heat transfer treatments are shown in Figure 16. The results indicate that the TR and resulting fire would not trigger TR propagation in an open cluster, i.e. without casing. However, it should be noted that the numerical simulations here did not consider LIB side rupture and other electric connections, which may affect the TR propagation dynamics. Even in the absence of LIB side rupture, the influencing factors will be very different in enclosed/encased LIB clusters. In such situations, the ejected hot gases/sparks and the fire will impinge on the casing and spread into the narrow gaps, leading to enhanced transfer of heat energy to the adjacent ones as well as those further away [38,39]. Work is ongoing to use the developed modelling approaches to analyse TR propagation in enclosed LIB clusters. In terms of LIB safety, it is also worth consider the fire hazards in conjunction with the integrity of LIB modules, etc. For example, whether it is viable to use weak casing for some applications so the fire can easily penetrate rather than spreading inside the module to trigger TR in the neighbouring cells, leading to much large fires.

Conclusions
Numerical and experimental studies have been conducted for high energy density 4.8 Ah 21700 cylindrical NMC LIBs at 100 % SOC to characterise the TR evolution and resulting fire. The kiey findings can be summarised as follows: • The LIB was characterised under thermal abuse conditions in EV-ARC and cone calorimetry tests. The mean peak values of the heat release rate, cell surface and internal temperatures were experimentally found to be 3.6 kW, 753 °C and 1080 °C, respectively.
• An analytical model has been developed to predict the cell internal pressure from vent opening onwards. The model uses the measured cell internal temperature, canister pressure and overall gas generation in EV-ARC tests as input data. The predicted cell internal pressure profile following venting is inline with literature findings for similar cell types. The calculated internal pressure had a peak of value 4.56 bar during TR.
• The predicted cell internal pressure and measured temperature along with the mean gas compositions were used as boundary conditions for the subsequent 3-D CFD simulations of the TR induced fire. The predicted heat release rate achieved reasonably good agreement with the cone calorimeter measurements.
• Predictions were also conducted for open cell clusters, revealing that TR is unlikely to propagate to the neighbouring LIBs in an open cluster in the absence of LIB side rupture.
It should, however, be noted that the situation would be different for enclosed clusters and with electric connections, for which the impingement of the ejected hot gases/sparks on the casing and their direct contact with other LIBs will augment the impact of the thermal hazards.
The present CFD modelling approach can serve as a useful tool to assess the thermal and environment hazards due to LIB fires and aid design optimization of mitigation measures in practical enclosed LIB clusters/modules.