Characterizing and predicting 21700 NMC lithium-ion battery thermal runaway induced by nail penetration

,


Introduction
Lithium-ion batteries (LIBs) are increasingly used in electrified vehicles and energy storage systems [1]. However, LIBs are made of flammable/reactive materials. The electrolyte is not only flammable but also has relatively low flashing points. The materials for the electrodes are also reactive, e.g. the de-lithiated cathodes have high oxidability and lithiated graphite can react with electrolyte under elevated temperatures. The reactivities of these component materials render LIBs intrinsically hazardous under thermal, mechanical and electric abuse. Mechanical abuse, such as crash, penetration, and bend, can lead to the fracture of the separator or the deformation of the electrode, internal short circuit (ISC) and Joule heat from the penetration. Nail penetration has been widely employed to simulate the internal short process [2][3][4][5] due to mechanical abuse. It is also required by LIB safety standards in different countries [6].
A recent review paper of Liu and Jia [2] provided a good overview of the main evolutionary processes following nail penetration, i.e. mechanical deformation, ISC, thermal runaway (TR), and explosion/fire. Numerous experimental [4,[7][8][9][10][11][12] and numerical investigations have also been conducted to gain insight about LIB behaviour during nail penetration. Some of these will be briefly reviewed below in the order of experiments, modelling of single cells and cell blocks before introducing the present study.
Mao et al. [7] analysed nail penetration using 18,650 LIB cells for different states of charge (SOC), penetration positions, depths and speeds. They found that TR reaction was more severe when the battery was penetrated at the centre with faster penetration speed and the severity also increased with the SOC. Wang et al. [8] investigated 1Ah pouch cell for the nail induced ISC both experimentally and numerically. Wilke et al. [9] experimentally found that the use of phase change composite reduced the maximum temperature of the penetrated 18,650 LIBs by 60 0 C. Feng et al. [4] investigated the mechanisms of penetration induced TR propagation within a large format LIB module. More recently, Huang et al. [12] inserted a small temperature sensor in the cell to facilitate in-situ sensing in nail penetration tests. They observed multiple peaks of local temperature at ISC spot before TR, which could not be captured in conventional nail penetration tests. The findings also suggested that increasing the contact resistance or anode material resistance during ISC can effectively enhance LIB safety by delaying TR.
Some experiments were also conducted to characterise the electrochemical parameters associated with nail abuse, including electrical and thermal properties of individual material in the cell, short circuit current and cell terminal voltage, etc. It is, however, difficult to measure the short-circuit current in the experiments due to the fluctuations in cell terminal voltage associated with the underlying electrical, chemical and electrochemical processes induced by the penetration process. Wang et al. [9] and Zhao et al. [11] measured the voltage drop in penetration tests for some pouch cells. Huang et al. [12] performed nail penetration tests of 18,650 cells from the radial direction and measured the voltage in a series of tests under different SOCs, penetration speeds and depths. The terminal voltage of the cell and short-circuit current dropped significantly once the nail was penetrated [9,11,12]. The measurement of voltage is even more challenging for axial nail penetration in cylindrical cells due to connectivity between the voltage/current collector and the conducting nail.
To the best of our knowledge, measurements of the terminal voltage evolution with time for nail penetration induced TR in cylindrical 21700 LIBs have not yet been reported. In general, poor reproducibility of measurements is often experienced in laboratory tests due to uncontrolled shorting resistance and contact resistance during the process of nail penetration through the cell. This would render it challenging for detailed modelling of the electrochemical processes due to unavailability of the key electrochemical parameters.
In the meantime, the dynamic electrochemical-thermal process can also be captured by simplified electrochemical-thermal models with reasonable accuracy [13]. Most reported numerical simulations also neglected the nail induced cell deformation and only considered the Joule heat generated due to penetration [8,14]. Zhao et al. [13] developed a three-dimensional (3-D) multiscale electrochemical-thermal model for pouch cells. Their parametric study highlighted strong coupling of the cell thermal response and electrochemical behaviour, identifying key influencing parameters as shorting resistance, nail diameter and thermal conductivity. Zavalis et al. [15] developed a 2-D model for prismatic cells. Fang et al. [16] developed a 3-D electrochemical-thermal model to study the ISC for a 1 Ah pouch cell and compared the predictions with the measured temperature distributions during shorting. More recently, Liang et al. [14] proposed a 3-D electrochemical-thermal model by estimating short-circuit area equivalent resistance from experimental tests of pouch cells. Yamanaka et al. [17] considered nail movement as well as the heat release rates due to abuse induced chemical reactions. The predictions achieved reasonably good agreement with their own measurements; and indicated that the nail penetration speed had more influence on the combustion risk than penetration location. Cheng et al. [18] simulated nail penetration into a 5.25 Ah pouch cell by combining an electrochemical model with a thermal abuse reaction model. They considered the localized heat source at the nail using the experimentally determined contact resistance, which varied by 1-2 orders of magnitude. The results indicated that such fluctuation considerably influence the temperature rise. In particular, the heterogeneity in contact resistance led to 93% increase in temperature rise in comparison with homogeneous contact resistance.
Unlike single cells, relatively few simulations were conducted for TR propagation in cell groups. One example is the work of Kurzawski et al. [19], who developed a TR propagation model based on thermochemical source terms. They considered a stack of five LiCoO2-graphite 3 Ah pouch cells with and without metal plates placed between each cell to investigate cell-to-cell failure propagation where the end cell of a stack was forced into TR through a nail. Their predictions were consistent with measurements over a range of SOCs. The use of the metal plates between the cells were found to be an effective mitigation measure to increase system heat capacity and prevent TR propagation.
It should also be noted that most previous numerical simulations [8,[13][14][15][16][17][18] addressed nail penetration in pouch cells, for which the electrode and separator layers are stacked together. During penetration, the nail creates a hole in the electrodes and separator at a single piercing point. Hence, the current flow and flux distribution are relatively easy to model. A cylindrical 18,650 or 21,700 LIB cell consists of sheet-like anodes, separators, and cathodes that are sandwiched, rolled up into a "jelly roll structure" and packed into a cylindrical can. During penetration, the nail creates many holes in the same sheet, rendering the modelling efforts more challenging. Krishnadash et al. [20] numerically investigated an integrated design of phase change material and cooling plate to prevent the propagation of nail induced TR in cell groups of 18,650 LIB. Using a combination of Newman pseudo electrochemical model, short-circuit model and thermal abuse model, their 2-D simulations captured the intense heat generation. The predictions indicated that the integrated cooling approach kept the coolant temperature below its boiling point to avoid the undesired situation of coolant boiling. However, the focus of the study was on the mitigation effects. While the model includes one additional heat source due to the localized increase of resistance at the interface between the nail and the cell, the effect of penetration directions and locations were not considered.
In the present study, both numerical and experimental investigations are conducted for nail penetration through axial and radial directions of 21,700 cylindrical LIBs. The predictive tool is based on the the coupling between local joule heating at the nail site, heat addition by reaction kinetics in a lumped cell model and computational fluid dynamics (CFD) for heat dissipation. The ISC during penetration is represented by the maximum molar rate of lithium transfer per volume in terms of driving cell voltage.
The key innovation includes:.
-Pioneering CFD predictions of nail penetration induced TR in 21,700 cylindrical cells incorporating anisotropic thermal conductivities of the jelly roll stack as well as convective and radiative heat exchange with the surrounding. -Newly conducted tests to establish the different characteristics of TR evolution induced by radial and axial penetrations and provide data for model validation. -Parametric studies using the validated predictive tool to fill experimental gaps by quantifying the effects of ISC, contact resistance between the nail and cell materials, convective heat transfer coefficient and cell surface emissivity on TR evolution.

Numerical model
The predictive tool has been developed within the frame of in-house modified open-source CFD code OpenFOAM. The cell is treated as a solid region with anisotropic thermal conductivities while the surrounding environment or/and natural convective cooling are considered. The heat generation within the cell due to ISC and nail induced irreversible chemical reactions are calculated using the Arrhenius equations derived from the current accelerating rate calorimetry (ARC) tests and added as source terms for the energy equation.

Energy balance equation
The cell is represented by a solid structure formed with homogeneous layers. Some of the chemical and electrical energy stored in the cell were assumed to convert into thermal energy during TR evolution [22]. The ratio of how much electrical energy converts to heat is related to the thermal behaviour of the cell, in which heat conduction was assumed to be dominant. The jelly roll stack is homogenised and assumed to be one lumped block with anisotropic thermal conductivities while other properties are assumed to be uniform. As OpenFOAM uses Cartesian coordinates, conversion to local cylindrical coordinates is needed and the energy balance equation in cylindrical coordinates can be written as [23]:.
where, ρ is the density, kg/m 3 ; Cp is the specific heat capacity of the cell (J/kgK); λ R , λ ϕ and λ Z are the thermal conductivities (W/mK) in radial, cylindrical and axial directions. Q exo is heat due to the irreversible reaction of TR (W), Q conv is the heat source of the convection (W), Q rad is the heat source of the radiation (W) and Q cond is the source of heat conduction through nail (W). For computational efficiency, heat dissipation to the surrounding environment was accounted for by assuming constant convective and radiative heat sources. Convective heat flux to the cell surface was evaluated following the Newtons law of cooling:.
where h conv is convective heat transfer coefficient (W/m 2 K), A cell is the surface area of the cell (m 2 ), T cell is the surface temperature of the cell (K), T ∞ is the ambient temperature (K). The radiative heat flux out to the boundary was evaluated by:.
where, ε is the average effective emissivity and σ is the Stefan-Boltzmann constant.

Source terms for the energy equation
Thermal abuse reactions between the cell materials such as a positive electrode, negative electrode and electrolyte occur in various temperature ranges and release a great deal of heat and gas, which may lead to TR [24]. The volumetric heat generated can be evaluated by [21,25]:.
where, Q ISC is the exothermic heat due to ISC, Q SEI is the heat generation due to decomposition of SEI layer, Q An is the heat due to decomposition of Anode, Q Ca is the heat due to decomposition of Cathode and Q Mix is the heat due to decomposition of electrolyte and other remaining materials including the binder and it also incorporates heat generated by the earlier combustion during this last stage.

Heat generation due to ISC Q ISC
During short circuits, lithium intercalated in graphite, LiC 6 , represents a fuel that reacts electrochemically with metal oxide like CoO 2 and results into LiCoO 2 and C 6 products [13]. This releases a chemical energy internally with a non-Arrhenius rate that can be approximated by the following relationship between chemical composition and voltagedriven discharge:.
Here, W ISC denotes the density of the reactant, H ISC the reaction heat. The total energy for the reaction (product of W ISC and H ISC ) is calculated further in Equation (7), c C6,total is the molar concentration of graphitic C 6 sites in the battery volume, assuming a homogenized approach; X LiC6 is the mole fraction of lithiated C 6 sites in the graphite. X LiC6 was set as 1 for 100% SOC [19]. A ISC represents the maximum molar rate of lithium transfer per volume during a short circuit event. The rate of the short circuit is defined in terms of a driving voltage and a short circuit resistance that represents a nail penetration event. The maximum rate A ISC is inversely proportional to a short-circuit resistance [19]:.
where F is the Faraday constant, which is defined as the charge in coulombs (C) of 1 mol of electrons, 96,485C/mol and n C6,total is the moles of graphite in a cell. R ISC , is the short-circuit resistance. It must be calculated for each specific battery type and the event. In the present study, R ISC was assumed to be 0.005 Ω for the initial predictions. E max is the maximum voltage of the cell (3.7 V). A typical 21,700 NMC cell with 4.8 Ah capacity contains 10.4 gm of graphite (i.e. 0.8667 mol of n C6,total ) [26]. The energy for the reaction was calculated using the following equation [19]:.

Heat generation due to the abuse reactions
The Arrhenius equation for the decomposition reactions of SEI, anode, cathode and electrolyte i.e. Q i (i = SEI, An-Ele, Ca-Ele, and Mix) can be expressed as follows [27,28]:.
) c i and. dci dt = − κ i where W i denotes the density of the reactant, H i the reaction heat, κ i the rate constant of the reaction, c i the normalized amount of the reactant, A i the pre-exponential factor for the reaction, E i the activation energy for the reaction, and k b is the Boltzmann's constant. The term, W i H i, was calculated based on the change in temperature at different reaction stages expressed as W i H i = C p m cell ΔT i /V cell . Here, V cell is the volume of the cell. During the solution process, the heat generated by the four different reactions is divided equally into the computational volume of the cell as source terms for the energy equation. The predicted temperature field is used to feed further loop of the solutions represented by Eq (8). The quantity of heat being added to all the computational volumes is same for a single time step, but the cells at outer region are affected more by heat exchange with the surrounding than the internal cells, resulting in temperature gradients.
The temperature measurements in the ARC is combined with the average specific heat capacity of the cell to compute energy generation. The reaction kinetics for the reactions of SEI, anode, cathode and electrolyte were calibrated by comparing with the ARC measurements. During the self-heating stage, the temperature of the cell increases exponentially because of heat generation from abuse reactions which are exothermic and far greater than the heat that can be dissipated through the battery thermal management system. The temperature rise of the LIB under an adiabatic environment is determined by the following equation following Jhu et al. [29]:.
where, x is the degree of reactions, n is the reaction order. Eq. (10) is derived from Eq. (9) by taking natural logarithm on both and reducing to simplified form as follows,. Fig. 1 (a) depicts the currently measured temperature response and temperature rate profiles under the adiabatic environment in the ARC test for the present LIB following the same procedure in our previous publication [30]. From 132.45 • C onwards, the temperature rise accelerates irreversibly till the maximum cell surface temperature. This is treated as the onset temperature T onset of the self-heating stage. Plotting ln(dT/dt) versus 1/T as shown in Fig. 1 (b), the E i and A i are obtained from the slope and intercept of the fitting curve, respectively. The released heat during the TR stage is calculated as [29,31]:.
where, m battery is the mass of the sample LIB; Cp is the total heat capacity for the sample LIB, which is around 900 (J/kg⋅K) that measured by the EV-ARC, T TR and T max are the thermal runaway temperature which was 195.58 • C determined by the EV-ARC tests for the present cell and the maximum cell surface temperature during the TR (K), respectively.
It should be noted that there are some overlapping of the exothermic reactions in some temperature ranges. The exact temperature ranges for the decomposition of SEI layer, anode, cathode and the electrolyte and other remaining materials including binder are dependent on battery specifications [22]. Differential scanning calorimetry (DSC) tests [32] is needed to determine kinetics parameters of each separate exothermic reaction in Eq. (4) and (8). Due to limited information about the chemical compositions of the 21,700 cell, the activation energy and the reaction exponents, at the decomposition of SEI, anode and cathode were taken from references [27,28]. In these stages the decomposition of electrolyte and other materials was neglected as even if their decomposition already started at these stages, the heat generation rates were very small in comparison with that due to SEI, anode and cathode. The effect of change in the thickness of SEI layer was not considered here [27]. As noted by previous researchers [33], the relatively large amount of heat generated during the last stage was partially due to the ignition of released combustible gases and/or flammable electrolyte. For simplicity, the combined heat generation due to the decomposition reactions of electrolyte and other remaining materials including binder and any remains of the anode and cathode as well as some due to combustion of the released flammable gases was considered together as Q Mix in Eq. (4). As shown in Fig. 1 (a), the maximum cell surface temperature recorded in the tests were well above the autoignition temperature of the released flammable gases which range from 450 0 C for some hydrocarbons, 580 0 C for hydrogen and 609 0 C for carbon monoxide. It is likely that Q Mix also includes the contribution from combustion following ignition of the released flammable gases. The thermo kinetics parameters are estimated by fitting the EV-ARC measures temperature versus time. The normalized amount of all the reactants (c i in Eq. (8)) is assumed to be 1 at the initial step and subsequently their rates of change follow Eq. (8).
The heat produced by exothermic reactions due to the decompositions of the SEI layer, anode, and cathode under adiabatic condition is calculated per unit volume of 21,700 cells based on the available data for 18,650 NMC cells [31] due to similarities of the cell chemistry. The reaction kinetics are taken from the available literature for the SEI, anode and Cathode material [27,35]. The parameters for the exothermic reaction due to decomposition of electrolyte and other binding material are calculated from Fig. 1. Although such adoption may introduce some errors due to possible differences in the electrodes between the two types of cells, the influence is expected to be trivial due to the anticipated small changes in the materials for the electrodes. The heat produced by the remaining exothermic reaction and the earlier combustion which started during this stage is calculated from Eq. (11). The thermo reaction kinetics and related parameters used for the nail penetration abuse model are listed in. Table 1. The activation energy terms in the reactions due to the decomposition of SEI layer, anode, and cathode were taken from references [27,34]. During the last reaction of the electrolyte decomposition, a large amount of heat is generated due to the combustion of released combustible gases and/or flammable electrolyte [35]. The heat produced by exothermic reactions of SEI layer, anode, and cathode under adiabatic condition was calculated per unit volume of 21,700 cells based on the available data for 18,650 cells in [31]. Although such adoption may introduce some errors due to possible differences in the electrodes between the two types of cells, the influence is expected to be trivial due to the anticipated small changes in the materials for the electrodes. The heat produced by the remaining exothermic reaction and the earlier combustion which started during this stage was calculated from Eq. (11).

Heat generation at the nail-tip
At the penetration spot, the localized heat source was modelled by setting a heat source boundary condition considering the joule heating effect [8]:.
where,R n = R nail + R ct , and.R nail =

Lpen σ nail A nail
In Eq. (12), I nail is the current during nail penetration process. R n , R nail and R ct (Ω) are the shorting resistance at the nail site, the resistance of the penetrated nail and contact resistance, respectively. Here, nail site refers to the small area of the battery in contact with the nail following penetration. L pen (m) is the penetration length which is the length of the nail portion that is inside the cell. These are fixed as 10.05 mm in radial tests and 15 mm in axial tests according to the experiments to be described in the next section. A nail (m 2 ), and σ nail are the cross-sectional area and electrical conductivity of the penetrated nail. The shorting resistance at the nail site was determined by the intrinsic resistance of the nail material and the contact resistance considering imperfect contact between the penetrated nail and the cell components, as shown in Eq. (12) [13]. As the contact resistance is a complex function of nail diameter, penetration speed, cell and nail material, it is difficult to be quantified accurately. For the initial calculations, the contact resistance was assumed to be 100 m Ω. Subsequently, parametric studies will be conducted to examine the sensitivity of the predictions on the contact resistance.

The solution procedure
The OpenFOAM code is based on the finite volume approach and variable time steps can be used. To facilitate the present study, the heat generation within the cell due to internal short circuit (Eq.5) and decomposition of cell components (Eq. (8)) as well as the Joule heat generation within the nail tip (Eq.12) were added as source terms in the energy equation of the in-house modified version. The spatially averaged temperatures of all the discretized finite volumes in the solid-cell region were calculated at every time step and used to evaluate the volumetric heat generation rate by the exothermic abuse reactions. The evaluated volumetric heat was added to the cell centres of each finite volume and used in the solution of the energy equation for the solid region. This approach resulted in a significant reduction in the computational time in comparison with solving the decomposition reactions using the temperature value of each discretized finite volume and adding the heat into the energy equation.
The surface temperature is calculated in the middle of the cylindrical cell surface. The rapid change in the temperature during TR is captured using adaptive time step. The cell properties are approximated as homogenized, where the initial mass of the individual species, density, and specific heat are evaluated from a tear-down of the cell components as documented in our previous publication [30]. Anisotropic thermal conductivities are used in the local cylindrical coordinates for the jelly roll. The cell parameters used for numerical analysis are listed in Table 2.

Experiments
The 21,700 cylindrical LIBs were subjected to mechanical failure via nail penetration. The capacity and nominal voltage of the cell are 4800 mAh and 3.7 V, respectively. The cells were fully charged to 100% SOC before the test.
Two categories of nail penetration tests were performed using the facility at the UK Health and Safety Executive's Science and Research Centre. In Type R tests, the nail was penetrated in the radial direction while in Type A tests it was penetrated in the axial direction. The experimental setup of the nail penetration is shown in Fig. 2. In Type R tests, the cell was placed horizontally onto a piece of fireboard with its diameter shaped into it and slotted between two metal rods which prevented the cell from moving during penetration. The fireboard was supported by a metal bracket that could move along the y-axis of the base plate. The base plate consisted of a layer of aluminium and stainless steel. The cells position was adjusted to ensure that the centre of the cell was below the nail which was directly above as shown in Fig. 2 (a). In Type A tests, the cell was secured vertically onto two pieces of fireboard that made an L-shape using two cable ties. The vertical piece of the fireboard had the diameter of the cell shaped into it so that the cell sat flush. This was supported by a metal bracket that could move along the y-axis of the base plate. The top layer of the base plate was aluminium to provide heat protection, and the bottom layer was stainless steel to provide stability during failure. The cells position was adjusted to ensure that the centre of the cell was just offset from the nail which was directly above as shown in Fig. 2 (b).
The tests were performed under a constant penetration speed of ~ 70 mm/s and the depth of impact for radial tests was 10.05 mm till the central mandrel and 15.0 mm in the axial tests. The evolution of this process was captured using three high definition (HD) cameras (HD-SDI CCTV, NiteDevil) recording at 25 frames per second and one high-speed camera (Phantom Miro LC320s, AMETEK) recording at 1000 frames. One of the HD cameras was used to capture detailed evolution of the penetration process at close distance while the other two cameras were covering the wide angles to capture visible extent of the sparks and fire behaviour. The high-speed camera was set-up to record side-view. The  surface temperature of the cell was measured using a Type-K thermocouple. Measurements were taken at every 0.2 s during the tests and recorded by a data logger (FlexLogger, Version 2021 R1.1). A fan was turned on in the middle stage of the tests to exhaust the emitted smoke. A series of images were extracted from the video footages recorded by the cameras and shown in Fig. 3, in which frames (a-h) show the evolution from ignition to extinguishment during a Type R test while frames (i-p) show the evolution for a Type A test. During the initial stages shown in Fig. 3(a) and (i), sparks were observed to eject out from the cell venting caps. Due to ISC, large amounts of products were generated which led to the increase of the cell internal pressure. During the later stages shown in Fig. 3(b) and (j), more sparks sprouted out from the hole punched by the nail. Afterwards, the spark ejection gradually weakened to almost full stop. The existence of weakened sparks was observed during the starting of the flame in Fig. 3(c) and (k). The ignition of the released flammable gases due to the abuse decomposition reactions was likely triggered by the sparks. But the possibility of spontaneous ignition of the emitted gases was also possible. This was followed by the jet flame as shown in Fig. 3(d) and (l). The flame lasted approximately 4 s in the Type R tests and 18 s in the Type A tests. One jet flame was observed on top of the cell in the Type A tests. For the Type R tests due to radial penetration, the gases generated by the abuse reactions ejected from both the safety vent as well as the rupture of the side of the cell near the top. The co-existence of jet flames from both side rupture and safety vent led to the relatively short flame duration in the Type R tests. Fig. 3(e) and (m) show the flames at their maximum size in both tests. Burning of very low intensity flames were observed for approximately another 2 s in Fig. 3(g) and (o). This might be due to the remaining slow release of the combustible decomposition products from the electrolyte. The final frames in Fig. 3(h) and (p) show that the cells were red hot as the internal heat generation was much greater than heat dissipation to the surrounding. Table 3 summarises the results of the nail penetration tests, including the starting cell temperature at the beginning of the test, the maximum temperature during TR and cell mass before and after the tests. These are the digital times. Even in the presence of fire, the cell surface temperature decreases after TR following the completion of the exothermic decomposition reactions and the very limited heat feedback from the fire to the cell surface.

Model validation
The predicted and measured temperatures in all the tests are compared in Fig. 4. The predicted temperature profiles for the Type R tests are in a reasonably good agreement with the measurements shown in Fig. 4 (a). However, relatively large discrepancies are seen in Fig. 4 (b) between the predicted and measured values for the Type A tests. Heat generated due to the penetration of the nail is mainly due to the short circuiting of anode and cathode. It is transient in nature. Hyung et al. [36] experimentally observed that the cell voltage dropped significantly within 8 ms, and a violent explosion occurred 200 ms after the nail penetration, which led to further rapid rise in the cell temperature. It is also noted that while reasonable repeatability was achieved in the Type R tests, there are more variations in the measured cell surface temperatures in the Type A tests. This is partly because the penetration location on top of the cell in the axial direction in each test was not the same, resulting in variation of the specific jelly roll section being directly affected by the nail and the onset of the abuse reactions. Due to orientation of the "jelly roll", both the maximum cell surface temperature and the rising rate of the temperature were found to be higher in the radial (within 10 s) than the axial nail penetration case (within 20 s).
To help understanding the reasons behind the discrepancies, the measured temperature rise from "Tstart" to "Tmax" was plotted against the duration of the maximum temperature during the tests in Fig. 5. Here the duration was measured from starting of the nail entering the cell to the time when the cell reached its maximum temperature in the tests. As Table 3 Measurements of cell surface temperature in the nail penetration tests.  Fig. 4. Comparison between the predicted and measured cell surface temperatures. mentioned earlier, the type, material and speed of the nail as well as the depth of penetration was kept constant at 10.05 mm and 15.0 mm for the Type R and A tests, respectively. Lower maximum temperature rise was observed in the Type A( Fig. 5 (b)) than that of the Type R ( Fig. 5 (a)). This may be attributed to the fact that the actual short-circuiting due to nail penetration in the jelly roll was influenced by the orientation of the penetration. In the axial nail penetration, the effect of the short circuit was partly offset by the increasing surface area over which the current could flow, as evidenced by the X-ray radiography captured during the nail penetration tests for 18,650 cylindrical cells by Finegan et al. [37]. For the same reason, it also took longer for the cell surface to reach its maximum temperature in the Type A tests. This is consistent with the experimental observations in [37]. Although liner fits are included in Fig. 5 for both types of penetrations, this is only to indicate that the limited data points seem to indicate an almost liner reduction of the cell surface temperature with time. However, there is a relatively large variation in the temperature measurements and the results are hence not conclusive.
To analyse the effect of mass loss of the cell during the penetration tests, the maximum temperature versus the percentage mass loss, which was defined as the difference between the final mass and initial mass cell mass divided by the later in the experiment, are plotted in Fig. 6. No obvious pattern was noticed in both Type R and Type A tests. It seems in both types of tests, larger mass losses induced by the penetration would result in relatively lower temperature rise, but the decrease in the temperature rise was not monotonic. Also, relatively larger mass losses were observed in the Type A tests for relatively lower temperature rise. More intense sparks followed by flames were observed for shorter duration in the Type R tests as shown in Fig. 3. Fig. 7 shows the damage to the exterior of the cell during the tests. The red marked circle depicts the hole made by the nail. Although there were more passages for the sparks to eject out in the Type R tests through the nailed induced hole as well as the safety vent openings, the Type A tests experienced relatively higher mass loss as shown in Fig. 6. It took longer for the cell to reach its maximum temperature. Table 4 summarises the observed duration for maximum cell surface temperature, temperature rise and the percentage mass loss during the tests. Due to its short-circuiting characteristics and higher mass loss, the spark ejection followed by the flame was observed for longer duration in Type A tests, implying the presence of exothermic reactions for longer duration and higher mass loss than the radial case. Fig. 8 depicts the evolution of the predicted heat generation rate due to ISC and the decompositions of SEI, anode, cathode and electrolyte and other remaining material denoted asQ Mix . In both cases, all the reactions follow the same pattern for the heat generation rate because they are dominated by the weighted average temperature of the cell main body. Q Mix increased sharply at the time of the maximum cell surface temperature, indicating the start of combustion as discussed earlier. Fig. 9 shows the predicted total heat generated in the entire event of TR due to nail penetration in both tests. The heat generated is largely dominated by the internal short circuit. Further contribution of the cathode and electrolyte decomposition reactions including the contribution from combustion after the ignition of the released flammable gases are more in both tests comparing to that from SEI and anode decomposition. The heat generation inside the nail (calculated from the Eq (12), Q nail ), is very small compared to all these heat of reactions (Q ISC , Q SEI , Q Anode , Q Cathode , and Q Mix ). The calculated Q nail is 9.77 J in the in the penetrated nail area of 10.05 mm for the case of radial penetration and 14.66 J in the penetrated area of 15 mm in the axial case, respectively. Hence, they are neglected while comparing the heat of reactions within the cell due its own components. Fig. 10 shows the temperature distribution inside the cell before, during and after TR. The temperature inside the cell is higher than that on the cell surface until the cell temperature becomes the same as the ambient temperature. Fig. 10 (c) shows the radial distribution of the temperature on the top surface of the cell. A temperature difference of ~ 50 • K was observed in Fig. 10 (c). The higher temperature inside the cell is caused by the prevailing abuse reactions, which dominate inside the cell until the cell surface temperature is the same as the ambient temperature.
The temperature calculated at the top, middle on the cylindrical surface, bottom and nail tip of the 21,700 LIB cells is shown in Fig. 11. Temperature evolution during the process is highly dependent on the location. The nail tip shows a higher temperature compared to other locations during TR. During the cooling, a non-penetrated portion of the   nail acts as an extended surface and enhances the nail cooling. The maximum nail tip and cell top temperature were predicted as 1137 • C and 711 • C. This is similar to the measurements of Finegan et al. [37], which revealed a temperature difference of 310 • C between the nail tip and cell top maximum temperatures for the tested 18,650 LIB cells. Our predicted value is slightly higher as the 21,700 LIB cells which can be due a combination of differences in cell materials and cell internal structure. It should also be pointed that the current numerical simulation did not include the venting of the gases during TR, and hence the amount of energy that was expelled by the vented gases was not accounted for. This would lead to higher predicted maximum temperatures in the battery cell. However, the venting of the gases is not only influenced by the thermal decomposition but also the manufacturing process. Relatively large variation exists on the mass fraction of the ejected components for different LIB chemistry and even the same LIB in different tests. In this sense, numerical investigations without accounting for this factor can still lead to valuable insight on other processes. The underlying electrochemical and physical processes during nail penetration occur at different temporal and spatial scales. Strictly speaking, these electrochemical and physical parameters need to be determined for accurate simulation of the penetration induced ISC,   which is intrinsically a 3D phenomenon and influenced by the nail material, angle, speed and depth of the penetration and the nail diameter. In the present numerical simulations, the geometry of the current collector tabs and the spirally wound jelly roll inside the cells were incorporated in the model to determine the contact resistance and internal short circuit. The nail diameter was also considered to determine the total short-circuiting area of the penetration and the cell voltage discharge during nail penetration.

Parametric studies
Due to the constraints of the experiments, it is not always practical or even possible to investigate the effects of various parameters while numerical simulation has more flexibility. Following validation, the developed modelling approach was used to conduct parametric studies to examine the sensitivity of the predictions to some modelling parameters, which cannot be defined accurately. Additionally, whether in experiments or practical situations, some of these parameters are also varied from case to case. Understanding the sensitivity of the TR evolution to these parameters will hence provide valuable information about potential measures to mitigate the situation through battery thermal management systems. Again, more variation in the maximum cell surface temperature was observed for the different axial tests compared to radial tests. Hence, only radial penetration was considered in the parametric study using numerical analysis.

Effect of the contact resistance between the nail and cell
The contact resistance of the nail with the jelly roll may vary considerably in different penetration tests. The nail/bare metal current collector has the smallest contact resistance while the nail/pure active material film has the highest contact resistance [33]. As the details of the internal arrangement and some key properties are constrained to the battery manufacturer, a homogeneous contact resistance was used to analyse its effect on heat generation by calculating the cell surface temperatures during penetration. Fig. 12 (a) depicts the predicted temperature variations during the penetration assuming the contact resistance changing from 0.05 to 0.5 Ω. The results show that the maximum cell surface decrease by approximately 100 • C when the contact resistance increases from 0.05 to 0.5. Fig. 12 (b) shows the corresponding temperature variation in the nail tip. The contact resistance directly affects the heat generation in the nail tip (see Eq.10). The predicted temperature of the nail tip decreases from 1200 to 700 • C as the contact resistance increases from 0.05 to 0.5. The current flows through the nail decreases with the increase of the contact resistance, leading to the decrease of the Joule heating and accordingly lower temperature at the nail site due to lower heat generation and higher dissipation because of higher contact resistance.

Effect of internal short-circuiting resistance (R ISC ) of the cell
Internal short circuit occurs when the positive and negative tabs come into contact as the penetrated nail generates a great amount of heat in a very short time [11]. The short circuit resistance depends on shorting intensity. Apart from heat generation, it is challenging for experiments to provide more insight into the short-circuiting characteristics like this. The short-circuit resistance may vary with materials and contact conditions between the electrodes, current collectors, and the penetrated nail [38]. As the battery capacity increases, the internal resistances of the cell decreases [11]. Due to lack of information about the internal short circuit resistance in the 21,700 LIBs, parametric study was conducted to investigate its influences on TR evolution. Fig. 13 shows the effect of internal short circuit resistance (R ISC ) on the temperature evolution during TR induced by nail penetration. Considering Eq. (6), the smaller the R ISC the faster the reaction rate (A ISC ) and evolution to TR. The occurrence of ISC at the nail site results in large shorting current and Joule heat, which dominates the total heat generation in the cell, leading to fast evolution to TR. The induced Joule heat can only trigger LIB TR if the R ISC decreases to a substantially low level [39]. From Fig. 13, it is estimated that the threshold internal short circuit resistance for the current 21,700 LIB with 4800 mAh lies between 0.005 Ω and 0.006 Ω as TR was only triggered below the R ISC = 0.005 Ω. It is obvious from Eq.(6) that further reducing the resistance to R ISC = 0.001 Ω would accelerate the speed of evolution to TR. A similar inference has been given by Zhao et al. [11], whose experiments suggested that the threshold internal short circuit resistance for 5000 mAh pouch type cells to be in the range of 0.006 Ω to 0.01 Ω.

Effect of convective heat transfer coefficient
In LIB thermal management, efficient heat dissipation from the cell plays a vital role. Analysis is hence carried out on the effect of convective and radiative heat transfer to TR evolution. To study the effect of natural [11,13,19,40,41] and forced convection cooling [42][43][44] around the cell, a convective boundary condition is assumed at all boundaries of the LIB with constant convective coefficients (Eq.2). Fig. 14 [44]. The predicted surface temperature differs from 50 • C to 400 • C from natural convection to forced convection. This confirms the importance of using the appropriate convective heat transfer coefficient in numerical simulations. In the meantime, it also highlights the considerable cooling effect which can be achieved through enhanced convective cooling to dissipate the heat. Fig. 14 (c) zoom in on the surface temperature evolution in the first 10 s of the radial penetration, the temperature increase rapidly following nail penetration and the cell quickly enters TR in a few seconds. Subsequently, the cell surface temperature changes much slowly under the influence of external natural or forced convection.

Effect of surface emissivity
As shown earlier in Fig. 3, all tests were conducted without plastic wrap around the cell. Hence, typical values for the cell surface emissivity reported in the literature [11,21,27] are considered for the parametric study. No significant change is observed for the maximum cell surface temperature during TR. The role of emissivity only affects post-TR heat dissipation. The surface emissivity needs to be specified to simulate radiative heat losses (Eq.3) from the top, side and bottom surfaces of the cell. Fig. 15 shows the cell surface temperature calculated during the axial nail penetration using the surface emissivity of 0.65, 0.7, 0.75, 0.8 [27], 0.85, 0.9 [11,21] and 0.95. It can be seen from Fig. 15(a) that the temperature differs by up to 50 • C in the case of natural convection (h = 2 W/m 2 K) which reflects the dominance of emissivity. The effect of emissivity can be neglected in forced convection cooling (h = 71.7 W/ m 2 K), as depicted in Fig. 15 (b). These results also suggest that the influence of surface emissivity is less significant than the convective heat transfer coefficient.

Concluding remarks
Combined numerical and experimental studies have been conducted for TR evolution induced by radial and axial nail penetration into 4.8 Ah 21,700 cylindrical cells at 100% SOC. A predictive tool has been developed within the frame of OpenFOAM. The cell is treated as a solid region with anisotropic thermal conductivities while both convective and radiative heat exchange with the surrounding are considered. The heat generation within the cell due to ISC and irreversible chemical reactions are added as source terms in the energy equation. Measurements were conducted for the evolution of cell surface temperatures and the overall mass loss during the tests during both radial and axial penetrations. Snapshots from the video footages illustrate the ejection of sparks from the safety vent and nail induced holes as well as the subsequent flame evolution. The numerical predictions well captured the changing trends of the measured cell surface temperatures. Reasonably good agreement has been achieved between the predicted and measured cell surface temperatures for radial penetration, but relatively larger discrepancies exist for the axial cases. This was partly because the unrecorded variation of the penetration location from tests to tests and the associated variation of the specific jelly roll in contact with the nail.
The validated predictive tool has been used to conduct parametric studies to analyse the effects of internal short circuit resistance, contact resistance between the nail and cell materials, convective heat transfer coefficient and cell surface emissivity on the predictions. The results indicate that:.
-TR was only triggered if the internal short circuit resistance RISC is below 0.005 Ω, and further reduction of RISC to 0.001 Ω was shown to considerably accelerate the evolution to TR. -Heat dissipation through convective heat transfer to the surrounding has relatively little influence on the maximum cell surface temperature reached during TR for both radial and axial penetration cases.
The main influence of heat dissipation is post-TR cooling. The predicted surface temperature differs by as much as 350 • C from natural to forced convection, indicating the potential of forced convective cooling to prevent TR propagation in cell groups/modules. -The effect of cell surface emissivity is almost negligible in the case of forced convective cooling.
Although the predictive tool is developed and validated for a specific type of 21,700 LIB, the modelling approach is generic and can be easily adapted to other cell types.
Disclaimer. Aspects of the work described in this paper were undertaken at the Health and Safety Executive (HSE) Science and Research Centre. Its contents, including any opinions and/or conclusions expressed, are those of the authors alone and do not necessarily reflect HSE policy.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.