What is battery thermal runaway, and can it be predicted via simulation?

In lithium-ion batteries, when the internal temperature exceeds approximately 130°C, the SEI layer begins to decompose. Above 200°C, exothermic reactions between the cathode and electrolyte can lead to runaway. Using FEM-based coupled electrochemical-thermal analysis, the propagation time can be predicted, which aids in the design of thermal barriers.

What equations are used in the heat generation model for battery thermal runaway?

The governing equations primarily consist of the Arrhenius-type reaction rate equation, k = A·exp(-Ea/RT), and the energy conservation equation using heat generation Q = I²R + I·T·(dU/dT). A standard four-stage reaction model—SEI decomposition, anode-electrolyte reaction, cathode decomposition, and electrolyte decomposition—is typically employed.

How is thermal propagation between cells modeled?

Modeling treats a single cell's thermal runaway as the trigger. The temperature propagation throughout the entire module is then calculated using 3D transient FEM analysis, which accounts for heat conduction, radiation, and convection between cells. Parametric analysis, varying the thickness and material of thermal barrier components, is used for fire propagation prevention design.

Can simulations substitute for UL 2580 or GB 38031 testing?

While not a complete substitute, virtual testing is utilized for pre-screening test conditions, optimizing thermal insulation design, and identifying worst-case scenarios. Physical testing remains mandatory for standard certification, but simulation significantly contributes to reducing the number of required tests.

Battery Thermal Runaway Simulation

Category: 連成解析 > 電気化学連成 | 更新 2026-04-12

Theory and Physics

🧑‍🎓

Is battery thermal runaway the thing that causes explosions? Can it be predicted with simulation?

🎓

In simple terms, when the internal temperature of a lithium-ion battery exceeds about 130°C, the SEI (Solid Electrolyte Interphase) layer begins to decompose. This exposes the negative electrode lithium, which reacts directly with the electrolyte, generating heat and further increasing the temperature. Above 200°C, oxygen release from the positive electrode and exothermic reactions of the electrolyte can run away, potentially exceeding 800°C in a matter of seconds.

🧑‍🎓

What, 800°C in seconds!? That would be dangerous if it happened in a battery installed in a car...

🎓

Exactly. If thermal runaway in one cell propagates heat to adjacent cells, the fire can spread chain-reaction style through the entire module → entire pack. This is "thermal propagation," the main cause of EV fires. That's precisely why it's crucial to use CAE to predict propagation time and optimize the thickness of insulation barriers and cooling design.

Thermal Runaway Mechanism

🧑‍🎓

Could you explain the "stages" of temperature increase in a bit more detail?

🎓

The temperature rise progresses through a 4-stage chain reaction. In simulation models, it's standard to represent each of these four stages with an Arrhenius equation.

Stage	Temperature Range	Reaction	Heat Generation Estimate
Stage 1	80–120°C	SEI layer decomposition	~250 J/g
Stage 2	120–200°C	Negative electrode-electrolyte reaction (exposure of intercalated Li)	~1700 J/g
Stage 3	180–250°C	Positive electrode decomposition (differs significantly between NMC, LFP)	~800 J/g
Stage 4	200°C–	Electrolyte decomposition (vaporization/combustion of organic solvent)	~400 J/g

🧑‍🎓

So, the decomposition of the SEI layer acts as a trigger, setting off a chain of reactions. Is LFP considered safer than NMC because Stage 3 positive electrode decomposition is milder?

🎓

Exactly. LiFePO₄ (LFP) has strong Fe-O bonds and a high oxygen release temperature. Therefore, heat generation from the positive electrode is low, making it less prone to thermal runaway. However, it's not "absolutely never runs away"; cases exist where external short circuits or overcharging force progression to Stage 2 and beyond. In simulations, we model this by changing the Arrhenius constants for each positive electrode material.

Governing Equations

🧑‍🎓

What are the governing equations for thermal runaway simulation?

🎓

The basis is the energy conservation equation (unsteady heat conduction equation) with added heat generation terms from chemical reactions and Joule heating:

$$ \rho c_p \frac{\partial T}{\partial t} = \nabla \cdot (k \nabla T) + Q_{\text{rxn}} + Q_{\text{ohm}} $$

Here, $\rho$ is density, $c_p$ is specific heat, $k$ is thermal conductivity. The first term on the right side is heat conduction, $Q_{\text{rxn}}$ is chemical reaction heat generation, and $Q_{\text{ohm}}$ is Joule heating.

Arrhenius Reaction Rate Model

🧑‍🎓

What's inside $Q_{\text{rxn}}$? I often hear about the Arrhenius equation...

🎓

We describe the reaction rate for each stage $i$ using the Arrhenius equation, multiply by the reaction heat $H_i$, and sum them up. This is $Q_{\text{rxn}}$:

$$ Q_{\text{rxn}} = \sum_{i} H_i \cdot A_i \exp\!\left(-\frac{E_{a,i}}{RT}\right) c_i^{n_i} $$

The meaning of each parameter is as follows:

$A_i$: Frequency factor (pre-exponential factor) [1/s]
$E_{a,i}$: Activation energy [J/mol]
$R$: Gas constant = 8.314 J/(mol·K)
$c_i$: Dimensionless concentration of reactant (0–1, decreases as consumed)
$n_i$: Reaction order
$H_i$: Reaction enthalpy [J/m³]

🧑‍🎓

The $\exp(-E_a/RT)$ part is key, right? As temperature $T$ increases, the reaction rate increases exponentially, generating more heat, which further increases temperature... That's exactly the runaway loop!

🎓

Yes, that's the essence of "thermal explosion" according to Semenov theory. The moment heat dissipation can't keep up, temperature skyrockets due to positive feedback.

Stage	$E_a$ [kJ/mol]	$A$ [1/s]	$H$ [J/g]
SEI decomposition	135	1.67×10¹⁵	257
Negative electrode-electrolyte	110–140	2.50×10¹³	1714
Positive electrode decomposition (NMC)	170–200	1.75×10⁹	792
Electrolyte decomposition	140–160	5.14×10²⁵	420

These values are typically obtained by fitting experimental data from DSC (Differential Scanning Calorimetry) or ARC (Accelerating Rate Calorimetry).

Battery Heat Generation Model

🧑‍🎓

What about $Q_{\text{ohm}}$, the Joule heating? Is the heat generation during normal operation different from during runaway?

🎓

The heat generation model during normal operation is called the Bernardi model. This is also important as an initial condition for thermal runaway simulation:

$$ Q_{\text{cell}} = I^2 R_{\text{int}} + I \cdot T \frac{dU_{\text{OCV}}}{dT} $$

First term $I^2 R_{\text{int}}$: Joule heating (resistive heating). Always positive, larger with higher current. Internal resistance $R_{\text{int}}$ is a function of SOC and temperature.
Second term $I \cdot T \cdot dU_{\text{OCV}}/dT$: Entropic heating (reversible heat). Sign changes with charge/discharge. $dU_{\text{OCV}}/dT$ is the temperature coefficient of OCV (Open Circuit Voltage).

For example, during fast charging (2C–3C), the first term (Joule heating) dominates, and cell temperature can rise to 60–80°C. If cooling is insufficient, it can approach the trigger temperature for SEI decomposition.

🧑‍🎓

I see, so we combine the heat generation model for normal operation with the Arrhenius reaction model for runaway.

Inter-Cell Thermal Propagation Model

🧑‍🎓

After one cell runs away, how is the fire spreading to adjacent cells modeled?

🎓

For the inter-cell propagation model, we consider all three heat transfer modes:

$$ q_{\text{prop}} = q_{\text{cond}} + q_{\text{conv}} + q_{\text{rad}} $$

Conduction $q_{\text{cond}}$: Conduction through direct contact between cells or via barrier materials. Most dominant. Controlled by the thermal conductivity $k$ and thickness $d$ of the barrier material.
Convection $q_{\text{conv}}$: Forced convection inside the module due to venting (ejection) of high-temperature gas. Vent gas temperature can reach 600–1000°C.
Radiation $q_{\text{rad}}$: Stefan-Boltzmann radiation proportional to $\sigma T^4$. Cannot be ignored above 800°C.

In practice, we set one "trigger cell" (e.g., nail penetration or internal short circuit) and evaluate the time until the temperature of adjacent cells exceeds the thermal runaway initiation temperature as the "propagation time." The design criterion is whether this propagation time can ensure occupant escape time (typically 5 minutes or more).

🧑‍🎓

It seems like changing the barrier material thickness just a little could significantly alter the propagation time. That's a perfect problem for simulation.

🎓

Yes. In one automotive manufacturer's case, changing a ceramic fiber-based barrier from 5mm to 8mm increased the propagation time from about 40 seconds to over 120 seconds. That 80-second difference can determine occupant safety. This kind of parametric analysis is FEM's forte.

Coffee Break Casual Talk

The "Trigger" of Thermal Runaway – What Happens at 130°C?

The key to understanding thermal runaway in lithium-ion batteries is the "chain" of temperature and reactions. SEI film decomposition begins at 80–120°C, separator melting and internal short circuits occur at 130–150°C, oxygen release from the positive electrode happens at 180–200°C—these exothermic reactions trigger each other in succession, and temperature can exceed 800°C in seconds. However, the timing of heat generation strongly depends on manufacturing variations and degradation state of the battery, so "runaway starts at the same temperature for every battery" is not true. Modeling this individual variation is one of the most difficult challenges in thermal runaway simulation.

Physical Meaning of Each Term

Heat storage term $\rho c_p \partial T / \partial t$: Represents the thermal inertia of battery materials (positive electrode, negative electrode, electrolyte, current collectors). For jelly-roll structures, anisotropy exists, so different effective heat capacities may be set for the in-plane and thickness directions.
Heat conduction term $\nabla \cdot (k \nabla T)$: Cylindrical cells exhibit strong anisotropy, with thermal conductivity differing by over 40 times between the in-plane direction ($k_r \approx$ 30 W/(m·K)) and the radial direction ($k_z \approx$ 0.7 W/(m·K)). Pouch cells also have large differences between in-plane and thickness directions.
Chemical reaction term $Q_{\text{rxn}}$: Four-stage Arrhenius-type reaction model. The positive feedback of temperature rise → reaction acceleration → further temperature rise is the essence of thermal runaway.
Electrical heat generation term $Q_{\text{ohm}}$: Described by the Bernardi model ($I^2R + I T \, dU/dT$) during normal operation. During an internal short circuit, localized high-current Joule heating through the short-circuit resistance dominates.

Assumptions and Applicability Limits

Homogenization assumption: The electrode layered structure inside the jelly-roll is treated as a homogeneous anisotropic material. Local reaction non-uniformity inside the cell is ignored.
Arrhenius-type reaction rate: Each stage is approximated by a simple first-order reaction. Actual reactions have multiple parallel paths, but an effective model can be obtained by fitting to DSC/ARC data.
Reactant concentration assumption: Tracks consumption of dimensionless concentration $c_i$ (0–1). Sometimes expressed as unconverted fraction.
Internal short circuit model: For reproducing nail penetration tests, the resistance value and contact area of the short-circuit location greatly affect results. Experimental calibration is essential.
Applicability limits: Fluid dynamics after gas venting (vent gas behavior) must be solved separately with CFD or given as a simplified convective heat transfer coefficient.

Dimensional Analysis and Unit System

Variable	SI Unit	Typical Value / Notes
Activation energy $E_a$	J/mol	SEI: 135 kJ/mol, Positive electrode: 170–200 kJ/mol
Frequency factor $A$	1/s	Order of 10⁹–10²⁵. Varies greatly by stage
Reaction enthalpy $H$	J/kg or J/m³	Note conversion between volume and mass basis
Internal resistance $R_{\text{int}}$	Ω	SOC and temperature dependent. Around 1–100 mΩ
Inter-cell thermal conductivity	W/(m·K)	Barrier material: 0.03–0.5, Aluminum casing: 150–200

Numerical Methods and Implementation

🧑‍🎓

I'm starting to understand the theory. How do you actually solve these equations on a computer?

🎓

Since it's an unsteady electrochemistry-thermal coupled problem, the standard approach is to use FEM (Finite Element Method) for spatial discretization and implicit schemes (Backward Euler or Crank-Nicolson method) for time discretization. The strong nonlinearity of the Arrhenius term is numerically tricky; this is the unique difficulty of thermal runaway simulation.

FEM Formulation

🧑‍🎓

What does the FEM formulation look like? Is FEM for heat conduction similar to structural analysis?

🎓

Up to approximating the temperature field with shape functions $N_i$ and transforming to the weak form, it's the same as standard heat conduction FEM. Applying the Galerkin method yields:

$$ [C]\left\{\dot{T}\right\} + [K]\{T\} = \{Q(\{T\})\} $$

$[C] = \int_\Omega \rho c_p [N]^T [N] \, d\Omega$: Heat capacity matrix
$[K] = \int_\Omega [B]^T k [B] \, d\Omega$: Thermal conductivity matrix
$\{Q\}$: Heat generation vector (integrated $Q_{\text{rxn}} + Q_{\text{ohm}}$)

The key point is that the right-hand side $\{Q\}$ depends strongly on temperature $\{T\}$. Because of the Arrhenius term $\exp(-E_a/RT)$, the heat generation becomes a nonlinear function of temperature. Therefore, an iterative method (Newton-Raphson) is essential.

Time Integration Scheme

🧑‍🎓

Temperature changes rapidly once runaway starts, right? How do you decide the time step?

🎓

This is the point practitioners struggle with most. Before thermal runaway, $\Delta t = 1\text{–}10$ seconds is sufficient, but after runaway begins, temperature gradients become steep, so $\Delta t$ must be reduced to 0.001–0.01 seconds to keep up.

Scheme	Stability	Accuracy	Recommended Scenario
Backward Euler (1st order implicit)	Unconditionally stable	$O(\Delta t)$	Robust but has large diffusion. Suitable for initial screening
Crank-Nicolson (2nd order implicit)	Unconditionally stable	$O(\Delta t^2)$	Focus on accuracy. Suitable for gentle pre-runaway heating
BDF2 (2nd order backward differentiation)	A-stable	$O(\Delta t^2)$	Strong for stiff systems. Effective during runaway transition

In practice, adaptive time step control is essential. If the temperature change $|\Delta T|$ exceeds a threshold (e.g., 5°C/step), automatically halve $\Delta t$, and increase it during gentle phases. Software like COMSOL Multiphysics has this feature available by default.

Handling Nonlinear Coupling

🧑‍🎓

You mentioned converging with the Newton-Raphson method. Are there issues like difficulty converging because the Arrhenius term is an exponential function?

🎓

Sharp question. The sensitivity of the Arrhenius term to temperature is extremely high. Because the Jacobian (tangent stiffness matrix) includes the temperature derivative of the reaction rate:

$$ \frac{\partial Q_{\text{rxn}}}{\partial T} = \sum_i H_i A_i c_i^{n_i} \frac{E_{a,i}}{RT^2} \exp\!\left(-\frac{E_{a,i}}{RT}\right) $$

In the runaway region, $\partial Q/\partial T$ becomes very large, worsening the condition number of the Jacobian. Countermeasures include:

Subdividing load steps: Automatically add substeps upon detecting runaway initiation.
Line search method: Apply damping to Newton update amounts to prevent overshoot.
Clamping reactant concentration: Force $c_i$ to stay within [0, 1] to prevent unphysical values.

Mesh Strategy

🧑‍🎓

What should we be careful about when meshing? Batteries have a jelly-roll structure with hundreds of layers, right?

🎓

In practice, meshing every single layer of the jelly-roll individually is computationally impossible. Therefore, we use the following multi-stage approach:

Scale	Model	Typical Element Count	Mesh Type
Electrode layer (μm)	1D Electrochemical Model (P2D Model)	100–1,000	1D line elements
Single cell (mm)	Homogeneous anisotropic 3D FEM	10,000–50,000	Hexahedral