3.1. Model and Implementation
The model selection for this study included a lumped thermal model that included heat accumulation, heat transfer by radiation, heat transfer by convection, and heat generation, as seen in Equation (1). Heat transfer by conduction was neglected as the cell was suspended in the ARC, minimizing heat losses by means of conduction. This method did not entirely remove heat transfer by conduction; however, it was assumed that the heat transfer contribution by conduction was small enough that it could be neglected.
Here,
is the specific heat capacity of the test cell,
is the initial mass of the test cell,
is the modeled cell temperature,
is the heat generated by the test cell during the TR event,
is the surface area of the cell that is available for heat transfer,
is the emissivity of the cell,
is the Stefan–Boltzmann constant,
is the ambient system or ARC temperature, and
is the convective heat transfer coefficient. A lumped thermal model approach is not always an accurate assumption but is a sufficient assumption if the Biot number is less than 0.1, as shown in Equation (2).
In the equation, is the convection coefficient of heat transfer, is the characteristic length of the cell, and is the conduction coefficient of heat transfer of the cell. For experiments conducted in the ARC, the convection coefficient was small enough that the conduction coefficient was much larger, and, thus, the Biot number requirement to utilize a lumped thermal model was satisfied. To void the lumped thermal model assumption, the Biot number must be larger than 0.1, according to Equation (2). For an 18650 cell, this was not violated until the convection coefficient rose above , which was two orders of magnitude larger than that experienced in the ARC of , wherein the Biot number associated with the natural convection coefficient was 0.001. The calculation was not repeated for the 20 Ah pouch-style cell because the characteristic length of the pouch cell is smaller than the radius of the 18650 cell, which would result in an even larger convection coefficient being required to violate Equation (2). Because the lumped sum assumption is valid for this instance, it can reasonably be assumed that the temperature of the cell is not dependent on position since the entire cell will be approximately at the same temperature, and, therefore, the cell temperature is only dependent on time. Due to this situation, the partial differential expressed in Equation (1) can be written as an ordinary differential, .
It was understood in the literature that the majority of heat generated during a TR failure event originates from side reactions occurring between the SEI, positive and negative electrodes, and electrolyte, and only a small amount of relative energy is contributed by the stored electrochemical energy in the cell [
4,
5]. The heat generation is defined using Equation (3):
where
is the heat generated by the decomposition of the SEI,
is the heat generated by the reactions occurring between the positive electrode and electrolyte,
is the heat generated by the reactions occurring between the negative electrode and electrolyte,
is the heat generated by the decomposition of the electrolyte, and
is the heat generated from the internal short circuit during the TR event.
Despite the complexity of the aforementioned heat generation term, a simplified model was presented in this study, patterned after the model developed by Chen et al. where the total heat generation was separated into two heat generation terms and was regressed from empirical data [
7]. The relative distribution of heat generation among the various reaction sites in the cell will remain unknown; however, the total heat generation was determined and was used to accurately predict the temperature response of a TR failure event. The simplified heat generation term is defined in Equation (4).
is defined as the heat generated from the TR overheat event from the first critical temperature to just before the separator melting point, also defined as . is defined as the heat generated from the TR overheat event, from to the second critical temperature. There was an experimental difference between the first and second sections of the TR event. During stage , a side reaction began to occur, producing oxygen and causing irreparable damage to the cell; however, the cell was not self-heating at this point. Heat was generated from these reactions, but the reactions occurred slowly enough that the cell heating still needed to be subsidized for the temperature to continue to rise. In contrast, during stage , micropores in the separator began to develop, allowing shorting to occur and causing the temperature to rise much faster. The increased temperature accelerated the rates of side reactions occurring at the positive electrode, negative electrode, and SEI active sites. The combination of these two effects caused heat generation at this stage to be self-heating, thus showing the separation of stages and .
The calculation of
and
is defined in Equations (5)–(8).
In the equation,
(
i =
I and
II) is the mass of the reactant material,
is the enthalpy of the reaction,
is the modeled or predicted temperature,
is the normalized amount of reactant ranging from 0 to 1 (meaning that 100% of the reactant existing in the cell was present at the beginning of the TR event and 0% of the reactant consumed still existed at the end of the failure event),
is the regressed frequency factor for stage
,
is the apparent activation energy for stage
,
is the universal gas constant,
is the reaction conversion, defined as the inverse of
, and
is the reaction order, which was determined by Chen et al. to be a second order [
7]. During the TR event, we assume that the reaction conversion began at 0 and progressed to 1 when the maximum temperature was reached.
is the regressed frequency factor for stage
, and
is the apparent activation energy for stage
. The method for discrete calculation is presented in
Appendix B.
As mentioned previously, the model was performed in Macro-Enabled Microsoft Excel, and and were computed by initially assuming the time step just before the first critical temperature was reached; the reaction progression was zero, the reaction progressed proportionally to the temperature rise, and concluded at the maximum temperature. The initial conditions of , , and are given as , , and , respectively. and were implicitly solved simultaneously, using a Macro written in VBA.
Enthalpy is further defined in Equation (9):
where
is the thermal inertia used to compensate for the deviation from adiabatic conditions during TR failure events in the ARC,
is the specific heat capacity of the cell, and
is the adiabatic temperature rise. For stage I,
(indicating a near absolute adiabatic environment), and
. During stage
,
because the conversion of the reactants at that point is not complete in the model, and
[
7].
3.2. Experimental Results and Model Comparison
The initial stages of the study included initiating TR failure events by means of overheating in an EV+ ARC for a suite of Samsung 30Q 18650 cells. Empirical data were necessary to obtain the reaction kinetic parameters and calculate the total heat generation for the thermal model. The TR event profile for the Samsung 30Q 18650 cells under study is shown in
Figure 2a–c. Please note that the cell temperature depicted in
Figure 2a encompassed two distinct temperature measurements on the same cell during the thermal runaway event. This shows the deviation from adiabatic behavior seen during the thermal runaway event.
Empirical regression of the test data was derived by combining Equations (A1)–(A3) to produce Equation (10) (derivation in
Appendix A).
As discussed before,
is the temperature at which the separator began to deform and melt. This temperature separated the concavity of the TR profile curve. At the beginning of stage
,
was still approximately zero. There was some reaction conversion that occurred in stage
; however, it was very small compared to the reaction progression that occurred in stage
. Making this adjustment to Equation (10) and taking the natural log of both sides, Equation (11) was produced, which was used for the calculation of the kinetic parameters for the electrochemical portion of the model.
From Equation (11), TR empirical data was plotted as
vs.
. From this plot, shown in
Figure 3 and
Figure 4, a linear line was approximated wherein the frequency factor for stage
I and stage
II, respectively, was calculated with Equation (12), and the activation energy for stage
I and stage
II, respectively, was calculated with Equation (13).
Using
Figure 4a and Equations (11)–(13), the Stage
II kinetic parameters for the electrochemical portion of the model were determined for the Samsung 30Q 18650 cells and are presented in
Table 1.
Using these parameters and Equations (1) and (4)–(8), the TR failure event was modeled and compared with empirical data for model effectiveness. The thermal model results are shown in
Figure 5a–c. The model was shown to successfully predict maximum temperature and energy release within 1.98% relative error. It is important to note that as shown in
Figure 5c, the first critical temperature was accurately predicted by the TR model; however, the second critical temperature was predicted to occur 27.4 °C earlier than the measured second critical temperature. This relative error was attributed to a lack of empirical data points collected between the second critical temperature and the maximum TR temperature. Introducing large variations in temperature in small time steps caused instabilities in the model. To improve the critical temperature predictions, the model should be trained on empirical data with an increased time step; however, the critical temperature prediction was not greater than the measured critical temperatures. This was important so as to not predict an unsafe maximum operating temperature.
The TR failure cell-level testing was performed on a larger 20 Ah pouch-style cell in a similar manner to the testing performed for the Samsung 30Q 18650-style cells. The cell-level testing results are shown in
Figure 6a–c. Again, it is important to note that the cell temperature in
Figure 6a encompasses two distinct temperature measurements on the same cell during the TR failure event. The TR failure profile was similar to the Samsung 30Q 18650 style cell; however, the maximum temperature reached 502 °C with a first and second critical temperature of 108 °C and 236 °C, respectively. Using this empirical data, the aforementioned method for regressing and calculating the kinetic parameters for use in the electrochemical portion of the model was used for the 20 Ah pouch-style cell. The
vs.
plot is shown in
Figure 3b and
Figure 4b, along with the linear regression, equation, and coefficient of determination.
Using
Figure 3b and
Figure 4b and Equations (11)–(13), the stage II kinetic parameters for the electrochemical portion of the model were determined for the 20 Ah pouch-style cell and are presented in
Table 2.
From these kinetic parameters and Equations (1) and (4)–(8), the TR failure event was modeled and compared with empirical data to establish model effectiveness. The thermal model results are shown in
Figure 7a–c.
It is important to note the differences between the regressed pre-exponential factor and activation energy from the Samsung 30Q 18650 style cell and the 20 Ah pouch-style cell. As shown in
Table 1 and
Table 2, the apparent activation energy was fairly similar for both cells, with the pouch-style cell being 28.4 kJ/mol greater than the cylindrical-style cell; however, the regressed pre-exponential factor for the pouch-style cell was three orders of magnitude greater than for the cylindrical-style cell. Both of these increased parameters led to a greater reaction rate for the pouch-style cell compared to that of the cylindrical cell. There are a number of differences between the cell types that primarily lead to the increased reaction rates of cells. One of these differences, as stated previously, could be that the characteristic length of the pouch cell is smaller than the characteristic length of the Samsung 30Q 18650 cells. Since the pouch cell had greater capacity, a larger amount of electrode surface area was available for reacting with the electrolyte and the SEI and provided lower thermal resistance for heat transfer by conduction through the cell.
Similar to the prediction of the TR failure of the cylindrical cells, the model was shown to successfully predict maximum temperature and energy release within 1.34% relative error. As shown in
Figure 7c, the first critical temperature was accurately predicted by the TR model; however, the second critical temperature was predicted to occur 77.5 °C earlier than the measured second critical temperature. The model was more unstable for the larger 20 Ah cell than the previous cylindrical cell, primarily because the total TR duration was much shorter than for the cylindrical cells, which led to fewer empirical data points being collected during the TR event and a larger change in temperature with a shortened time step, leading the model to encounter some instability. It is estimated that this error can be remedied by increasing the rate at which data are collected during the TR for which the model was trained. Unfortunately, the ARC was limited in the rate at which data were recorded, and modifications would need to be made to the ARC or an alternate method would need to be used to collect data at a faster rate during the TR event. Although there was a large discrepancy between the measured and predicted second critical temperatures, the predicted critical temperature was less than the measured temperature, which allowed the user to be certain that an unsafe critical temperature would not be predicted by the model and, instead, a much more conservative critical temperature would be predicted.
The modes of heat dissipation were also monitored during the TR failure event for both cell formats and the results are shown in
Figure 8a,b. In both cases, the large majority of heat transfer before and after the TR event was due to radiation, and only a small portion of the heat dissipated was by convection. This metric will facilitate module and system design due to the large heat accumulation and TR propensity for propagation among cells within a module or system.
Figure 9 shows the ratio of radiation to convection heat transfer during a thermal runaway event as well as a function of the convection coefficient.
It was shown that at low convection coefficients, like those experienced in the ARC, the ratio of heat transfer during the TR event was larger than 1, indicating a greater amount of heat transfer by radiation than by convection. However, as the convection coefficient became larger, the ratio of heat transfer remained less than 1, indicating that the convective heat transfer was greater than the radiative heat transfer, even during the TR event. These figures show the need to increase the total amount of heat dissipation by methods of heat transfer other than radiation, such as convection by some sort of cooling system or by conduction through a heat sink or phase change material. The thermal model was performed using a convection coefficient of to determine what effect the convective heat transfer had on the maximum temperature of the TR failure event. Despite a much higher convection coefficient, the predicted maximum temperature did not change from that recorded in the tests with natural convection. This phenomenon is likely due to the rapid duration of the TR failure event. Despite the ratio of radiation to convection being less than 1 with a higher convection coefficient, the duration at which heat was allowed to dissipate at extreme temperatures was very small. The total amount of energy removed from the cell by means of convection was still small compared to that in the rest of the test. This brings into question the feasibility of increasing the convection coefficient to improve heat transfer through a module for the purpose of preventing or mitigating TR propagation. Although forced cooling could improve the module’s performance during normal operation, during a failure event, the maximum temperature will not be significantly affected by increasing the convection coefficient. This may be due to the extreme temperatures and rapid TR duration of TR propagation mitigation strategies, such as optimizing the case material in the cell and module billet to increase the emissivity; thus, the radiative heat transfer during a failure event will prove more successful in preventing propagation among cells than increasing forced convection heat transfer.