Simulation Study of Thermal Runaway in Lithium-ion Battery Modules under Multi-Physics Coupling

Thermal runaway (TR) remains the paramount safety challenge hindering the widespread application of high-energy-density lithium-ion batteries. While significant research has focused on the thermal abuse behavior of single cells, a critical knowledge gap exists regarding the cascading failure mechanisms within series-connected modules, where complex thermal-electrical-chemical coupling dictates the propagation dynamics. This study systematically investigates the initiation and propagation characteristics of thermal runaway in a module comprising four 18650 NCA lithium-ion cells connected in series under three primary abuse conditions: internal short circuit (ISC), overcharge, and overdischarge. A three-dimensional thermal-electrical coupling model is developed within the ANSYS Fluent framework, incorporating a four-equation kinetic model to describe the exothermic chain reactions. By simplifying the cell geometry and accounting for anisotropic thermal conductivity, alongside the Bernardi heat generation model and the ampere-hour integration method, a robust simulation framework is established. The results elucidate the critical role of thermal coupling between cells in dictating the speed and severity of failure propagation. Under ISC, thermal runaway initiates at the fault location and spreads to adjacent cells primarily through axial heat conduction via busbars. Overcharge simulations reveal a 58% reduction in the time to trigger thermal runaway as the charging rate increases from 3C (10.2A) to 7C (23.8A), accompanied by accelerated reaction kinetics for the solid electrolyte interphase (SEI) decomposition and electrolyte reactions. Conversely, during overdischarge to critically low system voltages, the peak module temperature exhibits a substantial decrease, dropping by approximately 160°C when the voltage is lowered from 2V to 0V, highlighting a different failure progression dominated by voltage collapse and limited current flow. This research provides fundamental insights into the multi-field coupling behavior leading to catastrophic failure in lithium-ion battery packs, offering a theoretical basis for enhanced safety design and proactive thermal management strategies.

The relentless demand for higher energy density in lithium-ion battery systems, now exceeding 300 Wh/kg, is intrinsically linked to heightened thermal instability risks. Thermal runaway—a condition of uncontrolled temperature rise fueled by exothermic chemical reactions—poses a severe threat. Statistical analyses indicate that a majority of field failures originate from internal or external short circuits and electrical abuse scenarios like overcharge. Experimental investigation of these phenomena in modules is costly, hazardous, and offers limited insight into internal states. Therefore, high-fidelity multi-physics simulation has become an indispensable tool for probing the complex interplay of heat generation, transport, and chemical kinetics that leads to catastrophic failure. While models for single cells are well-established, simulating a pack requires capturing the thermal interactions between cells, which can dramatically accelerate failure propagation. This study addresses this need by constructing and validating a module-level simulation to unravel the thermal runaway mechanism under coupled electrical and thermal abuse.

1. Model Development and Simulation Setup

1.1 Geometric Model and Key Assumptions

The commercial 18650 cell features a complex internal jelly-roll structure. To make the module-level simulation computationally feasible while preserving the essential physics, the cell is simplified to a homogeneous cylinder with dimensions of 18 mm in diameter and 65 mm in height. The positive tab is modeled as a small cylinder, and the negative tab as a rectangular block, connected by busbars. The module of four cells in series is assembled accordingly. The following simplifying assumptions are made to formulate the problem:

  1. The effective material properties (density, specific heat) of the lithium-ion battery are homogeneous and independent of temperature and state of charge (SOC).
  2. The effective thermal conductivity is anisotropic, differentiated into axial (z-direction), radial, and tangential components, with axial and tangential conductivities being equal. These are also treated as constant.
  3. Heat generation within the cell volume is uniform; local hot spots from microstructure are not resolved.
  4. Internal convection within the electrolyte is negligible due to its limited mobility, and radiation heat transfer is omitted due to its minor contribution.

The computational domain is discretized using an O-grid meshing technique for the cylindrical cells. A mesh independence study confirmed that a mesh with over 197,000 elements provides a temperature field accuracy with a maximum relative error below 1.5%, ensuring a balance between fidelity and computational speed.

1.2 Thermophysical and Electrochemical Parameters

The effective properties for the homogenized lithium-ion battery model are derived from its constituent materials (positive electrode, negative electrode, separator, electrolyte, casing). The average density ($\rho$) is calculated from the cell mass and volume:
$$\rho = \frac{m_a}{V_a}$$
where $m_a$ is the cell mass (kg) and $V_a$ is the cell volume (m³).

The average specific heat capacity ($c$) is obtained via a mass-weighted average:
$$c = \frac{1}{m_a} \sum c_i m_i$$
where $c_i$ and $m_i$ are the specific heat (J/(kg·K)) and mass (kg) of component $i$.

The anisotropic thermal conductivities are calculated using an electrical resistance analogy, considering the series and parallel paths of heat flow through the layered structure. For a direction $j$ (x, y, z):
$$k_j = \frac{\sum_i k_i \delta_{j,i}}{L_j} = \frac{k_p \delta_{j,p} + k_n \delta_{j,n} + k_s \delta_{j,s}}{L_j}$$
where $k_p$, $k_n$, $k_s$ are the thermal conductivities of the positive electrode, negative electrode, and separator, respectively; $\delta_{j,p}$, $\delta_{j,n}$, $\delta_{j,s}$ are their thicknesses in direction $j$; and $L_j$ is the total cell dimension in that direction. The key parameters for the NCA18650B cell and the calculated effective properties are summarized in Tables 1 and 2.

Table 1: Key Specifications of the NCA18650B Lithium-ion Battery
Parameter Value
Nominal Capacity 3400 mAh
Nominal Voltage 3.6 V
Charge Cut-off Voltage 4.2 V
Discharge Cut-off Voltage 2.4 V
Maximum Continuous Charge Current 3.4 A (1C)
Maximum Continuous Discharge Current 6.8 A (2C)
Internal Impedance < 100 mΩ
Mass 48.5 g
Table 2: Calculated Effective Thermophysical Properties of the Battery Cell
Parameter Value
Average Density, $\rho$ 2766.3 kg/m³
Average Specific Heat, $c$ 1075.9 J/(kg·K)
Axial (Z) Thermal Conductivity, $k_z$ 29.85 W/(m·K)
Radial Thermal Conductivity, $k_r$ 1.47 W/(m·K)
Tangential Thermal Conductivity, $k_{\theta}$ 29.85 W/(m·K)

The heat generation rate within the lithium-ion battery during normal operation is described by the Bernardi model, which accounts for irreversible Joule heating and reversible entropic heat:
$$q = \frac{I}{V}\left[(U – U_1) – T\frac{dU_1}{dT}\right]$$
where $I$ is the current (A, positive for discharge), $V$ is the cell volume (m³), $U$ is the open-circuit voltage (V), $U_1$ is the terminal voltage (V), and $T$ is the temperature (K). For abuse conditions leading to thermal runaway, this model is coupled with a detailed kinetic model.

The cell’s state of charge (SOC) is tracked using the ampere-hour integration method:
$$S(t) = S_0 – \frac{\int_0^t I dt}{3600 C_0}$$
where $S(t)$ is the SOC at time $t$, $S_0$ is the initial SOC, and $C_0$ is the nominal capacity in Ampere-hours (Ah).

1.3 Thermal Runaway Kinetics and Boundary Conditions

To simulate the exothermic reactions during thermal abuse, a four-equation kinetic model is implemented. This model captures the major chain reactions that drive thermal runaway in a lithium-ion battery:

  1. SEI Decomposition: $$ \frac{dC_{SEI}}{dt} = -A_{SEI} \exp\left(-\frac{E_{SEI}}{RT}\right) C_{SEI}^{m_{SEI}} $$
  2. Negative Electrode with Electrolyte: $$ \frac{dC_{ne}}{dt} = A_{ne} \exp\left(-\frac{t_{SEI}}{t_{SEI,ref}}\right) \exp\left(-\frac{E_{ne}}{RT}\right) C_{ne}^{m_{ne}} $$
  3. Positive Electrode with Electrolyte: $$ \frac{d\alpha}{dt} = A_{pe} \exp\left(-\frac{E_{pe}}{RT}\right) \alpha^{m_{pe,1}} (1-\alpha)^{m_{pe,2}} $$
  4. Electrolyte Decomposition: $$ \frac{dC_{e}}{dt} = -A_{e} \exp\left(-\frac{E_{e}}{RT}\right) C_{e}^{m_{e}} $$

The variables $C_{SEI}$, $C_{ne}$, and $C_e$ represent the normalized remaining reactant fraction (from 1 to 0), while $\alpha$ represents the reaction progress (from 0 to 1). The total abuse heat generation rate is the sum of the heat released from each reaction:
$$ q_{abuse} = H_{SEI}\rho_{SEI}\left|\frac{dC_{SEI}}{dt}\right| + H_{ne}\rho_{ne}\left|\frac{dC_{ne}}{dt}\right| + H_{pe}\rho_{pe}\left|\frac{d\alpha}{dt}\right| + H_{e}\rho_{e}\left|\frac{dC_{e}}{dt}\right| $$
where $H_i$ and $\rho_i$ are the reaction enthalpies and densities of the respective components.

The initial temperature for all simulations is set to 25°C. A transient pressure-based solver is used with the standard k-ε turbulence model for any convective cooling (set to natural convection in this study). A time step of 1 second is used for accuracy.

1.4 Abuse Condition Simulation Setup

Three distinct abuse scenarios are simulated to trigger thermal runaway in the lithium-ion battery module:

  1. Internal Short Circuit (ISC): A small spherical region with a very low resistance ($5 \times 10^{-8} \Omega$) is defined at the geometric center of Cell #2 to simulate a localized internal short. The four-equation kinetic model is activated to simulate the ensuing exothermic reactions.
  2. Overcharge: The module is charged at constant current rates significantly above the maximum specification: 3C (10.2A), 5C (17A), and 7C (23.8A), until thermal runaway is initiated. The charging continues without a voltage cutoff, forcing the cells into deep overcharge.
  3. Overdischarge: The module is discharged at 1C (3.4A) to system voltages far below the specified cut-off of 9.6V (4 cells * 2.4V). The final system voltages are set to 2V, 1V, and 0V to simulate severe overdischarge conditions.

2. Simulation Results and Analysis

Thermal runaway is identified by key indicators: temperature exceeding a critical threshold (often >140°C), a rapid temperature rise rate (>1°C/s), abnormal voltage behavior, and the progression of the internal kinetic reactions.

2.1 Thermal Runaway Initiated by Internal Short Circuit

The simulation of the ISC condition reveals a clear cascade failure mechanism. The localized short in Cell #2 causes intense Joule heating, rapidly elevating its temperature beyond the SEI decomposition threshold. This triggers the exothermic chain reactions, leading to a full thermal runaway in Cell #2 within approximately 35 seconds, with a peak temperature exceeding 300°C and a maximum heating rate of 8.2°C/s.

The anisotropic thermal conductivity plays a crucial role in the failure propagation. The high axial conductivity (29.85 W/(m·K)) facilitates efficient heat transfer along the busbars and tabs. Consequently, Cell #1, directly connected via Busbar 1, is the first to be affected by the intense heat from Cell #2, leading to its own thermal runaway. Heat then propagates through Busbar 2 to Cell #3, and finally to Cell #4 via Busbar 3 and residual radiation/convection, resulting in the complete involvement of the module. The temperature evolution shows a sharp spike during the active reaction phase for each cell, followed by a cooling trend once the reactants are exhausted.

2.2 Thermal Runaway Initiated by Overcharging

The overcharge simulations demonstrate a strong dependence of thermal runaway severity and timing on the charging current (C-rate). As the C-rate increases from 3C to 7C, the time to trigger thermal runaway decreases by 58%. Furthermore, the peak temperature achieved during the runaway event is significantly higher at higher C-rates.

Analysis of the internal reaction kinetics provides deeper insight. The normalized reactant fractions for the SEI decomposition and electrolyte decomposition reactions are plotted over time. The results indicate that at a lower overcharge rate (3C), the temperature rise is more gradual. This slower heating delays the onset of the SEI decomposition reaction and slows down the overall reaction conversion rate. The electrolyte decomposition is also less severe.

Conversely, at a high overcharge rate (7C), the rapid joule heating quickly pushes the lithium-ion battery into the temperature regime where all exothermic reactions accelerate exponentially. The SEI layer decomposes faster, and the reaction between the cathode and the electrolyte proceeds more completely and rapidly. This leads to a faster total reaction conversion, releasing a massive amount of heat in a shorter time, which explains the earlier trigger time and higher peak temperature.

2.3 Thermal Runaway Initiated by Overdischarging

The overdischarge scenario presents a distinct behavior compared to ISC and overcharge. When the module is forced to discharge to deeply low voltages (0V, 1V, 2V system voltage), the peak temperature observed during the simulation is inversely related to the final voltage. Specifically, lowering the cut-off voltage from 2V to 0V results in a peak module temperature reduction of approximately 160°C.

This counter-intuitive result can be explained by the electrical behavior during deep overdischarge. As the cell voltage is driven far below its normal operating range, severe polarization occurs, and the effective internal resistance rises dramatically. Although detrimental side reactions (like copper dissolution from the anode current collector) occur, the sharply reduced current flow (due to the high effective resistance and low driving voltage) limits the total Joule heating. The heat generation from these side reactions is often time-dependent and may not accumulate sufficiently to trigger a classic, violent thermal runaway before the system voltage collapses entirely. Thus, the failure mode shifts from a thermal explosion to one characterized by rapid capacity fade, voltage collapse, and potentially a milder temperature rise.

3. Conclusion

This study presents a comprehensive simulation-based analysis of thermal runaway initiation and propagation in a series-connected lithium-ion battery module. By developing a coupled thermal-electrical model integrated with a four-equation chemical kinetics model, the complex failure mechanisms under internal short circuit, overcharge, and overdischarge conditions were elucidated. Key findings include:

  1. The anisotropic thermal properties of the cell and the thermal coupling via busbars are critical determinants of failure propagation speed and sequence following an internal short circuit.
  2. Overcharge severity is exponentially sensitive to the C-rate. Higher currents drastically reduce the time to failure and increase the peak runaway temperature by accelerating the kinetics of key exothermic reactions like SEI and electrolyte decomposition.
  3. Deep overdischarge leads to a different failure signature, where voltage collapse and limited current flow can result in a lower peak temperature compared to less severe overdischarge, highlighting the shift in dominant failure mechanisms.

The established simulation framework provides a powerful tool for probing the safety limits of lithium-ion battery packs. Future work will involve coupling this model with state-of-health (SOH) estimation algorithms and exploring its application to next-generation battery chemistries like sodium-ion and solid-state batteries. Furthermore, this model can directly inform the design of advanced thermal management systems and fault detection algorithms, ultimately contributing to the development of safer, more reliable energy storage systems.

Scroll to Top