Modern energy storage systems rely on densely packed lithium-ion battery modules, where a single cell failure can trigger cascading thermal runaway propagation (TRP) leading to catastrophic accidents. In this study, we investigate the thermal runaway (TR) behavior of 280 Ah LiFePO₄ energy storage cells from single-cell triggering to module‑level cascade propagation, and propose a physics‑constrained fault inversion localization framework. By combining controlled thermal abuse experiments, a data‑driven multiphysics model with forgetting factor recursive least squares (FFRLS) parameter identification, and a hybrid genetic algorithm‑grey wolf optimizer (GA‑GWO) enhanced time difference of arrival (TDOA) methodology, we achieve accurate prediction of TRP paths and robust real‑time source localization of the initiating cell. The key novelty lies in establishing a SOC‑dependent activation energy map that corrects the conventional Arrhenius kinetics, thus enabling precise simulation of energy storage cell cascade failures under varying state of charge (SOC). We further demonstrate that the proposed fault inversion approach maintains a localization accuracy above 94 % across different SOC levels, with an average directional dual‑threshold error (DD‑TEM) below 0.03, effectively eliminating cross‑layer misjudgment. This work provides a systematic framework for designing safer energy storage cabinets and for early warning of battery failures.
Experimental Investigation of Single‑Cell Thermal Runaway
To establish physicochemical boundary conditions for model development, we conducted thermal abuse experiments on 280 Ah prismatic LiFePO₄/graphite energy storage cells (dimensions 174 mm × 208 mm × 72 mm, weight 5435 ± 5 g, nominal voltage 3.2 V, operating voltage window 2.5‑3.65 V). All cells were from the same production batch with initial state of health (SOH) of 100 %. They were charged to target SOC levels (25 %, 50 %, 75 %, 100 %) using a constant‑current‑constant‑voltage protocol (IEC 62660‑1) with a NEWARE CT‑4008‑5V20A‑A cycler.
The cell was fixed in a steel fixture and heated by a 1000 W constant‑power heater attached to the cell’s large face. Five K‑type thermocouples (diameter 1 mm, measurement error ±2.5 °C) were placed at the center of the large face (Tₐ), top left and bottom right corners (Tb, Tc), the safety vent (Td), and 5 cm above the vent (Te). Temperature data were recorded by an ICPCON I‑7018 module, and cell voltage was logged by the cycler.
The experiment revealed a clear SOC‑dependent TR behavior. Table 1 summarizes the key triggering parameters. At 100 % SOC, the cell reached thermal runaway fastest (trigger time 469 s), exhibited a peak temperature of 378.3 °C, and a maximum temperature rise rate of 26.5 °C/s. In contrast, at 25 % SOC, the TR event was delayed to 2683 s, with a peak temperature of only 236.8 °C and a much slower heating rate (2.2 °C/s). The voltage collapse was abrupt for SOC ≥50 %, while for 25 % SOC it showed a gradual decay.
| Parameter | SOC 25 % | SOC 50 % | SOC 75 % | SOC 100 % |
|---|---|---|---|---|
| Vent opening time (s) | 1407 | 916 | 378 | 338 |
| Voltage drop time (s) | 1496 | 1430 | 984 | 327 |
| TR trigger time (s) | 2683 | 1698 | 1010 | 469 |
| Peak temperature time (s) | 2721 | 1799 | 1219 | 847 |
| Temperature at voltage drop (°C) | 136.5 | 136.6 | 114.9 | 61.2 |
| Temperature at vent opening (°C) | 107.9 | 101.7 | 70.3 | 73.8 |
| Temperature at TR trigger (°C) | 201.8 | 157.2 | 120.1 | 79.4 |
| Peak temperature (°C) | 236.8 | 238.4 | 350.9 | 378.3 |
| Max temperature rise rate (°C/s) | 2.2 | 3.9 | 16.4 | 26.5 |
These experimental data provided the calibration target for the subsequent multiphysics model, especially for identifying the activation energy of side reactions as functions of SOC.
Modeling Methodology: FFRLS‑Enhanced Multiphysics Simulation
Governing Equations for Thermal Runaway Side Reactions
The thermal runaway of a lithium‑ion energy storage cell is driven by four exothermic side reactions: SEI decomposition, negative‑electrolyte reaction, positive‑electrolyte reaction, and electrolyte decomposition. Each reaction rate follows the Arrhenius law:
$$ R_i = A_i \, c_i^{m_i} \, \exp\left(-\frac{E_{a,i}}{RT}\right) $$
where \(R_i\) is the reaction rate (s⁻¹), \(A_i\) the pre‑exponential factor (s⁻¹), \(c_i\) the dimensionless concentration of reactive species, \(m_i\) the reaction order, \(E_{a,i}\) the activation energy (J mol⁻¹), \(R\) the gas constant (8.314 J mol⁻¹ K⁻¹), and \(T\) the temperature. The heat generation from each reaction is:
$$ Q_i = H_i \, W_i \, R_i $$
where \(H_i\) is the specific heat release (J kg⁻¹) and \(W_i\) the reactive material loading (kg m⁻³). Table 2 lists the parameters used for the 280 Ah cell, based on literature and manufacturer data.
| Reaction | \(H_i\) (J kg⁻¹) | \(W_i\) (kg m⁻³) | \(A_i\) (s⁻¹) | \(c_0\) | \(m_i\) |
|---|---|---|---|---|---|
| SEI decomposition | 7.21 × 10⁵ | 413 | 1.67 × 10¹⁵ | 0.15 | 1 |
| Negative‑electrolyte | 8.99 × 10⁵ | 413 | 2.50 × 10¹³ | 0.75 | 1 |
| Positive‑electrolyte | 2.53 × 10⁵ | 925 | 6.67 × 10¹³ | 0.04 | 1 |
| Electrolyte decomposition | 1.60 × 10⁵ | 500 | 5.14 × 10²⁵ | 1.00 | 1 |
Energy Conservation in the Battery Domain
The temperature field inside the cell is governed by the heat conduction equation with internal heat generation:
$$ \rho C_p \frac{\partial T}{\partial t} = \nabla \cdot (\lambda \nabla T) + \frac{Q_{\text{side-tot}}}{V_{\text{cell}}} – \frac{Q_{\text{diss-tot}}}{V_{\text{batt}}} $$
where \(\rho\) is density (2118 kg m⁻³), \(C_p\) is specific heat capacity (1088 J kg⁻¹ K⁻¹), \(\lambda\) is anisotropic thermal conductivity (\(\lambda_x=1.5\), \(\lambda_y=\lambda_z=18\) W m⁻¹ K⁻¹), \(Q_{\text{side-tot}}\) is the total side‑reaction heat, \(Q_{\text{diss-tot}}\) accounts for convection and radiation losses, \(V_{\text{cell}}\) is the jellyroll volume, and \(V_{\text{batt}}\) is the total volume including casing. The dissipation term is:
$$ Q_{\text{diss-tot}} = h A_{\text{batt}} (T – T_{\text{amb}}) + \varepsilon \sigma A_{\text{batt}} (T^4 – T_{\text{amb}}^4) $$
with \(h=8.7\) W m⁻² K⁻¹, \(\varepsilon=0.8\), and \(T_{\text{amb}}=23\) °C.
FFRLS‑Based SOC‑Dependent Activation Energy
Conventional models assume constant activation energies, which leads to significant errors when SOC varies dynamically, especially at high SOC (>75 %). To overcome this, we propose a data‑driven parameter identification using the Forgetting Factor Recursive Least Squares (FFRLS) algorithm. The activation energy for each side reaction is expressed as a quadratic function of SOC:
$$ E_{a,i}(\text{SOC}) = k_{i,0} + k_{i,1} \cdot \text{SOC} + k_{i,2} \cdot \text{SOC}^2 $$
The FFRLS algorithm recursively updates the coefficient vector \(\boldsymbol{\theta} = [k_{i,0},\, k_{i,1},\, k_{i,2}]^\text{T}\) by minimizing the prediction error between simulated and experimental temperature evolution. The update rules are:
$$ \boldsymbol{\theta}(k+1) = \boldsymbol{\theta}(k) + \mathbf{K}(k+1) \, e(k+1) $$
$$ \mathbf{K}(k+1) = \frac{\mathbf{P}(k) \boldsymbol{\phi}(k+1)}{\lambda + \boldsymbol{\phi}^\text{T}(k+1) \mathbf{P}(k) \boldsymbol{\phi}(k+1)} $$
$$ \mathbf{P}(k+1) = \frac{1}{\lambda} \left[ \mathbf{I} – \mathbf{K}(k+1) \boldsymbol{\phi}^\text{T}(k+1) \right] \mathbf{P}(k) $$
where \(\boldsymbol{\phi}(k) = [1,\, \text{SOC}(k),\, \text{SOC}^2(k)]^\text{T}\) is the regression vector, \(\lambda=0.98\) is the forgetting factor, and \(e(k+1)\) is the temperature residual. Table 3 gives the identified activation energies for different SOC levels.
| SOC | SEI decomposition | Negative‑electrolyte | Positive‑electrolyte | Electrolyte decomposition |
|---|---|---|---|---|
| 25 % | 1.36 × 10⁵ | 1.42 × 10⁵ | 1.80 × 10⁵ | 1.86 × 10⁵ |
| 50 % | 1.28 × 10⁵ | 1.35 × 10⁵ | 1.71 × 10⁵ | 1.78 × 10⁵ |
| 75 % | 1.19 × 10⁵ | 1.25 × 10⁵ | 1.62 × 10⁵ | 1.69 × 10⁵ |
| 100 % | 1.14 × 10⁵ | 1.17 × 10⁵ | 1.56 × 10⁵ | 1.64 × 10⁵ |
The model was implemented in a 3D finite‑element framework using a non‑structured mixed mesh (approximately 946 139 elements after grid‑independence verification). The simulated temperature evolution at the center thermocouple location (Tₐ) matched the experimental data with a maximum peak‑temperature deviation of only 2.59 % across all SOC levels, demonstrating the effectiveness of the FFRLS calibration.
Module‑Level Thermal Runaway Propagation Analysis
Using the validated single‑cell model, we constructed a module comprising 52 cells (4 rows × 13 columns) arranged in series, with 0.5 mm aerogel insulation pads between cells and a water‑cooled bottom plate. The module was placed inside a cabinet‑like enclosure. We simulated TRP by triggering the centrally located cell #33 at t = 100 s via a localized heat source. Four distinct propagation phases were observed for SOC = 100 % and 75 %: single‑cell destabilization, neighbor diffusion, spatial propagation, and global failure. For SOC ≤ 50 %, no propagation occurred because the heat generated by the failing cell was insufficient to trigger adjacent cells.
The thermal propagation path was strongly anisotropic. Heat transfer in the y‑direction (along the electrode plane, facing the large sides of cells) was much faster than in the x‑direction (through the cell thickness). This is because the thermal conductivity in the y‑direction (\(\lambda_y = 18\) W m⁻¹ K⁻¹) is an order of magnitude higher than that in the x‑direction (\(\lambda_x = 1.5\) W m⁻¹ K⁻¹). For the 100 % SOC module, the average trigger interval between adjacent cells along the x‑axis was 73 s, while along the y‑axis it was only 66 s. For 75 % SOC, the corresponding intervals were 98.2 s and 89.5 s, respectively. The peak temperature reached 424.5 °C in the 100 % SOC module and 412.5 °C in the 75 % SOC module.
These findings establish a critical SOC threshold (somewhere between 50 % and 75 %) below which cascade propagation is self‑extinguishing. The path‑dependent heat flow also provides a physical constraint for the subsequent fault localization: the temperature rise at cabinet‑mounted sensors follows a deterministic sequence determined by the spatial arrangement and anisotropic conductivity of the energy storage cells.
Fault Inversion Localization: GA‑GWO Enhanced TDOA
Framework Overview
We developed a two‑stage fault inversion method to identify the precise location of the initial thermal runaway cell within an energy storage cabinet. The cabinet contains six diagnostic zones (each a sub‑region of the module array), with eight temperature sensors evenly distributed on the four side walls (two per side).
Stage 1 – Coarse localization. When a real thermal runaway event occurs, the sensor signals are denoised via wavelet thresholding, and the time delays relative to the first‑triggered sensor are extracted. These delays are matched against a pre‑computed database that contains the expected delay signatures for a fault initiated at the centroid of each zone. The zone with the minimum least‑squares error (\(L_{i,\text{loss}}\)) is selected as the candidate region. In our tests, the correct zone was always identified, with a maximum error of 0.945 for other zones and a minimum error of 0.163 for the true zone.
Stage 2 – Fine localization. Within the identified zone, we use a hybrid GA‑GWO algorithm to jointly optimize the (x, y, z) coordinates of the fault cell, its SOC, and SOH. The fitness function is based on the sum of squared differences between the measured and simulated sensor trigger delays:
$$ \text{Fitness}(\mathbf{x}) = \frac{1}{n} \sum_{j=1}^{n} \left( t_{i,j}^{\text{sim}} – t_{j}^{\text{meas}} \right)^2 $$
The GA component (crossover and mutation) expands the search diversity, while the GWO component (leadership hierarchy of α, β, δ wolves) refines the global optimum. We performed a sensitivity analysis to select the best hyperparameters: population size N = 30, maximum iterations = 75, linear decay of the convergence factor. The algorithm terminates when the improvement over 5 consecutive generations is below 0.01 %.
Error Metric
To evaluate localization accuracy, we defined the Directional Dual‑Threshold Error Metric (DD‑TEM). This metric accounts for the anisotropic propagation characteristics by applying different thresholds in the x (thickness) and y (in‑plane) directions:
$$ \text{DD} = \left( \sum_{i \in \{x,y,z\}} \frac{|p_i – a_i|}{\tau_i} \right) + \lambda \cdot \mathbb{I}\{|p_z – a_z| > \tau_z\} $$
where \(p\) is the predicted cell center and \(a\) the actual center; \(\tau_x=174\) mm, \(\tau_y=208\) mm, \(\tau_z=72\) mm correspond to the cell dimensions; \(\lambda=10\) penalizes cross‑layer errors (misjudgment in the z‑stack direction). A DD value below 0.03 indicates correct localization; above 0.1 indicates unacceptable error.
Localization Results
We benchmarked our GA‑GWO‑TDOA method against three alternatives: pure TDOA, PSO‑driven inversion, and single‑dimension GWO (optimizing only coordinates). Tests were performed for SOC levels of 25 %, 50 %, 75 %, and 100 %, with Gaussian noise added to simulate realistic sensor perturbations. Each case was repeated five times. The results are summarized in Table 4.
| SOC | TDOA only | PSO | Single‑dim GWO | Proposed GA‑GWO‑TDOA |
|---|---|---|---|---|
| 25 % | 0.127 | 0.053 | 0.079 | 0.027 |
| 50 % | 0.111 | 0.051 | 0.094 | 0.023 |
| 75 % | 0.107 | 0.043 | 0.040 | 0.021 |
| 100 % | 0.088 | 0.048 | 0.057 | 0.019 |
The proposed method consistently achieved DD < 0.03, corresponding to a localization accuracy above 94 %, and never produced a cross‑layer (z‑axis) misjudgment. In contrast, pure TDOA often misidentified the layer, especially at low SOC, while single‑dimension GWO suffered from large errors when the temperature gradient was weak.
We also evaluated the impact of sensor count and placement on localization accuracy (Table 5). Eight sensors uniformly distributed around the four cabinet sides provided an optimal balance: 95 % accuracy with DD = 0.021, whereas reducing to four sensors on only two sides still gave 90 % accuracy (DD = 0.041). This demonstrates the robustness of the method under sensor‑limited or failure scenarios.
| Number of sensors | Placement | Accuracy (%) | Avg. DD‑TEM |
|---|---|---|---|
| 2 | Left + front | 80 | 0.057 |
| 4 | Four sides | 90 | 0.025 |
| 4 | Left + front | 90 | 0.041 |
| 8 | Four sides | 95 | 0.021 |
| 16 | Four sides | 95 | 0.020 |
Conclusion
In this work, we have established a comprehensive methodology for modeling and localizing cascading thermal runaway failures in large‑format LiFePO₄ energy storage cells. Three principal contributions are highlighted:
- Experimental characterization of SOC‑dependent TR thresholds. Thermal abuse tests on 280 Ah cells revealed that the peak temperature varies from 236.8 °C at 25 % SOC to 378.3 °C at 100 % SOC, and the maximum temperature rise rate increases by more than an order of magnitude. These data provide essential boundary conditions for modeling.
- A data‑driven multiphysics model with FFRLS parameter identification. By incorporating SOC‑dependent activation energies into the Arrhenius kinetics, our model achieves a peak‑temperature prediction error below 3 % across all SOC levels. The module‑level simulations further identify the critical SOC threshold for propagation and quantify the anisotropic heat flow in the energy storage cell array.
- A robust GA‑GWO‑TDOA fault inversion framework. The two‑stage approach first narrows the fault to a cabinet zone using TDOA database matching, then precisely locates the triggering cell by joint optimization of spatial coordinates and state parameters. The method delivers 94 %+ accuracy with a DD‑TEM error below 0.03, even under sensor noise and reduced sensor counts. It successfully eliminates cross‑layer misjudgments that plague conventional methods.
Our findings directly support the design of safer energy storage cabinets: they indicate that maintaining SOC below a critical threshold (approximately 60 % for the studied cell) can inherently suppress cascade propagation, and that the proposed fault localization algorithm can be integrated into battery management systems for real‑time early warning. Future work will focus on extending the model to include thermal barrier optimization and multi‑cell internal short‑circuit scenarios, as well as validating the localization method on full‑scale energy storage racks.