In the pursuit of clean and renewable energy solutions, electrochemical energy storage systems have become pivotal. Among these, the LiFePO4 battery stands out due to its high energy density, long cycle life, and cost-effectiveness. However, as the scale of energy storage systems expands and application scenarios diversify, the issue of thermal runaway in LiFePO4 batteries poses a significant bottleneck to safety and reliability. In particular, overcharging is a common trigger for thermal runaway in lithium-ion battery energy storage stations. This phenomenon can lead to lithium dendrite growth, separator penetration, internal short circuits, and ultimately, catastrophic failure. Therefore, understanding and mitigating thermal runaway in LiFePO4 batteries is crucial for advancing safe energy storage technologies.
This study employs simulation methods to delve into the thermal runaway behavior of LiFePO4 batteries under overcharge conditions. By establishing a three-dimensional electrochemical-thermal coupling model, we investigate the evolution of various battery components during thermal runaway and perform multi-parameter analysis and optimization. Our research aims to provide theoretical insights for preventing and suppressing overcharge-induced thermal runaway in LiFePO4 batteries, thereby enhancing their safety profile in real-world applications.
To simulate the complex processes within a LiFePO4 battery, we developed an integrated model that combines electrochemical reactions, heat generation, and thermal runaway kinetics. The model is based on the Newman framework, which incorporates mass conservation, charge conservation, and energy conservation principles. The governing equations are as follows.
Solid-phase charge conservation is described by:
$$
\nabla \cdot \mathbf{i}_s + \nabla \cdot (-\sigma \nabla \phi_s) = 0
$$
where $\nabla$ is the Laplace operator, $\mathbf{i}_s$ is the solid-phase current density vector, $\phi_s$ is the solid-phase potential, and $\sigma$ is the solid-phase electrical conductivity.
Solid-phase mass conservation is given by:
$$
\frac{\partial c_s}{\partial t} = \nabla \cdot (D_s \nabla c_s)
$$
where $c_s$ is the lithium-ion concentration in the electrode material, and $D_s$ is the solid-phase lithium-ion diffusion coefficient.
Liquid-phase mass conservation is expressed as:
$$
\frac{\partial c_e}{\partial t} = \nabla \cdot (D_e \nabla c_e) + \frac{1 – t_+^0}{F} \nabla \cdot \mathbf{i}_e
$$
where $c_e$ is the lithium-ion concentration in the electrolyte, $D_e$ is the liquid-phase diffusion coefficient, $t_+^0$ is the transference number of lithium ions, and $F$ is Faraday’s constant.
Liquid-phase charge conservation is:
$$
\nabla \cdot \mathbf{i}_e + \nabla \cdot \left( -k \nabla \phi_e – \frac{2kRT}{F} \left(1 + \frac{\partial \ln f_\pm}{\partial \ln c_e}\right) (t_+^0 – 1) \nabla \ln c_e \right) = 0
$$
where $\mathbf{i}_e$ is the liquid-phase current density, $k$ is the ionic conductivity, $\phi_e$ is the liquid-phase potential, $R$ is the universal gas constant, $T$ is the temperature, and $f_\pm$ is the mean molar activity coefficient.
The charge transfer at the solid-liquid interface is governed by the Butler-Volmer equation:
$$
j = a i_0 \left[ \exp\left(\frac{(1-\alpha)F}{RT} \eta\right) – \exp\left(-\frac{\alpha F}{RT} \eta\right) \right]
$$
where $j$ is the local current density, $a$ is the specific surface area, $i_0$ is the exchange current density, $\alpha$ is the charge transfer coefficient, and $\eta$ is the overpotential.
Energy conservation accounts for heat generation from electrochemical reactions, irreversible polarization, ohmic heating, and side reactions, as well as heat dissipation through convection and radiation. The equation is:
$$
\frac{\partial (\rho c_p T)}{\partial t} = \nabla \cdot (\lambda \nabla T) + q_{\text{rea}} + q_{\text{act}} + q_{\text{ohm}} + q_s – q_{\text{dis}}
$$
where $\rho$ is the density, $c_p$ is the specific heat capacity, $\lambda$ is the thermal conductivity, $q_{\text{rea}}$ is the electrochemical reaction heat, $q_{\text{act}}$ is the activation polarization heat, $q_{\text{ohm}}$ is the ohmic heat, $q_s$ is the side reaction heat, and $q_{\text{dis}}$ is the heat dissipation rate. The convective heat dissipation is given by $q_{\text{dis}} = h A \Delta T$, where $h$ is the convective heat transfer coefficient, $A$ is the surface area, and $\Delta T$ is the temperature difference.
For thermal runaway modeling, we consider key exothermic side reactions: Solid Electrolyte Interphase (SEI) decomposition, negative electrode-electrolyte reaction, positive electrode-electrolyte reaction, and electrolyte decomposition. The total side reaction heat is:
$$
q_s = q_{\text{sei}} + q_{\text{ne}} + q_{\text{pe}} + q_{\text{e}}
$$
Each reaction rate follows an Arrhenius expression:
$$
r_x(T, c_x) = A_x \cdot c_x \exp\left(-\frac{E_{a,x}}{RT}\right)
$$
where $x$ denotes the reaction type, $A_x$ is the pre-exponential factor, $c_x$ is the normalized concentration of active material, and $E_{a,x}$ is the activation energy. The heat generated from each reaction is:
$$
q_x = h_x w_x r_x
$$
where $h_x$ is the heat released per unit mass, and $w_x$ is the mass content of active material. The concentration changes are described by:
$$
\frac{dc_{\text{sei}}}{dt} = -r_{\text{sei}}, \quad \frac{dc_{\text{ne}}}{dt} = -r_{\text{ne}}, \quad \frac{d\alpha}{dt} = r_{\text{pos}}, \quad \frac{dc_{\text{e}}}{dt} = -r_{\text{e}}
$$
Additionally, lithium plating on the negative electrode during overcharging is modeled by incorporating plating kinetics into the electrochemical model, where the local current density at the graphite surface is the sum of lithium intercalation and plating reactions.
The parameters for the electrochemical and thermal runaway models are summarized in the tables below. These values are critical for accurately simulating the behavior of LiFePO4 batteries under overcharge conditions.
| Parameter | Positive Electrode | Separator | Negative Electrode |
|---|---|---|---|
| Thickness $L$ (μm) | 90 | 20 | 80 |
| Solid-phase volume fraction $\epsilon_s$ | 0.513 | 0.360 | 0.580 |
| Liquid-phase volume fraction $\epsilon_l$ | 0.347 | 0.540 | 0.369 |
| Electrode particle radius $r$ (μm) | 1 | — | 15 |
| Max active material concentration $C_{s,\text{max}}$ (mol/m³) | 22,806 | — | 31,370 |
| Solid-phase Li diffusion coefficient $D_s$ (m²/s) | 8 × 10⁻¹⁶ | — | 1.08 × 10⁻¹⁴ |
| Initial electrolyte concentration $c_l$ (mol/m³) | 1,200 | 1,200 | 1,200 |
| Variable $x$ | Pre-exponential Factor $A_x$ (s⁻¹) | Activation Energy $E_{a,x}$ (J/mol) | Normalized Concentration $c_{x,0}$ | Energy per Unit Mass $h_x$ (J/g) | Active Material Content $w_x$ (g/m³) |
|---|---|---|---|---|---|
| SEI decomposition | 1.667 × 10¹⁵ | 1.3508 × 10⁵ | 0.15 | 257 | 6.100 × 10⁻⁵ |
| Positive electrode reaction | 6.66 × 10¹³ | 1.2254 × 10⁵ | 0.04 | 314 | 0.122 × 10⁻⁵ |
| Negative electrode reaction | 2.50 × 10¹³ | 1.3508 × 10⁵ | 0.75 | 1,714 | 6.100 × 10⁻⁵ |
| Electrolyte decomposition | 5.14 × 10²⁵ | 2.7400 × 10⁵ | 1.00 | 155 | 4.069 × 10⁻⁵ |
The simulation couples the one-dimensional electrochemical model with a three-dimensional thermal model. The electrochemical model computes reaction rates and heat generation, which are passed to the thermal model to calculate temperature distribution. The temperature feedback influences the electrochemical reactions and triggers the thermal runaway equations when critical temperatures are reached. This dynamic coupling allows for a comprehensive analysis of the LiFePO4 battery’s behavior under overcharge.
To validate our model, we compared simulation results with experimental data under adiabatic conditions at a 1C charging rate. The temperature profiles during heating and thermal runaway showed good agreement, confirming the model’s accuracy. Furthermore, prior studies using similar models for LiFePO4 batteries have reported consistent results, reinforcing the reliability of our approach.
We first analyzed the thermal runaway process of a LiFePO4 battery under a 0.5C overcharge scenario. The voltage and temperature evolution revealed distinct phases. Initially, the voltage increased gradually, and the temperature rise was moderate. After approximately 120 minutes, the temperature rise accelerated significantly, indicating the onset of overcharging. During this phase, lithium plating on the negative electrode increased internal resistance, leading to additional heat generation. When the battery temperature exceeded 400 K, an irreversible thermal runaway was triggered, characterized by a sharp temperature spike to about 540 K within a short period. The heat generation power surged to over 8,000 W/m³ during this event. After reaching the peak temperature, the side reaction intensities diminished due to material depletion, and the battery began to cool due to heat dissipation to the environment. However, the damage was already irreversible, compromising the safety and performance of the LiFePO4 battery.
The contributions of different battery components to heat generation during thermal runaway were examined. The SEI decomposition reaction initiated first, releasing heat early in the process. This was followed by reactions between the negative electrode and electrolyte, and subsequently, the positive electrode and electrolyte. The electrolyte decomposition became significant during the intense thermal runaway stage. The heat released from the positive and negative electrode reactions was substantial and prolonged, making them primary heat sources. The reaction extents and heating rates varied: SEI decomposition reached completion with a peak heating rate of 0.07 K/s; the positive electrode reaction achieved full conversion with a peak rate of 2.3 K/s; the negative electrode reaction had a lower extent of about 40% but a peak rate of 0.6 K/s; and the electrolyte decomposition completed with a peak rate of 0.17 K/s.
The temperature distribution within the LiFePO4 battery during different stages was also investigated. In the normal charging phase (0–7,200 s), heat generation was slow, and temperature distribution was relatively uniform, with slight internal heating. During overcharging (7,200–23,000 s), heat generation increased, leading to larger temperature gradients, though convection helped maintain some uniformity. In the full thermal runaway phase (after 23,000 s), intense internal heating dominated, resulting in high temperatures throughout the battery, with only superficial cooling at the surfaces due to convection.
Next, we explored the impact of operating conditions on thermal runaway in LiFePO4 batteries. By varying charging rates, convective heat transfer coefficients, and ambient temperatures, we identified strategies to delay thermal runaway and mitigate risks.
Under a constant ambient temperature of 298 K and heat transfer coefficient of 8 W/(m²·K), different charging rates (0.5C, 1C, and 2C) were applied. Higher charging rates led to earlier onset of thermal runaway due to increased internal heat generation from resistance and reaction kinetics. Specifically, at 2C, thermal runaway occurred at 3,050 s; at 1C, at 7,540 s; and at 0.5C, at 23,000 s. The peak thermal runaway temperatures were similar across rates, as the underlying mechanisms were consistent. However, the time to reach thermal runaway was significantly affected, emphasizing the importance of controlling charging rates in LiFePO4 battery systems.
We then analyzed the effect of convective heat transfer, modeled by varying the heat transfer coefficient (4, 8, and 12 W/(m²·K)) at a 0.5C charging rate and 298 K ambient temperature. A higher heat transfer coefficient delayed thermal runaway and reduced peak temperatures. For instance, at 4 W/(m²·K), thermal runaway occurred at 17,680 s with a peak temperature of 577 K; at 8 W/(m²·K), at 23,000 s with 540 K; and at 12 W/(m²·K), at 30,570 s with 498 K. Enhanced cooling thus extends the safe operating window and lowers the severity of thermal runaway in LiFePO4 batteries.
Ambient temperature variations (288 K, 298 K, and 308 K) were studied with a 0.5C charging rate and 8 W/(m²·K) heat transfer coefficient. Higher ambient temperatures accelerated thermal runaway onset due to faster reaction kinetics and reduced heat dissipation efficiency. Each 10 K increase in ambient temperature advanced thermal runaway by 250–350 s. Peak temperatures were slightly higher at elevated ambient temperatures (differences of 2–4 K). Post-thermal runaway cooling was also slower at higher ambient temperatures due to smaller temperature gradients.
These findings highlight that optimizing operating conditions—such as reducing charging rates, improving cooling systems, and managing ambient temperatures—can effectively delay thermal runaway and enhance the safety of LiFePO4 batteries.
To further improve LiFePO4 battery safety, we conducted an orthogonal experimental design to optimize electrode component parameters. The goal was to maximize thermal runaway time and battery capacity while considering multiple factors. The selected parameters included positive electrode particle diameter ($D_{\text{pos}}$), positive electrode thickness ($L_{\text{pos}}$), positive electrode porosity ($\epsilon_{\text{pos}}$), separator thickness ($L_{\text{film}}$), negative electrode particle diameter ($D_{\text{neg}}$), negative electrode thickness ($L_{\text{neg}}$), and negative electrode porosity ($\epsilon_{\text{neg}}$). Separator porosity was excluded due to its minimal impact based on preliminary analysis.
An L18(3⁶) orthogonal array was used, with three levels for each factor. The experiments simulated constant-current charging at 24 A/m². The response variables were thermal runaway time and specific capacity (Ah/m²). The results are summarized in the following tables.
| Run | $D_{\text{pos}}$ (m) | $L_{\text{pos}}$ (m) | $\epsilon_{\text{pos}}$ | $L_{\text{film}}$ (m) | $D_{\text{neg}}$ (m) | $L_{\text{neg}}$ (m) | $\epsilon_{\text{neg}}$ | Thermal Runaway Time (s) | Capacity (Ah/m²) |
|---|---|---|---|---|---|---|---|---|---|
| 1 | 2×10⁻⁷ | 7.0×10⁻⁵ | 0.30 | 1×10⁻⁵ | 0.5×10⁻⁵ | 6.5×10⁻⁵ | 0.32 | 10,012 | 20.06 |
| 2 | 2×10⁻⁷ | 9.0×10⁻⁵ | 0.35 | 2×10⁻⁵ | 1.5×10⁻⁵ | 8.0×10⁻⁵ | 0.37 | 7,979 | 23.67 |
| 3 | 2×10⁻⁷ | 1.1×10⁻⁴ | 0.40 | 3×10⁻⁵ | 2.5×10⁻⁵ | 1.0×10⁻⁴ | 0.42 | 7,769 | 26.40 |
| 4 | 1×10⁻⁶ | 7.0×10⁻⁵ | 0.30 | 2×10⁻⁵ | 1.5×10⁻⁵ | 1.0×10⁻⁴ | 0.42 | 8,248 | 16.87 |
| 5 | 1×10⁻⁶ | 9.0×10⁻⁵ | 0.35 | 3×10⁻⁵ | 2.5×10⁻⁵ | 6.5×10⁻⁵ | 0.32 | 5,399 | 23.13 |
| 6 | 1×10⁻⁶ | 1.1×10⁻⁴ | 0.40 | 1×10⁻⁵ | 0.5×10⁻⁵ | 8.0×10⁻⁵ | 0.37 | 12,462 | 29.00 |
| 7 | 2×10⁻⁶ | 7.0×10⁻⁵ | 0.35 | 1×10⁻⁵ | 2.5×10⁻⁵ | 8.0×10⁻⁵ | 0.42 | 5,303 | 16.73 |
| 8 | 2×10⁻⁶ | 9.0×10⁻⁵ | 0.40 | 2×10⁻⁵ | 0.5×10⁻⁵ | 1.0×10⁻⁴ | 0.32 | 13,488 | 25.79 |
| 9 | 2×10⁻⁶ | 1.1×10⁻⁴ | 0.30 | 3×10⁻⁵ | 1.5×10⁻⁵ | 6.5×10⁻⁵ | 0.37 | 6,581 | 27.07 |
| 10 | 2×10⁻⁷ | 7.0×10⁻⁵ | 0.40 | 3×10⁻⁵ | 1.5×10⁻⁵ | 8.0×10⁻⁵ | 0.32 | 9,225 | 20.05 |
| 11 | 2×10⁻⁷ | 9.0×10⁻⁵ | 0.30 | 1×10⁻⁵ | 2.5×10⁻⁵ | 1.0×10⁻⁴ | 0.37 | 8,350 | 23.67 |
| 12 | 2×10⁻⁷ | 1.1×10⁻⁴ | 0.35 | 2×10⁻⁵ | 0.5×10⁻⁵ | 6.5×10⁻⁵ | 0.42 | 11,519 | 26.40 |
| 13 | 1×10⁻⁶ | 7.0×10⁻⁵ | 0.35 | 3×10⁻⁵ | 0.5×10⁻⁵ | 1.0×10⁻⁴ | 0.37 | 13,917 | 18.48 |
| 14 | 1×10⁻⁶ | 9.0×10⁻⁵ | 0.40 | 1×10⁻⁵ | 1.5×10⁻⁵ | 6.5×10⁻⁵ | 0.42 | 5,564 | 21.59 |
| 15 | 1×10⁻⁶ | 1.1×10⁻⁴ | 0.30 | 2×10⁻⁵ | 2.5×10⁻⁵ | 8.0×10⁻⁵ | 0.32 | 7,483 | 30.47 |
| 16 | 2×10⁻⁶ | 7.0×10⁻⁵ | 0.40 | 2×10⁻⁵ | 2.5×10⁻⁵ | 6.5×10⁻⁵ | 0.37 | 4,795 | 18.27 |
| 17 | 2×10⁻⁶ | 9.0×10⁻⁵ | 0.30 | 3×10⁻⁵ | 0.5×10⁻⁵ | 8.0×10⁻⁵ | 0.42 | 12,050 | 21.67 |
| 18 | 2×10⁻⁶ | 1.1×10⁻⁴ | 0.35 | 1×10⁻⁵ | 1.5×10⁻⁵ | 1.0×10⁻⁴ | 0.32 | 9,584 | 31.45 |
The mean responses for each factor level were computed to identify optimal conditions. For thermal runaway time, the best combination to maximize time was $D_{\text{pos}} = 2 \times 10^{-6}$ m, $L_{\text{pos}} = 7.0 \times 10^{-5}$ m, $\epsilon_{\text{pos}} = 0.35$, $L_{\text{film}} = 3 \times 10^{-5}$ m, $D_{\text{neg}} = 0.5 \times 10^{-5}$ m, $L_{\text{neg}} = 1.0 \times 10^{-4}$ m, and $\epsilon_{\text{neg}} = 0.37$, yielding 13,917 s. For capacity, the best combination was $D_{\text{pos}} = 1 \times 10^{-6}$ m, $L_{\text{pos}} = 1.1 \times 10^{-4}$ m, $\epsilon_{\text{pos}} = 0.30$, $L_{\text{film}} = 2 \times 10^{-5}$ m, $D_{\text{neg}} = 2.5 \times 10^{-5}$ m, $L_{\text{neg}} = 8.0 \times 10^{-5}$ m, and $\epsilon_{\text{neg}} = 0.32$, yielding 30.467 Ah/m².
Considering both safety and performance, a balanced optimal configuration was determined: $D_{\text{pos}} = 2 \times 10^{-6}$ m, $L_{\text{pos}} = 9 \times 10^{-5}$ m, $\epsilon_{\text{pos}} = 0.40$, $L_{\text{film}} = 2 \times 10^{-5}$ m, $D_{\text{neg}} = 5 \times 10^{-6}$ m, $L_{\text{neg}} = 1 \times 10^{-4}$ m, and $\epsilon_{\text{neg}} = 0.32$. This design resulted in a thermal runaway time of 13,488 s and a specific capacity of 25.786 Ah/m², which is comparable to typical LiFePO4 battery capacities.
Range analysis and ANOVA were performed to assess factor significance. For thermal runaway time, the order of influence from most to least was: negative electrode particle diameter > negative electrode thickness > negative electrode porosity > positive electrode thickness > separator thickness > positive electrode particle diameter > positive electrode porosity. Negative electrode particle diameter and thickness were statistically significant factors. For capacity, the order was: positive electrode thickness > negative electrode porosity > negative electrode thickness > separator thickness > negative electrode particle diameter > positive electrode particle diameter > positive electrode porosity, with positive electrode thickness and negative electrode porosity being significant.
These results indicate that increasing positive electrode thickness, separator thickness, and negative electrode thickness generally delays thermal runaway by enhancing thermal mass and heat conduction paths. In contrast, larger particle diameters and higher negative electrode porosity tend to accelerate thermal runaway due to increased reaction rates and reduced stability. Positive electrode porosity showed minimal effect. For capacity, thicker positive electrodes increase active material content, boosting capacity, while higher negative electrode porosity reduces it slightly.
In conclusion, this study provides a detailed analysis of thermal runaway in LiFePO4 batteries through simulation and optimization. Key findings include: (1) The irreversible thermal runaway threshold for LiFePO4 batteries is approximately 400 K. (2) Thermal runaway initiates with SEI decomposition, followed by vigorous reactions at the positive and negative electrodes, which are primary heat sources. (3) Operating conditions such as lower charging rates, reduced ambient temperatures, and enhanced convective cooling can significantly delay thermal runaway and mitigate its severity. (4) Orthogonal optimization of electrode parameters reveals that negative electrode characteristics (particle diameter and thickness) are critical for thermal runaway time, while positive electrode thickness and negative electrode porosity dominate capacity. An optimized design was identified to balance safety and performance. These insights offer valuable guidance for designing safer LiFePO4 battery systems, contributing to the advancement of reliable energy storage solutions. Future work could extend this model to battery pack levels or incorporate real-time monitoring strategies for early thermal runaway detection.