A Comprehensive Model and Experimental Investigation of Venting Dynamics and Mass Loss in Lithium-ion Batteries under Thermal Runaway

1. Introduction

1.1. Motivation

In the new era of the highly developed society, the demand for energy has been increasing with the emergence of industrialization and the integration of information technology. A global shift towards more sustainable and environmentally friendly energy solutions is ongoing, particularly in the area of energy storage.

Lithium-ion batteries (LIBs) have appeared as a key technology to meet these demands due to their high energy density, efficiency, long cycle life, and environmental compatibility. They have already become the backbone of many electronic devices and are increasingly being used in electric vehicles (EVs) and grid-scale energy storage systems, contributing significantly to the reduction of greenhouse gas emissions. Despite their widespread adoption, LIBs have also brought safety issues such as thermal runaway (TR), which can lead to serious safety risks such as fires and explosions. [1,2,3,4,5,6]

Thermal runaway is a critical phenomenon where a battery cell undergoes an uncontrollable increase in temperature, triggering a cascade of exothermic reactions that rapidly liberate substantial energy in a short duration. This can result in a significant increase in temperature and pressure, which can lead to an explosion or fire [1,2,3,4,7,8,9,10,11,12,13,14,15,16,17,18]. During TR, the increase in battery temperature not only causes electrolyte evaporation [19], but also triggers reactions that produce a large volume of gases [20]. This combination can result in a sharp increase in internal pressure, leading to expansion of the battery casing.

The venting process is a critical safety feature designed to mitigate the pressure build-up within a battery cell during TR [19,21]. As the internal pressure reaches a critical threshold in a battery equipped with a safety vent valve, this valve opens, releasing a mixture of electrolyte vapors, reactive gases, and solid particles into the surrounding environment [22]. In the absence of the vent valve, a feature often omitted in pouch cell designs, the mixture typically leaks and bursts from the sealant edges of the battery [21]. Together with the previous evaporation stage, this whole process consumes energy and partly cools the battery [19]. However, it often cannot stop the process of TR, where the cell temperature continues to rise steeply. The high temperature and substantial proportion of flammable gas in the resulting mixture often lead to a second venting eruption in the form of explosions and fires, appearing as jet flames [23]. The cell loses the most material during this process and shows a peak of its temperature. After consuming or losing the major mass of its active material, the TR reaches an end, where the battery is cooled by the ambient environment.

To optimize the cell and battery design in terms of safety and to gain a further understanding, models are useful tools. Therefore, accurate modeling of the venting process and mass loss is crucial for predicting the behavior of LIBs during TR and for developing strategies to lower the risks of causing greater hazard.

1.2. Relevant Literature

Several studies have made significant contributions to the understanding and prediction of TR [2,7,8,9,10,11,12,13,14,15,16,17,18,24,25,26,27,28] and venting in LIBs [5,7,19,21,23,29,30,31,32,33,34,35]. A review of the literature related to the modeling of TR and venting with their purpose and researched parameters is shown in Table 1. Researchers first focused on revealing the mechanism of TR by experimental studies and numerical modeling. Richard and Dahn [8,9] investigated the thermal stability of lithium-ion cells, focusing on lithiated mesocarbon microbead (MCMB) material in various electrolytes under adiabatic conditions. Utilizing an accelerated rate calorimetry (ARC), the study explored the effects of lithium content, electrode surface area, electrolyte type, and initial heating temperature on the thermal stability of these materials. Based on the acquired data, the exothermic reaction kinetics of the solid electrolyte interphase (SEI) layer’s decomposition and regeneration during TR on the negative electrode has been characterized, leading to the development of a simple lumped model to predict TR events. Following that, Hatchard et al. [12] have introduced a more comprehensive one-dimensional thermal model incorporating the TR reactions in the positive electrode, which was based on oven exposure testing and was used for assessing the effectiveness of cell designs in heat dissipation. Kim et al. [13] marked a significant advancement in TR modeling by moving beyond one-dimensional models to a comprehensive three-dimensional model, and integrating electrolyte and binder decomposition. Their study offered detailed formulations for several potential exothermic reactions within the battery with geometrical features, enabling simulations that reveal how heat propagates through a cell during TR. These research works have mainly focused on the heat generation and temperature change of the battery during TR process, providing a primary understanding to the TR mechanisms.

With the deepened study of the TR process, more attention has been paid to more detailed failures from TR. During TR, the rising temperature of the battery not only causes the evaporation of the electrolyte [19] but also triggers reactions that produce a large volume of flammable or toxic gases, which contributes to the increase of the internal pressure of the cell and can lead to combustion or explosion under higher temperature [20]. This mechanism has been extensively studied by experiments and quantification work by various researchers [20,36,37,38]. To first determine the composition of the main generated gases from TR, Golubkov et al. [20] conducted a series of TR experiments on 18650 LIBs using various cathode materials: a combination of Nickel Manganese Cobalt Oxide (NMC) and Lithium Cobalt Oxide (LCO), NMC, and Lithium Iron Phosphate (LFP), and quantified the gases with gas chromatography (GC) too. Subsequently, they further explored the TR characteristics of two commercially available types of 18650 LFP and Lithium Nickel Cobalt Aluminium Oxides (NCA) cells and estimated the mass inventory based on teardown results [36]. Shen et al. designed experiments using an adiabatic explosion chamber (AEC) under an inert atmosphere to test prismatic LIBs with LFP and different NMC cathode materials [37]. Besides the analysis of gas components, they also conducted the calculation of flammability limit based on the generation of ${\mathrm{H}}_2$, and thus evaluated the hazard degree of TR from LIBs with different cathodes. These experimental studies have clarified the types and quantities of gases that can be generated during TR in specific types of LIBs, thereby enabling advanced modeling of TR events.

Coman et al. first developed a lumped model for TR in 18650 LIBs, taking into account the evaporation of the electrolyte and the ejection of contents from the jelly roll [19]. They provided equations to describe the development of internal pressure and venting but also pointed out the difficulty of describing the mass loss directly with equations because of its unpredictability. Therefore, they estimated the mass loss rate based on statistical results from previous experiments. Mao et al. conducted experimental studies in both open environments and combustion chambers to investigate the combustion behavior of 18650 LIBs with NMC 532 cathodes. Their research focuses on the correlation between heat release, mass loss during TR, and the state of charge (SOC) of the LIBs [39]. Based on this and the estimation of the generated gas components from other experimental studies, they subsequently developed a lumped model to characterize the jet flow and fire dynamics of these LIBs under TR [23]. Similarly, Kong et al. developed a numerical model based on the assumption of gas generation to describe the venting dynamics and combustion behavior of 18650 LIBs with NCA cathodes. However, their model did not focus on solid particle ejection and mass loss [7]. Wang’s team addressed this gap with their study by developing a multi-scale model for TR that integrates the ejection of solid particles and mass loss [22]. They assumed that the mass loss from solid particle ejection is proportional to the total mass loss of the battery and conducted experiments in an open environment for venting flow analysis and in a sealed chamber to characterize the ejected particles.

The above literature review clearly indicates that the current trend in the field of TR modeling has shifted from detailing heat generation mechanisms to predicting TR behavior. While thermal behavior, such as temperature change and heat release, remain primary focuses, there is growing attention to venting and combustion behavior, including gas generation and release, particle ejection, venting flow dynamics, and mass loss rate in recent years. To characterize and quantify the venting behavior of cylindrical LIBs and sodium-ion batteries (SIBs) during TR events, Fedoryshyna et al. developed a custom test bench equipped with shear force sensors and a high-speed camera to capture detailed data on venting dynamics, including gas velocity, mass flow rate, and temperature variations [33]. They tested 21700 and 18650 LIBs with NMC811/SiC chemistry and 18650 SIBs with ${\mathrm{NaMn}}_{1/3}{\mathrm{Fe}}_{1/3}{\mathrm{Ni}}_{1/3}{\mathrm{O}}_2$/HC chemistry, all at a 100 % SOC. TR was triggered through continuous heating, and measurements focused on gas ejection velocity, recoil force, venting gas temperature, and mass loss during TR. Gillich et al. also developed a mechanical measurement approach to analyze venting behavior during TR in 18650 LIBs [34]. Using a 3-axis force sensor and high-frequency data logging, the setup enables precise measurements of recoil force and gas flow rates during TR in both high-power and high-energy 18650 cells.

Among these, predicting particle ejection and mass loss rate poses the greatest difficulties due to the randomness of ejection and explosion, which also brings challenges to the reliability of the experiments for validation [40]. Therefore, the impact of the experiment setup on the TR and venting process has also become a study focus, including the triggering method of TR, the SOC, the cell format etc. Willstrand et al. investigated how test setups, SOC, and trigger methods influence the TR characteristics and gas release in large-format lithium-ion cells (157 Ah) [40]. They conducted 37 tests on large-format prismatic cells in both closed and open setups with six trigger methods, including heating, local heating, nail penetration, overcharge, and external fire, and with varied SOC levels. Their test results demonstrated a significant impact from various TR trigger methods and different SOC levels, particularly on gas production rate, mass loss, and peak temperatures. They also noted that the trigger methods did not notably influence the total gas generation or gas composition, and the heating rate had minimal effect on the TR onset temperature.

1.3. Contribution and structure

This study presents a novel modeling approach to simulate the TR and venting process in LIBs, focusing on the accurate prediction of temperature development, gas generation, and mass loss. The model is validated by TR experiments using a type of high-energy 18650 cell (NMC811/SiC) in an open environment, triggered with heaters providing constant heating power, similar to the setup used in most of the referenced studies involving venting process [7,19,23,39]. By incorporating the venting process and mass loss behavior, the study extends the understanding of the TR phenomenon, which aligns with the current research trends in this area, offering valuable insights and recommendations for future research to improve LIB safety and support the broader adoption of sustainable energy solutions.

The article is structured as follows: Section 2 presents the model approaches and experiment setup for validation. The simulation and experiment results will be demonstrated and discussed in section 3. A conclusion will be then given in section 4.

2. Methodology

In this study, a novel modeling approach is developed to simulate the TR and venting process in LIBs, focusing on the temperature development, gas generation, and mass loss. Experiments using a type of high-energy 18650 cell - LG^®-INR18650-MJ1 are conducted to parameterize and validate the model. An overview of the setup conditions and methodology of this study is shown in Figure 1, which includes the schematic illustration representing the two different heating methods adopted in our experiments, the number of the experiments conducted under different heating conditions, and an example to explain the onset temperature - the point after which the cell temperature increases uncontrollably, and the peak temperature. All measurements discribed in this study have been conducted at MAP Battery Tests GmbH facilities.

2.1. Model Approach

Given the complexity of thermal reactions and venting dynamics during TR, advanced modeling approaches were implemented in MATLAB^® to simulate the behavior of a LIB under TR. The flow chart of the modeling work done in this research is shown in Figure 2 with the main inputs and outputs. The model developed in this study comprises three distinct sub-models: (I) ’Thermal Runaway Sub-model’, (II) ’Internal Pressure Sub-model’, and (III) ’Venting Sub-model’. Each sub-model has a specific focus and their outputs serve as inputs for the others, establishing an interconnected framework. These outputs are continually compared with the experimental results to validate and refine the model, ensuring that the final outputs are satisfactorily aligned with observed data.

2.1.1. Thermal runaway sub-model

The flow chart of the TR sub-model is shown in Figure 3. The initiation of TR typically results from heat generated by triggered reactions at both anode, cathode and electrolyte, which can be differentiated into separate reactions for discussion. In this work, the LIB with the NMC cathode was chosen to be modeled. Based on the findings of Richard and Dahn and Kriston et al. [8,9,41], we assumed the following reactions for the used cell type in the anode: SEI decomposition (A1), SEI regeneration (A2), Li-electrolyte reaction (A3), and Anode-binder decomposition (A4), while in cathode and electrolyte, the reactions are: autocatalytic decomposition of NMC (C1), electrolyte reactions (C2), NMC binder decomposition (C3), and further decomposition of NMC with oxygen reactions (C4)

A1 - SEI decomposition 

When the temperature of a lithium-ion cell reaches around 90 °C, the SEI layer formed on the surface of the anode can decompose and release heat, which can contribute to triggering other reactions. The reaction has been modeled and validated by Richard et al [8] through the following Arrhenius equations Eqs. (1) and (2):

$${\displaystyle \frac{\mathrm{d}{x}_{\mathrm{SEId}}}{\mathrm{d}t}}=-A_{\mathrm{SEId}}{x}_{\mathrm{SEId}}^{0.5}exp\left(-{\displaystyle \frac{E_{\mathrm{SEId}}}{\mathrm{R}\mathrm{T}}}\right)$$

(1)

$${\dot{Q}}_{\mathrm{SEId}}=m_{\mathrm{an}}{h}_{\mathrm{SEId}}\left(-{\displaystyle \frac{\mathrm{d}{x}_{\mathrm{SEId}}}{\mathrm{d}t}}\right)$$

(2)

Equation (1) shows the rate of decomposition of the metastable part of the SEI layer as a function of time $\mathrm{d}{x}_{\mathrm{SEId}}/\mathrm{d}t$, which depends on the frequency factor $A_{\mathrm{SEId}}$, the activation energy $E_{\mathrm{SEId}}$, and the temperature T. Equation. (2) relates the amount of heat released during the SEI decomposition $Q_{\mathrm{SEId}}$ to the mass of anode material $m_{\mathrm{an}}$, the heat of the reaction $h_{\mathrm{SEId}}$, and the rate of decomposition $\mathrm{d}{x}_{\mathrm{f}}/\mathrm{d}t$, which is negative due to the decreasing amount of metastable SEI. Before decomposition, the initial mole fraction of the metastable part of the SEI layer is $x_{\mathrm{f}0}=0.1$[8]. The gas constant $\mathrm{R}$ was set to 8.314 ${\mathrm{Jmol}}^{-1}{\mathrm{K}}^{-1}$ in the calculations.

A2 - SEI regeneration 

The SEI regeneration process occurs after the decomposition of the SEI layer with a peak at around 200 °C. During this process, intercalated lithium reacts with the electrolyte, and deposits on the degenerated SEI film, resulting in the regeneration of a new film, which was modeled using Arrhenius-type equations as follows [8]:

$${\displaystyle \frac{\mathrm{d}{x}_{\mathrm{SEIr}}}{\mathrm{d}t}}=-A_{\mathrm{SEIr}}{x}_{\mathrm{SEIr}}exp\left(-{\displaystyle \frac{z_{\mathrm{i}}}{z_0}}\right)exp\left(-{\displaystyle \frac{E_{\mathrm{SEIr}}}{\mathrm{R}\mathrm{T}}}\right)$$

(3)

$${\displaystyle \frac{\mathrm{d}{z}_{\mathrm{i}}}{\mathrm{d}t}}=A_{\mathrm{SEIr}}{x}_{\mathrm{SEIr}}exp\left(-{\displaystyle \frac{z_{\mathrm{i}}}{z_0}}\right)exp\left(-{\displaystyle \frac{E_{\mathrm{SEIr}}}{\mathrm{R}\mathrm{T}}}\right)$$

(4)

$${\dot{Q}}_{\mathrm{SEIr}}=m_{\mathrm{an}}{h}_{\mathrm{SEIr}}\left(-{\displaystyle \frac{\mathrm{d}{x}_{\mathrm{SEIr}}}{\mathrm{d}t}}\right)$$

(5)

Here the rate of change of the intercalated lithium amount $x_{\mathrm{SEIr}}$ is related with the amount of the lithium in the SEI layer per unit surface area $z_{\mathrm{i}}$, which can be used to present the thickness of the SEI diffusion layer. Given that the growth of the thickness of the SEI layer impacts the transport of lithium through this layer and consequently affects the rate of consumption of intercalated Li, a factor of $\exp \left(-z_{\mathrm{i}}/z_0\right)$ is incorporated. This factor serves to decrease the consumption rate of intercalated Li as the SEI layer expands. The variable $z_0$ represents the initial amount of Li in the SEI per unit surface area, with a value of 0.15 ${\mathrm{m}}^{-2}$. The heat generated during the SEI regeneration process is modeled by Eq. (5) using the heat of the reaction $h_{\mathrm{SEIr}}$ [8].

A3 - Li-electrolyte reaction 

As the TR goes further, the Li-electrolyte reaction takes place at elevated temperatures, even before the completion of A2. During this reaction, the regenerated SEI decomposes, and lithium directly reacts with the electrolyte. This process is faster and more energetic than A2, consuming the intercalated lithium in the anode and generating hydrocarbons. Based on the specific electrolytes used, these subsequent reactions may happen [41]:

$$2\mathrm{Li}+{\mathrm{C}}_3{\mathrm{H}}_4{\mathrm{O}}_3\left(\mathrm{EC}\right)\to \mathrm{L}{\mathrm{i}}_2\mathrm{C}{\mathrm{O}}_3+{\mathrm{C}}_2{\mathrm{H}}_4$$

(I)

$$2\mathrm{Li}+{\mathrm{C}}_3{\mathrm{H}}_6{\mathrm{O}}_3\left(\mathrm{DMC}\right)\to \mathrm{L}{\mathrm{i}}_2\mathrm{C}{\mathrm{O}}_3+{\mathrm{C}}_2{\mathrm{H}}_6$$

(II)

$$2\mathrm{Li}+{\mathrm{C}}_5{\mathrm{H}}_{10}{\mathrm{O}}_3\left(\mathrm{DEC}\right)\to \mathrm{L}{\mathrm{i}}_2\mathrm{C}{\mathrm{O}}_3+{\mathrm{C}}_4{\mathrm{H}}_{10}$$

(III)

As A3 begins before A2 is finished, both reactions consume intercalated lithium simultaneously, indicating that A2 and A3 should be coupled according to Kriston et al. [41], whose experiment has shown an abrupt high heat flow during the process of A2, fitting this coupled model represented by the following equation, where frequency factor $A_{\mathrm{LiE}}$ and activation energy $E_{\mathrm{LiE}}$ are used:

$${\displaystyle \frac{\mathrm{d}{x}_{\mathrm{SEIr}}}{\mathrm{d}t}}=-A_{\mathrm{SEIr}}{x}_{\mathrm{SEIr}}exp\left(-{\displaystyle \frac{z_{\mathrm{i}}}{z_0}}\right)exp\left(-{\displaystyle \frac{E_{\mathrm{SEIr}}}{\mathrm{R}\mathrm{T}}}\right)-A_{\mathrm{LiE}}{x}_{\mathrm{SEIr}}exp\left(-{\displaystyle \frac{E_{\mathrm{LiE}}}{\mathrm{R}\mathrm{T}}}\right)$$

(6)

The generated heat can still be calculated through Eq. (5), for the value of $h_{\mathrm{SEIr}}$ used in this work has already included the A3 reaction.

A4 - Anode-Binder decomposition 

Similarly, an Arrhenius-type single component decomposition is proposed as an approximation to represent anode-binder decomposition, with $x_{\mathrm{AnB}}$ being the dimensionless concentration of binder and other components from the previous decomposition process, and $A_{\mathrm{AnB}},E_{\mathrm{AnB}},h_{\mathrm{AnB}}$ as the thermal triplets for calculations [41]:

$${\displaystyle \frac{\mathrm{d}{x}_{\mathrm{AnB}}}{\mathrm{d}t}}=-A_{\mathrm{AnB}}{x}_{\mathrm{AnB}}exp\left(-{\displaystyle \frac{E_{\mathrm{AnB}}}{\mathrm{R}\mathrm{T}}}\right)$$

(7)

$${\dot{Q}}_{\mathrm{AnB}}=m_{\mathrm{an}}{h}_{\mathrm{AnB}}\left(-{\displaystyle \frac{\mathrm{d}{x}_{\mathrm{AnB}}}{\mathrm{d}t}}\right)$$

(8)

C1 - Autocatalytic decomposition of NMC 

During the process of TR, the NMC material undergoes a self-accelerating decomposition process at over 200 °C [12,42]. It is typical for all models to consider only one decomposition reaction and describe it through an Arrhenius-type equation [41,43]:

$${\displaystyle \frac{\mathrm{d}{x}_{\mathrm{NMCd}}}{\mathrm{d}t}}=-A_{\mathrm{NMCd}}{x}_{\mathrm{NMCd}}exp\left(-{\displaystyle \frac{E_{\mathrm{NMCd}}}{\mathrm{R}\mathrm{T}}}\right)$$

(9)

$${\dot{Q}}_{\mathrm{NMCd}}=m_{\mathrm{ca}}{h}_{\mathrm{NMCd}}\left(-{\displaystyle \frac{d{x}_{\mathrm{NMCd}}}{dt}}\right)$$

(10)

C2 - Electrolyte reactions 

As the temperature increases, the electrolyte evaporates at low temperatures initially as an endothermic reaction. Subsequently, the remaining electrolyte undergoes decomposition. Under specific temperatures and mixing conditions with the released oxygen, there is also a possibility of a combustion of the electrolyte.

The generation of oxygen can be modeled using an Arrhenius-type equation summarized by Kristion et al. [41]:

$${\displaystyle \frac{\mathrm{d}{x}_{\mathrm{O}2}}{\mathrm{d}t}}=-y_{\mathrm{O}2}{A}_{\mathrm{NMCd}}{x}_{\mathrm{NMCd}}exp\left(-{\displaystyle \frac{E_{\mathrm{NMCd}}}{\mathrm{R}\mathrm{T}}}\right),$$

(11)

which is highly related to the decomposition of NMC in C1. The factor $y_{\mathrm{O}2}$ is the stoichiometric amount of oxygen released with a value range from 0 to 1, which is chosen to be 0.1 in this work according to Kriston et al. [41]

Based on that, the consumption of the electrolyte through decomposition, evaporation, and combustion is coupled and modeled with their respective frequency factors A and activation energies E through the Eqs. (12) and (13), while the generated heat is calculated with the heat of each reactions h and the mass of electrolyte $m_{\mathrm{El}}$ using Eq. (14):

$${\displaystyle \frac{\mathrm{d}{x}_{\mathrm{EL}}}{\mathrm{d}t}}={\displaystyle \frac{\mathrm{d}{x}_{\mathrm{EL},\mathrm{eva}}}{\mathrm{d}t}}+{\displaystyle \frac{\mathrm{d}{x}_{\mathrm{EL},\mathrm{dec}}}{\mathrm{d}t}}+{\displaystyle \frac{\mathrm{d}{x}_{\mathrm{EL},\mathrm{com}}}{\mathrm{d}t}}$$

(12)

$${\displaystyle \frac{\mathrm{d}{x}_{\mathrm{EL}}}{\mathrm{d}t}}=-A_{\mathrm{eva}}{x}_{\mathrm{El}}exp\left(-{\displaystyle \frac{E_{\mathrm{eva}}}{\mathrm{R}\mathrm{T}}}\right)-A_{\mathrm{dec}}{x}_{\mathrm{El}}exp\left(-{\displaystyle \frac{E_{\mathrm{dec}}}{\mathrm{R}\mathrm{T}}}\right)-A_{\mathrm{com}}{x}_{\mathrm{El}}{x}_{\mathrm{O}2}exp\left(-{\displaystyle \frac{E_{\mathrm{com}}}{\mathrm{R}\mathrm{T}}}\right)$$

(13)

$$\begin{array}{cc}\hfill {\dot{Q}}_{\mathrm{El}}=& m_{\mathrm{El}}·(-A_{\mathrm{eva}}{x}_{\mathrm{El}}exp\left(-{\displaystyle \frac{E_{\mathrm{eva}}}{\mathrm{R}\mathrm{T}}}\right)·h_{\mathrm{eva}}\hfill \\ \hfill & -A_{\mathrm{dec}}{x}_{\mathrm{El}}exp\left(-{\displaystyle \frac{E_{\mathrm{dec}}}{\mathrm{R}\mathrm{T}}}\right)·h_{\mathrm{dec}}-A_{\mathrm{com}}{x}_{\mathrm{El}}{x}_{\mathrm{O}2}exp\left(-{\displaystyle \frac{E_{\mathrm{com}}}{\mathrm{R}\mathrm{T}}}\right)·h_{\mathrm{com}})\hfill \end{array}$$

(14)

Among these reactions, electrolyte combustion is influenced not only by oxygen release but also by the mixing ratio of the electrolyte and oxygen, as well as the corresponding ignition point. Due to this complexity, not all studies have accounted for electrolyte combustion [3]. In light of these challenges, electrolyte combustion was not considered in our modeling work, allowing for more accurate results.

C3 - NMC binder decomposition 

Based on the work of Kriston et al. [41], the NMC cathode binder material also decomposes during the TR process in the temperature range of 400 – 500 °C, which can be presented through the equations with the thermal triplets for this reaction $A_{\mathrm{CaB}},E_{\mathrm{CaB}},h_{\mathrm{CaB}}$:

$${\displaystyle \frac{\mathrm{d}{x}_{\mathrm{CaB}}}{\mathrm{d}t}}=-A_{\mathrm{CaB}}{x}_{\mathrm{CaB}}exp\left(-{\displaystyle \frac{E_{\mathrm{CaB}}}{\mathrm{R}\mathrm{T}}}\right)$$

(15)

$${\dot{Q}}_{\mathrm{CaB}}=m_{\mathrm{ca}}{h}_{\mathrm{CaB}}\left(-{\displaystyle \frac{\mathrm{d}{x}_{\mathrm{CaB}}}{\mathrm{d}t}}\right)$$

(16)

C4 - Further decomposition of NMC with oxygen reactions 

In the temperature range 500 °C – 600 °C, the released oxygen burns the carbon or other organic compounds, triggering a further decomposition of the NMC material [41]. The mechanism is similar to the C1 reaction and can also be modeled similarly:

$${\displaystyle \frac{\mathrm{d}{x}_{\mathrm{NMCf}}}{\mathrm{d}t}}=-A_{\mathrm{NMCf}}{x}_{\mathrm{NMCf}}exp\left(-{\displaystyle \frac{E_{\mathrm{NMCf}}}{\mathrm{R}\mathrm{T}}}\right)$$

(17)

$${\dot{Q}}_{\mathrm{NMCf}}=m_{\mathrm{ca}}{h}_{\mathrm{NMCf}}\left(-{\displaystyle \frac{\mathrm{d}{x}_{\mathrm{NMCf}}}{\mathrm{d}t}}\right)$$

(18)

Although oxygen release occurs during this process, the associated heat generation is incorporated into the terms $h_{\mathrm{CaB}}$ and $h_{\mathrm{NMCf}}$ in Eqs. (16) and (18). Therefore, it is not calculated separately in order to reduce computational time and maintain the stability of the model.

Temperature rise 

In our modeling approach, the cells are considered both thermally and chemically lumped, meaning that the temperature within a single cell is assumed to be uniform. To align with our experimental setup, the TR of the cell is triggered by an external heater with a fixed heating power until the cell reaches a specified temperature $T_{\mathrm{set}}$. Consequently, the temperature rising rate of the cell from the external heater can be expressed as:

$${\left.{\displaystyle \frac{\mathrm{d}\mathrm{T}}{\mathrm{d}\mathrm{t}}}\right|}_{\mathrm{cell}}={\displaystyle \frac{{\dot{Q}}_{\mathrm{heater}}}{c_{\mathrm{p},\mathrm{cell}}·m_{\mathrm{cell}}}}$$

(19)

The heat generated from the TR reactions is another factor contributing to the temperature rise. Since the heat generation of each reaction has already been discussed above, the heat flow from thermal runaway can be determined by summing the heat contributions from each reaction i from A1 to A4, and C1 to C4:

$${\dot{Q}}_{\mathrm{TR}}=\sum {\dot{Q}}_i={\dot{Q}}_{\mathrm{SEIr}}+{\dot{Q}}_{\mathrm{SEId}}+{\dot{Q}}_{\mathrm{AnB}}+{\dot{Q}}_{\mathrm{NMCd}}+{\dot{Q}}_{\mathrm{El}}+{\dot{Q}}_{\mathrm{CaB}}+{\dot{Q}}_{\mathrm{NMCf}}$$

(20)

Therefore, considering both the contribution of the external heater and internal reactions, the temperature rise of the cell can be concluded as:

$${\displaystyle \frac{\mathrm{d}\mathrm{T}}{\mathrm{d}\mathrm{t}}}=\left\{\begin{array}{cc}{\displaystyle \frac{{\dot{Q}}_{\mathrm{heater}}+{\dot{Q}}_{\mathrm{TR}}}{c_{\mathrm{p},\mathrm{cell}}·m_{\mathrm{cell}}}},\hfill & T\le T_{\mathrm{set}}\hfill \\ {\displaystyle \frac{{\dot{Q}}_{\mathrm{TR}}}{c_{\mathrm{p},\mathrm{cell}}·m_{\mathrm{cell}}}},\hfill & T>T_{\mathrm{set}}\hfill \end{array}\right.$$

(21)

2.1.2. Internal pressure sub-model

A large volume of gas production is also observed during the TR process, originating from the electrolyte evaporation [19] and the TR reactions described in the thermal runaway sub-model [20,41]. This results in a sharp increase in the internal pressure of the cell. In our work, the modeled cell is an 18650 type, equipped with a safety valve. When the internal pressure reaches a critical threshold, the safety valve opens, releasing a mixture of electrolyte vapors, reactive gases, and solid particles, thereby reducing the internal pressure [22,44]. This process is modeled with an internal pressure sub-model, with its workflow illustrated in Figure 4:

Gas generation 

As the temperature increases during TR with ongoing exothermic reactions, a significant amount of flammable or toxic gases can be generated as products of the TR reactions or side reactions, contributing to the increase of the internal pressure of the cell. This mechanism has been extensively studied through experiments and quantification work by various researchers [20,36,37,38]. Generally, the generation of these 6 types of gases are focused and studied: ${\mathrm{CO}}_2,\mathrm{CO},{\mathrm{H}}_2,{\mathrm{CH}}_4,{\mathrm{C}}_2{\mathrm{H}}_4,{\mathrm{C}}_2{\mathrm{H}}_6$[45]. The reaction mechanism within the cell is complex and filled with uncertainty. Consequently, it is reasonable to simplify the modeling by assuming that the generation of each gas during a particular reaction stage corresponds to a single primary chemical reaction. This approach facilitates a clearer modeling of the relationship between gas generation and the reaction extent [22]. The ${\mathrm{CO}}_2$ comes mainly from the decomposition of the SEI layer (A1), the reaction between the exposed lithium and electrolyte (A3), and the decomposition of electrolyte (C2) [20,41]. Hydrogen ${\mathrm{H}}_2$ is produced primarily from the decomposition of the binder material, specifically polyvinylidene fluoride (PVDF). [20]. Besides ${\mathrm{CO}}_2$, the generation of $\mathrm{CO},{\mathrm{CH}}_4,{\mathrm{C}}_2{\mathrm{H}}_4,$ and ${\mathrm{C}}_2{\mathrm{H}}_6$ can also be observed from the Li-electrolyte reaction (A3) [20,41]. The above-mentioned reaction mechanisms and their associated chemical equations are summarized in Table 2.

With the gas generation mechanisms established, the matrix shown ${\omega }_{a,b}$ in Eq. (22) can be employed to quantitatively assess the contributions of each TR reaction to the production of specific gas components:

$$\begin{array}{cc}\hfill & {\left(\begin{array}{c}\Delta x_{\mathrm{A}1}\\ \Delta x_{\mathrm{A}3}\\ \Delta x_{\mathrm{A}4}\\ \Delta x_{\mathrm{C}2}\\ \Delta x_{\mathrm{C}3}\end{array}\right)}^T\left(\begin{array}{cccccc}{\omega }_{{\mathrm{CO}}_2,1}& 0& 0& 0& 0& 0\\ {\omega }_{{\mathrm{CO}}_2,2}& {\omega }_{\mathrm{CO}}& 0& {\omega }_{{\mathrm{CH}}_4}& {\omega }_{{\mathrm{C}}_2{\mathrm{H}}_4}& {\omega }_{{\mathrm{C}}_2{\mathrm{H}}_6}\\ 0& 0& {\omega }_{{\mathrm{H}}_2,1}& 0& 0& 0\\ {\omega }_{{\mathrm{CO}}_2,3}& 0& 0& 0& 0& 0\\ 0& 0& {\omega }_{{\mathrm{H}}_2,2}& 0& 0& 0\end{array}\right)\hfill \\ \hfill & =(\Delta n_{{\mathrm{CO}}_2},\Delta n_{\mathrm{CO}},\Delta n_{{\mathrm{H}}_2},\Delta n_{{\mathrm{CH}}_4},\Delta n_{{\mathrm{C}}_2{\mathrm{H}}_4},\Delta n_{{\mathrm{C}}_2{\mathrm{H}}_6})\hfill \end{array}$$

(22)

In this matrix, rows (1) - (5) correspond to the TR reactions discussed, while columns (1) - (6) represent the gas components ${\mathrm{CO}}_2,\mathrm{CO},{\mathrm{H}}_2,{\mathrm{CH}}_4,{\mathrm{C}}_2{\mathrm{H}}_4,$ and ${\mathrm{C}}_2{\mathrm{H}}_6$. The amount of gas b produced in reaction a is quantified using the coefficients ${\omega }_{ab}$. Therefore, the amount of gas b produced over a specific time period $\Delta n_b$ can be calculated using the formula provided in Eq. (23):

$$\Delta n_b=\sum _{a=1}^5{\omega }_{ab}\Delta x_a$$

(23)

Golubklov et al., Shi et al., Shen et al. and Willstrand et al. have all conducted experiments on NMC LIBs [20,37,38,40]. The total amount of the gases produced in their studies correlates well with the capacity of the LIBs, as their results with 100 % SOC LIBs are detailed in Table 3:

As shown in Table 3, despite variations in the shape and specific elemental ratios within the cathode material, LIBs with the same cathode type tend to have similar ratios of total gas production to battery capacity. For example, NMC LIBs with 100 % SOC usually produce approximately 0.1 $\mathrm{mol}/\mathrm{Ah}$, whereas LFP LIBs yield around 0.025 $\mathrm{mol}/\mathrm{Ah}$. Furthermore, given that the gas generation mechanism is directly related to the TR reactions, which is theoretically the same for LIBs with identical materials, the amount fraction of each type of generated gas should also remain consistent. This consistency has been confirmed through experiments on LIBs with the same cathode material and identical elemental ratios. For example, the component of the generated gas from the experiments with NMC 9 0.5 0.5 LIBs is summarized in Table 4:

From these listed experiments conducted on NMC9 0.5 0.5 LIBs, a significant correlation emerges in the composition of the main gases generated, with $\mathrm{CO}$ being the most abundant, comprising approximately 40 % of the total. Additionally, ${\mathrm{CO}}_2$ and ${\mathrm{H}}_2$ are also prevalent, each accounting for over 20 % of the gases produced. These findings suggest that it is feasible to use data from previous experiments to predict the amount of gases generated by TR in LIBs with the same cathode type.

Internal pressure evolution 

With the understanding of gas generation and the calculation of gas generation amount, the internal pressure development can be modeled. Venting is initiated when the internal pressure of the cell surpasses a critical threshold, typically around 1900 $\mathrm{kPa}$ for 18650 LIBs [7]. Before that, the internal pressure $P_{\mathrm{cell}}$ is attributed to three principal factors: the pressure from electrolyte evaporation $P_{\mathrm{eva}}$, the pressure from gas generation $P_{\mathrm{gas}}$, and the original existed pressure $P_{\mathrm{amb}}$, which is represented in Eq. (24) [7]:

$$P_{\mathrm{cell}}=P_{\mathrm{eva}}+P_{\mathrm{gas}}+P_{\mathrm{amb}}$$

(24)

$P_{\mathrm{eva}}$ is quantified using Antoine’s equation Eq. (25) [46], which reflects how the pressure from evaporated electrolyte escalates as the cell’s temperature $T_{\mathrm{cell}}$ rises:

$$P_{\mathrm{eva}}={10}^{9.4338-{\displaystyle \frac{1413}{-44.25+T_{\mathrm{cell}}}}}\left[\begin{array}{c}\mathrm{Pa}\end{array}\right]$$

(25)

Simultaneously, pressure from gas generation within the cell $P_{\mathrm{gas}}$ is determined using the ideal gas law Eq. (26) [7]:

$$P_{\mathrm{gas}}={\displaystyle \frac{n_{\mathrm{gas}}R{T}_{\mathrm{cell}}}{V_{\mathrm{h}}}}$$

(26)

Here $n_{\mathrm{gas}}$ represents the moles of gas produced, $\mathrm{R}$ is the universal gas constant, $T_{\mathrm{cell}}$ is the temperature, and $V_{\mathrm{h}}$ is the volume of the headspace within the battery occupied by the vapor and generated gas.

After the internal pressure reaches the critical pressure, the safety valve is triggered to open, which releases the accumulated gases and breaks the vapor-liquid equilibrium condition inside the cell. Based on the work of Kong et al. [7], a set of ordinary differential equations is formulated to simulate the dynamic change of internal pressure over time.

As concluded with Eq. (27), the changing rate of the internal pressure $d{P}_{\mathrm{cell}}/dt$ is attributed to the changing rates of the partial pressure from electrolyte vapor $d{P}_{\mathrm{eva}}/dt$ and from generated gases $d{P}_{\mathrm{gas}}/dt$:

$${\displaystyle \frac{d{P}_{\mathrm{cell}}}{dt}}={\displaystyle \frac{d{P}_{\mathrm{eva}}}{dt}}+{\displaystyle \frac{d{P}_{\mathrm{gas}}}{dt}}$$

(27)

Recognizing that the mass of vapor changes due to the expulsion of the gas mixture and that the vapor temperature varies along with the cell’s temperature, the rate of change in the partial pressure attributable to electrolyte vapor is determined. This is accomplished through a differential representation of the ideal gas law, as expressed in Eq. (28) with the mole mass of the electrolyte vapor represented by $M_{\mathrm{eva}}$:

$${\displaystyle \frac{d{P}_{\mathrm{eva}}}{dt}}={\displaystyle \frac{R{T}_{\mathrm{cell}}}{M_{\mathrm{eva}}{V}_{\mathrm{h}}}}{\displaystyle \frac{d{m}_{\mathrm{eva}}}{dt}}+{\displaystyle \frac{m_{\mathrm{eva}}R}{M_{\mathrm{eva}}{V}_{\mathrm{h}}}}{\displaystyle \frac{d{T}_{\mathrm{cell}}}{dt}}$$

(28)

Under the assumption that the composition of the venting mixture remains proportional, the mass flow rate of vapor $d{m}_{\mathrm{eva}}/dt$ can be estimated using the system’s overall mass loss rate, $d{m}_{\mathrm{loss}}/dt$ and the mass fraction of electrolyte vapor in the venting mixture, ${\phi }_{\mathrm{eva}}$, as represented in Eq. (29):

$${\displaystyle \frac{d{m}_{\mathrm{eva}}}{dt}}=-{\phi }_{\mathrm{eva}}{\displaystyle \frac{d{m}_{\mathrm{loss}}}{dt}}$$

(29)

Similarly, the rate of change in the partial pressure attributable to the generated gases is also determined by applying a differential form of the ideal gas law Eq. (30):

$${\displaystyle \frac{d{P}_{\mathrm{gas}}}{dt}}={\displaystyle \frac{R{T}_{\mathrm{cell}}}{V_{\mathrm{h}}}}{\displaystyle \frac{d{n}_{\mathrm{gas}}}{dt}}+{\displaystyle \frac{n_{\mathrm{gas}}R}{V_{\mathrm{h}}}}{\displaystyle \frac{d{T}_{\mathrm{cell}}}{dt}}$$

(30)

Unlike the mass flow rate of electrolyte vapor, the rate of change in the amount of gases generated from TR reactions is influenced not only by venting, but also by the continuous production of gases from ongoing TR reactions. This dual-source contribution to the gas generation rate is captured in Eq. (31):

$${\displaystyle \frac{d{n}_{\mathrm{gas}}}{dt}}={\displaystyle \frac{d{n}_{\mathrm{gas}}}{dt}}{|}_{\mathrm{generation}}-{\displaystyle \frac{{\phi }_{\mathrm{gas}}}{M_{\mathrm{gas}}}}{\displaystyle \frac{d{m}_{\mathrm{loss}}}{dt}}$$

(31)

Here the gas generation amount $d{n}_{\mathrm{gas}}/dt{|}_{\mathrm{generation}}$ is calculated following the methodology outlined in former section of gas generation, and the mole mass of generated gases, $M_{\mathrm{gas}}$ is calculated with Eq. (32) as follow. Meanwhile, the mass fraction of the gases within the venting mixture is denoted as ${\phi }_{\mathrm{gas}}$.

$$M_{\mathrm{gas}}={\displaystyle \frac{{\sum }_{b=1}^6{M}_{\mathrm{gas},\mathrm{b}}·n_{\mathrm{gas},\mathrm{b}}}{n_{\mathrm{gas},\mathrm{generation}}}}$$

(32)

2.1.3. Venting sub-model

With the internal pressure and gas generation already acquired, the venting sub-model is then developed, which focuses on the venting dynamics, including the characterization of mass and energy losses during the venting process, with the workflow shown in Figure 5.

Venting dynamics calculation 

Based on the work from Kong et al. [7], to calculate the venting dynamics within this sub-model, the flow condition — whether subsonic or choked - must first be determined. The flow is considered subsonic if the ambient pressure $P_{\mathrm{amb}}$ is greater than a function of the internal pressure $P_{\mathrm{cell}}$ factored by the heat capacity ratio $\gamma$, represented in Eq. (33). Otherwise, the flow is regarded as choked.

$$P_{\mathrm{amb}}>P_{\mathrm{cell}}{\left({\displaystyle \frac{2}{\gamma +1}}\right)}^{{\displaystyle \frac{\gamma }{\gamma -1}}}$$

(33)

Upon establishing the flow regime, the Mach number M, venting pressure $P_{\mathrm{vent}}$, venting temperature $T_{\mathrm{vent}}$, and venting velocity $v_{\mathrm{vent}}$ can be calculated using the isentropic nozzle flow equations and pseudo diameter approach [7,19,22,47]. If the flow is determined to be subsonic the following equations are used:

$$M=\sqrt{{\displaystyle \frac{2}{\gamma -1}}\left[{\left({\displaystyle \frac{P_{\mathrm{cell}}}{P_{\mathrm{amb}}}}\right)}^{{\displaystyle \frac{\gamma -1}{\gamma }}}-1\right]}$$

(34)

$$P_{\mathrm{vent}}=P_{\mathrm{amb}}$$

(35)

$$T_{\mathrm{vent}}=T_{\mathrm{cell}}{\left({\displaystyle \frac{P_{\mathrm{amb}}}{P_{\mathrm{cell}}}}\right)}^{{\displaystyle \frac{\gamma -1}{\gamma }}}$$

(36)

If the flow is determined to be choked, the Mach number is set to 1 and the pressure and venting temperature are determined by:

$$M=1$$

(37)

$$P_{\mathrm{vent}}=P_{\mathrm{cell}}{\left({\displaystyle \frac{2}{\gamma +1}}\right)}^{{\displaystyle \frac{\gamma }{\gamma -1}}}$$

(38)

$$T_{\mathrm{vent}}=T_{\mathrm{cell}}\left({\displaystyle \frac{2}{\gamma +1}}\right)$$

(39)

These parameters are fundamental to determining the density of the gaseous venting mixture ${\rho }_{\mathrm{vent}},\mathrm{gas}$ with the Eq. (40), where $M_{\mathrm{gas}},\mathrm{mix}$ is the molar mass of the gas mixture.

$${\rho }_{\mathrm{vent},\mathrm{gas}}={\displaystyle \frac{P_{\mathrm{vent}}{M}_{\mathrm{gas},\mathrm{mix}}}{R{T}_{\mathrm{vent}}}}$$

(40)

Mass loss calculation 

The material ejected from a battery includes a gaseous venting mixture and solid particles. With the determination of the density of this combined ejection material ${\rho }_{\mathrm{vent}}$, the venting speed $v_{\mathrm{vent}}$ is calculated with Eq. (41):

$$v_{\mathrm{vent}}=M\sqrt{{\displaystyle \frac{\gamma P_{\mathrm{vent}}}{{\rho }_{\mathrm{vent}}}}}$$

(41)

The study from Wang et al. [22] and Finegan et al. [48] unfolds the ejection process of the particles in four distinct stages. (Stage I): Initially, the battery is at rest, with the jelly roll and safety valve intact. (Stage II): As TR progresses, substantial gas release raises the cell temperature, leading to a pressure differential that fractures the electrode materials. Concurrently, binder decomposition at elevated temperatures causes electrode particle debonding. These factors contribute to inter-particle fracture and the formation of gas pockets that drive the fragmented electrode materials toward the headspace. (Stage III): As these materials are propelled by the gas flow, they move through newly formed channels until they are expelled into the ambient environment, particularly accelerated near the safety valve. (Stage IV): This expulsion process creates channels within the gas pockets, leading to the observable collapse of the spiral-wound layers.

These mechanisms indicate that the mass flow rate of LIBs during venting is hard to describe by an exact principle, particularly concerning the ejection of solid particles. Nevertheless, experimental evidence suggests a similarity in the total mass loss rate for identical types of LIBs under well-controlled conditions [38]. The fraction of mass loss attributed to solid particle ejection has been the focus of several studies: Coman et al. [19] performed numerous experiments under consistent conditions and utilized the average total mass loss rate for their predictive modeling. Shen et al. [37] concentrated on quantifying the ratio of the mass of ejected solid particles to that of the generated gases. Recognizing the difficulty in precise measurement and formulation of the mass flow rate of ejected solid particles, Wang et al. [22] postulated that these rates are proportional to the total mass loss, forming the basis for their model. Accordingly, in the context of this work, we adopt the assumption that the mass flow rate of solid particles is proportional to the total mass flow rate during venting. Under this assumption, the relationship between the mass flow rate of the ejected solid particles and that of the gas mixture can be described with the ratio k:

$${\displaystyle \frac{d{m}_{\mathrm{particles}}}{dt}}=k{\displaystyle \frac{d{m}_{\mathrm{vent},\mathrm{gas}}}{dt}}$$

(42)

With $d{m}_{\mathrm{loss}}/dt=d{m}_{\mathrm{vent},\mathrm{gas}}/dt+d{m}_{\mathrm{particles}}/dt$, the partial mass flow rates are determined as:

$${\displaystyle \frac{d{m}_{\mathrm{particles}}}{dt}}={\displaystyle \frac{k}{k+1}}{\displaystyle \frac{d{m}_{\mathrm{loss}}}{dt}}$$

(43)

$${\displaystyle \frac{d{m}_{\mathrm{vent},\mathrm{gas}}}{dt}}={\displaystyle \frac{1}{k+1}}{\displaystyle \frac{d{m}_{\mathrm{loss}}}{dt}}$$

(44)

The total density of the venting mixture ${\rho }_{\mathrm{vent}}$ is then estimated with k as in Eq. (45):

$${\rho }_{\mathrm{vent}}={\displaystyle \frac{k+1}{\frac{k}{{\rho }_{\mathrm{cell}}}+\frac{1}{{\rho }_{\mathrm{vent},\mathrm{gas}}}}}$$

(45)

Therefore, with the density of venting mixture, the discharge coefficient $C_{\mathrm{d}}$ considering the non-ideal effect, and the vent area of the cell $A_v$, the total mass flow rate $d{m}_{\mathrm{loss}}/dt$ can be determined with Eq. (46) [7], while the partial mass flow rate can then be calculated combining the Eqs. (43) and (44).

$${\displaystyle \frac{d{m}_{\mathrm{loss}}}{dt}}=C_{\mathrm{d}}{A}_v{\rho }_{\mathrm{vent}}{v}_{\mathrm{vent}}$$

(46)

Energy loss calculation 

In the assessment of energy dissipation during the venting process associated with mass loss in LIBs, the sub-model incorporates several equations to describe the energy changes due to the ejection of various materials based on the work of Coman et al. [19].

The energy required for electrolyte evaporation is absorbed from the internal energy of the battery. This heat loss is quantified by ${\dot{Q}}_{\mathrm{eva}}$, as represented in Eq. (47) with the mass rate of evaporated electrolyte $-d{m}_{\mathrm{eva}}/dt$, and its latent heat of vaporization $L_{\mathrm{vap}}$.

$${\dot{Q}}_{\mathrm{eva}}=-{\displaystyle \frac{d{m}_{\mathrm{eva}}}{dt}}·L_{\mathrm{vap}}$$

(47)

Similarly, the rate of heat loss from gas venting, ${\dot{Q}}_{\mathrm{vent},\mathrm{gas}}$, is calculated with Eq. (48) by multiplying the mass rate of the venting gas $-d{m}_{\mathrm{vent},\mathrm{gas}}/dt$, the specific heat capacity of the gas $c_{p,\mathrm{gas}}$, and the temperature of the venting gas $T_{\mathrm{vent}}$, which reflects the heat carried away by the gases escaping from the battery.

$${\dot{Q}}_{\mathrm{vent},\mathrm{gas}}=-{\displaystyle \frac{d{m}_{\mathrm{vent},\mathrm{gas}}}{dt}}{c}_{p,\mathrm{gas}}{T}_{\mathrm{vent}}$$

(48)

For solid particle ejection, ${\dot{Q}}_{\mathrm{ej}}$ denotes the rate of heat loss, is calculated with Eq. (49), with the specific heat capacity estimated to be the specific heat capacity of the cell, $c_{p,\mathrm{cell}}$. This term signifies the heat associated with the ejected solid particles as they are expelled from the cell [22].

$${\dot{Q}}_{\mathrm{ej}}=-{\displaystyle \frac{d{m}_{\mathrm{particles}}}{dt}}{c}_{p,\mathrm{cell}}{T}_{\mathrm{vent}}$$

(49)

Collectively, the total rate of temperature change within the cell, $d{T}_{\mathrm{cell}}/dt$, is derived from the sum of the heat dissipation rates due to electrolyte evaporation, gas venting, and solid particle ejection with Eq. (50), providing an overarching view of how the cell’s temperature is influenced by these concurrent processes during a TR event.

$${\displaystyle \frac{d{T}_{\mathrm{cell}}}{dt}}={\displaystyle \frac{{\dot{Q}}_{\mathrm{eva}}+{\dot{Q}}_{\mathrm{vent},\mathrm{gas}}+{\dot{Q}}_{\mathrm{ej}}}{m_{\mathrm{cell}}{c}_{p,\mathrm{cell}}}}$$

(50)

2.2. Experiment Setup

In order to ensure that the developed model can reflect the complexity of real-world scenarios, experiments have been designed and conducted to validate the model in this study. The battery under test was the LG^®-INR18650-MJ1, featuring NMC811 chemistry and was charged to 100 % SOC for all tests. As for the trigger of the TR, multiple different constant heating powers - 20 W, 30 W, 50 W, 70 W and 90 W were applied to the battery via a 15 × 54 mm handmade NiCr wire heater fixed to the surface of the cell, either covering the bottom of the cell or vertically the full body of the cell, as shown in (a) and (b) in Figure 6. To monitor temperature changes with precision, four J-type thermocouples were positioned evenly across the battery surface, so the mean value from them could be taken to represent the cell temperature, as represented in (c) in Figure 6. Considering the high temperature and fast-changing rate during the TR process, mass losses were recorded using a high-resolution electronic balance only before and after the experiment, while a load cell SIEMENS^® SIWAREX WL260 SP-S AA was used during the experiment to capture the real-time mass loss resulting from the TR and venting process. The fixture system and load cell are represented in Figure 7.

A visualization system was also employed in our experiment, consisting of a high-speed camera Cyclone-65-70 manufactured by Optronis^®, which was capable of capturing 50 to 169 frames per second. This setup was enhanced by a light source arrangement: a laser emitting at 640 $\mathrm{nm}$ with a power of 10 $\mathrm{mW}$, or a spotlight depending on the effect, which was used to illuminate the battery and ensure clear visualization of the TR event as it unfolds, as shown in Figure 8.

2.3. Simulation Setup

The LIB sample employed in this study is the NMC811 INR18650MJ1 cell from LG^®, as chosen for the experiment, featuring a capacity of 3.5 Ah, a mass of around 46.5 g, a diameter of 18 mm, and a length of 65 mm. The cell’s electrolyte is a solvent comprised of ethylene carbonate (EC), dimethyl carbonate (DMC) and propylene carbonate (PC), while the anode is comprised of graphite containing minor silicon content, and the principle cathode constituent is ${\mathrm{Li}}_{\mathrm{x}}{\mathrm{Ni}}_{0.83}{\mathrm{Mn}}_{0.07}{\mathrm{Co}}_{0.1}{\mathrm{O}}_2$, as measured in work [49,50,51]. The main components of the sample used in this study align with the sample used in the experiments of Kriston et al. [41]. Therefore, some TR reaction kinetics and their modifications employed in this work have been primarily adapted from Kriston’s work.

The heat generated from the TR event originating from the anode, cathode, and electrolyte is calculated separately through the model. This calculation requires knowledge of the mass of the anode, cathode, and electrolyte, which could not be obtained from the manufacturer. Based on Golubkov’s in-lab measurements after disassembling and drying the cells [20], the mass of the anode is estimated to be 15 % of the cell mass, while the cathode accounts for 26 %. The mass percentage of the electrolyte, as discussed by Lebedeva et al. [52], can affect the performance of the cell and may vary. In this study, the mass of the electrolyte involved in reactions is estimated to be 10 % of the cell mass, with the mass ratio of the composition EC:DMC:PC is $6:3:1$ [49].

As for the trigger of the TR, the heating power and the cut-off temperature of the external heater are set to be the same with value of the experiment. A constant temperature boundary is applied to the outermost side of the cell, representing 12 °C air that transfers heat through convection with the cells, as the only heat transfer method with the surrounding environment considered. The temperature of the cell is considered to be homogenously distributed in this work, while the gases all fit the ideal gas law and are also homogenously distributed in the reaction volume and venting mixture.

The primary focus of parametrization work for TR sub-model lies in determining the so-called thermal triplets: frequency factor A, activation energy E, and heat of reaction h. These parameters have been investigated in numerous studies, employing commonly used methods such as accelerated rate calorimetry (ARC), differential scanning calorimetry (DSC), differential thermal analysis (DTA), or thermal gravimetric analysis (TGA) coupled with gas analysis [8,12,25,41,53]. In this work, Richard and Dahn’s data is employed for A1 and A2 reactions, while Kriston’s data is applied for A3, A4, and cathode reactions. The parameters are summarized in Table 5.

As for the prediction of gas generation, the generation amount parameters ${\omega }_{ab}$ used in matrix ${\omega }_{a,b}$ Eq. (22) are estimated based on the work from Shen et al. with the experiment result with the same LIB cathode material NMC811 [37]. With the approaches for TR reactions and the parameters considered in this work, the matrix ${\omega }_{a,b}$ is then modified to matrix Eq. (51):

$$\begin{array}{cc}& {\left(\begin{array}{c}\Delta x_{\mathrm{A}1}\\ \Delta x_{\mathrm{A}3}\\ \Delta x_{\mathrm{A}4}\\ \Delta x_{\mathrm{C}2}\\ \Delta x_{\mathrm{C}3}\end{array}\right)}^T\left(\begin{array}{cccccc}0.275& 0& 0& 0& 0& 0\\ 0.033& 0.080& 0& 0.190& 0.055& 0.002\\ 0& 0& 0.025& 0& 0& 0\\ 0.033& 0& 0& 0& 0& 0\\ 0& 0& 0.060& 0& 0& 0\end{array}\right)\hfill \\ & =(\Delta n_{{\mathrm{CO}}_2},\Delta n_{\mathrm{CO}},\Delta n_{{\mathrm{H}}_2},\Delta n_{{\mathrm{CH}}_4},\Delta n_{{\mathrm{C}}_2{\mathrm{H}}_4},\Delta n_{{\mathrm{C}}_2{\mathrm{H}}_6})\hfill \end{array}$$

(51)

The mass flow rate fraction ratio is obtained by fitting the simulation result with the experiment results, while the estimation of other parameters for the internal pressure sub-model and venting sub-model have taken reference from [7,54], where also 18650 type of LIB was studied. The principle parameters and the values are listed in the following Table 6.

3. Results and discussion

In this section, the results of the experiments and simulations with LG MJ1 batteries are presented and discussed. Outcomes from simulations are compared with the experimental results from our study and the existing literature, which evaluates the employed modeling approaches in terms of their precision and robustness in predicting the multifaceted behaviors of TR in LIBs.

3.1. Experiment results

3.1.1. Thermal runaway evolution

The entire TR process was recorded using the visualization system in the experiments. The phases of thermal runaway observed in an experiment with a 20 W vertically attached heater are presented in Figure 9, while the corresponding temperature-mass change curve is shown in Figure 10.

Here two distinct venting events were observed. The first venting event (Stage II in Figure 9) is noticeable at the early heating stage, occurring at approximately 120 °C. It is marked by the opening of the safety valve, which coincides with a minor mass loss, 14 % of the total battery mass, and a dip in temperature, reflecting the release of pressure and thermal energy from the cell. Following this venting episode, the battery enters a relatively stable phase (Stage III), where no further obvious venting behaviors are observed. However, the energy released by the venting is not able to stop the progression of TR, as the cell temperature continues to rise, while the generation of gases continues within the battery. Subsequently, the second venting event (Stage IV) happens just before the TR burst, at a stage where the battery already exhibits elevated internal pressure and temperature. This venting is considerably more intense, rapidly followed by an explosion (Stage V), where the TR process has reached a peak. The explosion ejected a large amount of particles, and the resulting impact caused the load cell to oscillate up and down, thereby recording highly fluctuating pressure values, making it impossible to read the data at that moment. Once the load cell stabilized, the recorded readings indicated a significant mass loss of 60 %. Following this phase, the TR process concludes swiftly, and the battery transitions into the cooling down stage (Stage VI), characterized by a rapid decrease in temperature.

The temperature at which a cell undergoes TR, known as the onset temperature, is a crucial parameter to identify, as it marks the point at which the cell’s temperature begins to rise uncontrollably. Understanding this temperature enables the development of safety strategies designed to prevent further damage and mitigate the risks associated with TR. Similarly, the peak temperature reached during a TR event is a vital parameter to record, as it indicates the maximum severity of the thermal event and its potential consequences. Knowing the peak temperature allows for the evaluation of safety margins, the testing of material thermal stability, and the optimization of thermal management systems to withstand extreme conditions. Therefore, both the onset and peak temperatures in each experiment in this study were recorded and summarized, as shown in Figure 11:

The part (a) of Figure 11 presents a box plot of the onset temperatures $T_{\mathrm{onset}}$ recorded from the experiments. The experiments using the full-body heating method exhibit similar median onset temperatures around 180 °C, with the onset temperatures generally being lower for experiments conducted at higher heating power (30 W). In contrast, the results from the bottom heating methods show significantly lower onset temperatures compared to the full-body heating methods, with the onset temperatures decreasing further as the heating power increases. These observations indicate that both localized heating at the bottom and the application of higher heating power accelerate the onset of TR, triggering TR at lower temperatures.

The peak temperature results, $T_{\max }$, are presented as a box plot in part (b) of Figure 11. The experiments using the full-body heating methods at 20 W and 30 W demonstrate higher peak temperatures, with median values around 500 °C and 450 °C, and show a wider range of temperature variation. In comparison, the bottom heating methods (50 W, 70 W, and 90 W) produce lower peak temperatures, with median values between 350 °C and 400 °C. This suggests that localized heating at the bottom produces less severe TR outcomes, characterized by lower maximum temperatures reached during TR events.

The results show that both the onset and peak temperatures are significantly affected by the heating method and power level, which aligns well with the work from Willstrand et al. [40], where the onset temperature resulted at around 150 °C, while the peak temperature varied from 450 °C to 900 °C. Full-body heating results in higher onset and peak temperatures, indicating a more intense and uniform thermal reaction throughout the battery. In contrast, bottom heating triggers TR at lower onset temperatures and results in lower peak temperatures, suggesting that TR events occur earlier but are less severe. This difference can be attributed to the uniform heat distribution provided by the full-body heating method, which takes longer to reach the critical temperature but causes the entire battery to undergo a more severe TR event simultaneously. Additionally, since the simulation uses a lumped model that approximates the heating method as equivalent to full-body heating, this section will also include a comparison and discussion of the simulation results with the experimental findings.

Besides temperature, the mass loss from a TR event is also a key parameter to quantify, as it indicates both the severity of the venting event and the extent of material degradation within the battery. This information helps in identifying the specific components most susceptible to degradation, allowing for the optimization of battery designs to reduce the likelihood of venting or to control the release of materials more safely. Additionally, mass loss data provides valuable input for refining predictive models of TR behavior by offering more detailed and accurate information. Consequently, the mass loss fractions observed in the experiments were carefully measured and summarized as a box plot in Figure 12.

The median mass loss fractions from all experiments, regardless of the heating method, are within the range of 70 % to 77 % regarding to the initial cell mass. The full-body heating method shows greater variability in mass loss, from approximately 60 % to 75 % with 20 W heating power, and from 65 % to 75 % with 30 W. This variability is likely due to the uniform heating of the entire cell, leading to varying venting dynamics. In contrast, the bottom heating method shows a more consistent mass loss fraction, with values ranging from 76 % to 77 % at 70 W and from 74 % to 76 % at 90 W, indicating that localized heating results in a more predictable pattern of material degradation and venting. This outcome is highly consistent with the findings reported by Fedoryshyna et al. and Gillich et al. [33,34], where the mass loss rate from TR tests with 18650 MJ1 cells ranged from 55 % to 85 %, with a median value of 75%. These relatively narrow mass loss fractions not only validate the reliability of our experimental results but also facilitate a more straightforward assumption of a specific mass loss rate in modeling work, which in turn reduces the complexity and unpredictability associated with venting behavior during TR events.

The standard deviation of the experiment results can be conclude in the following Table 7.

Here the standard deviations of the results from each heating method are similar. Although the standard deviation is smaller in experiments with higher heating power, we frequently observed that either a heating power exceeding 30 W or the bottom heating method could cause localized high temperatures. This would lead to a localized TR event, resulting in a significant amount of material being ejected through the newly formed breach. Consequently, the typical venting process we aimed to study could not be observed. Therefore, despite the relatively high standard deviation in the results, we selected the case with the 20 W full-body heating method for further simulation and study.

3.1.2. Venting analysis with load cell

The load cell used in our experiments acts as a pressure sensor, capturing data that reflects both the actual mass of the battery and the recoil force $F_{\mathrm{recoil}}$ generated during venting events. A schematic illustration of this force relationship is provided in part (a) of Figure 13. An accurate mass reading from the load cell is only possible when the system is in equilibrium. However, during venting events, particularly in the intense second venting phase, a substantial recoil force is generated due to the rapid ejection of materials, causing the recorded value to significantly exceed the gravitational force of the battery. To interpret the load cell data, we assume that the mass flow rate $\dot{m}$ remains constant during these brief venting intervals, enabling us to estimate both the battery’s mass and the corresponding recoil force. This estimation is crucial for determining the venting speed, which is calculated using the principle of conservation of momentum, as described by Eq. (52). The results of the estimation for the actual battery mass and venting speed are also represented in part (b) of Figure 13.

$$F_{\mathrm{vent}}=F_{\mathrm{recoil}}=\dot{m}{v}_{\mathrm{vent}}$$

(52)

From the data collected, the maximum venting speed was estimated to be 180 m/s. The load cell recorded a mass loss of 34.33 g (74 %), which closely matches the balance measurement of 34.31 g (74 %). This strong agreement between the load cell and balance measurements largely supports the assumptions made regarding the mass flow rate and recoil force during the venting process. However, the inherent randomness of TR and venting events causes significant oscillations in the recoil force, leading to occasional discrepancies between the load cell readings and the balance measurements, which is also shown in the result from the literature [33,34,40]. This suggests that while the load cell can be a useful tool, further refinement of the experimental setup is necessary to improve its reliability, which could be an important focus for future work.

3.2. Simulation results

3.2.1. Thermal runaway parameters

The simulation results offer a detailed perspective on the TR behavior of an LG MJ1 battery, highlighting the evolution of temperature, mass loss, and venting speed. Figure 14 illustrates the simulated progression of these parameters during TR in an LG MJ1 battery. The simulation was configured with a trigger method using a 20 W heater operating at 55 % efficiency, continuously heating the battery until it reached 200 °C, in order to replicate the experimental setup conditions. After calibrating the simulation with the experimental data, it was assumed that the mass flow rate of solid particles is 4.8 times that of the gas mixture.

As shown in Figure 14, the simulation identified two distinct venting phases, closely mirroring the experimental observations. From part (a), the first venting event is characterized by a smaller mass loss of 13 %, occurring at a lower temperature from 90 °C to 120 °C, which aligns with the early stage of TR shown in the experiment. The second, more significant venting phase is shown to start under conditions of higher temperature and pressure, followed by a burst of TR with a substantial mass loss of 61 %. From the venting speed result in part (b), this second venting also reached a maximum venting speed of approximately 180 m/s, consistent well with the estimated experimental data.

A detailed comparison between the simulation results and experimental data is presented in Figure 15, revealing a strong correlation, particularly in the temperature changes during the TR process. The peak temperatures recorded in the simulation and experiment are 478.3 °C and 473.0 °C, respectively, showing a difference of only 1 %, and resulted in the range of our other experiment results. Similarly, the onset temperatures are 183.8 °C in the simulation and 176.6 °C in the experiment, with a difference of 4 %, which also aligns greatly with our other experiment results. However, the experimental results exhibit time delays of approximately 250 seconds in key features, such as the first venting event, TR onset, and peak temperatures, compared to the simulation. These delays can be attributed to the simplifications inherent in the lumped model used in the simulation, such as the omission of internal heat transfer within the cell. Despite these discrepancies, the overall temporal progression from the first venting event to the onset of TR shows significant similarity, with 179 seconds in the experiment and 188 seconds in the simulation. This consistency demonstrates the predictive capability of the model, confirming its utility in simulating TR behavior.

The mass loss figures obtained from both the simulation and experimental methods further demonstrate a strong correlation, with both showing a mass loss of approximately 74 %. This value aligns closely with the mass loss percentages recorded in our other experiments conducted with fully charged LG MJ1 batteries. Such consistency across different tests strengthens the validity of the proportional assumption used in the simulation, specifically the assumption regarding the mass loss from solid particle ejection relative to the total mass loss.

These similarities between the experimental and simulated results validate the effectiveness of the simulation model. The model’s ability to replicate these critical aspects of the TR process indicates that it accurately reflects the dynamics occurring within the battery during TR events.

3.2.2. Gas generation analysis

As shown in part (a) of Figure 16, the gas generation result from the simulation indicates that ${\mathrm{H}}_2$, ${\mathrm{CO}}_2$, and $\mathrm{CO}$ are the primary gases produced, with different patterns in their timing and quantity of release represented in part (b), which fits to the gas generation matrix (22) implemented.

During the first venting phase, ${\mathrm{CO}}_2$ is predominantly generated. This early release aligns with the SEI-decomposition which starts at around 90 °C and produces ${\mathrm{CO}}_2$. It is also shown from the simulation result that while ${\mathrm{CO}}_2$ emerges first, the majority of other gases, including ${\mathrm{H}}_2$ and $\mathrm{CO}$, are released during the second venting event. This phase corresponds with the ongoing of TR reactions, especially the Li-electrolyte reaction, where the most gases are sourced from.

A comparative analysis of the gas generation results from this simulation and the experimental studies reported in the literature [37] is represented in Figure 17.

The experimental findings by Shen et al. [37] indicate higher levels of ${\mathrm{CO}}_2$ and CO production, along with a greater overall gas volume compared to the simulation outcomes from our study. Specifically, our simulations generated a total of 0.238 mol of gases, while the literature, assuming that gas generation is proportional to battery capacity (as discussed in the ’Internal pressure sub-model’ section), predicts a theoretical output of 0.350 mol. Although there is a noticeable discrepancy, the strong correlation between our simulation results and the experimental data supports the validity of the assumptions regarding the quantity of gas generated. Nevertheless, the noted differences also highlight potential limitations in the current modeling approach that may require further refinement to enhance the accuracy of the model:

• Simplified Reaction Modeling: There are just primary gas generation reactions considered in the modeling approach, with the omitting of the complex secondary reactions that occur at higher temperatures and can produce additional gases, particularly ${\mathrm{CO}}_2$ and CO. This simplification likely contributes to the lower total gas output in the simulation.
• Heating Power and TR Completion: The simulation possibly employed a lower heating power, insufficient to fully drive the battery through all the reactions in TR. This would naturally result in less gas generation as the TR reactions contributing to gas production would not fully develop or reach their peak activity.

These findings highlighted the need to adjust the simulation parameters by incorporating a broader range of chemical interactions and thermal effects. Such improvements could ensure that future simulations align more closely with empirical data, thereby providing more reliable predictions for safety assessments and battery design optimizations.

4. Conclusions

In this work, a comprehensive model was developed and validated to simulate TR and venting processes in 18650 NMC LIBs. The model successfully captures key elements of TR, such as temperature evolution, internal pressure changes, venting phases, and mass loss dynamics, and its validity was confirmed through experiments with LG^®-INR18650-MJ1 cells. The simulations have shown a strong correlation with the experimental results, particularly in replicating the double venting mechanisms, gas generation and mass loss characteristics observed during these events. The study also established the feasibility of modeling gas generation as proportional to cell capacity, and the mass loss rate from solid particle ejection as proportional to the total mass loss rate. These approaches not only simplify the modeling of mass loss behavior but also suggest a promising direction for future research efforts.

There are also areas where the model could be further refined by incorporating a more detailed treatment of internal heat transfer and expanding the range of chemical reactions considered. These enhancements would provide a more accurate representation of the TR process, particularly in predicting gas generation, and thereby improve the simulation of the venting process.

In conclusion, this work contributions to a deeper understanding of TR in LIBs by clarifying the gas generation and venting mechanisms involved. The findings presented here provide a foundation for future research focused on improving battery safety, which is essential for the widespread adoption of this key technology in sustainable energy applications.