Multiscale Simulation of Mechanically Induced Thermal Runaway in EV Battery Packs with Passive Safety Optimization

Chapter 1: Introduction and Research Framework

1.1 Background and Motivation

The global automotive industry is undergoing an unprecedented paradigm shift toward electrification. Against this backdrop, China’s new energy vehicle (NEV) sector has demonstrated remarkable technological leadership, with EV battery systems increasingly recognized as the core enabling technology. Nevertheless, despite the well-established environmental benefits and energy-efficiency advantages of electric vehicles, the critical challenge of thermal runaway in lithium-ion power batteries remains an unresolved bottleneck that threatens both public safety and market acceptance.

Official statistics and fire investigations indicate a year-on-year increase in electric vehicle fire incidents, rising from approximately 100 cases in 2021 to 141 cases in 2023. Among these incidents, the dominant contributing component was the EV battery system, with the principal root causes being thermal runaway of individual cells and electrical short circuits. While the proportion of fires caused by wiring faults has declined due to increasingly stringent electrical design standards in battery systems, the fundamental problem of cell-to-cell thermal runaway propagation has remained largely unsolved. This persistent safety challenge demands comprehensive investigation through advanced computational approaches.

1.2 Thermal Runaway Mechanisms in EV Battery Systems

In the context of EV battery engineering, thermal runaway phenomena are generally classified according to three categories of triggers, which can act independently or synergistically:

Mechanical abuse: External structural loads arising from vehicle collisions, bottom impacts, or battery pack crush events produce geometric deformation and internal component displacement, resulting in separator rupture and internal short circuit—the dominant pathway in many real-world EV battery fires.
Electrical abuse: Overcharge or over-discharge conditions induce irreversible phase transformations in electrode materials, leading to lithium dendrite growth that penetrates the separator, creating internal short-circuit channels.
Thermal abuse: Elevated ambient temperatures accelerate internal chemical reaction rates, causing self-heating that culminates in thermal runaway.

Critically, actual thermal runaway events in EV battery systems exhibit multi-factor coupling behavior. This complexity necessitates the development of high-fidelity multiphysics simulation frameworks capable of capturing the sequential chain of mechanical damage → electrical failure → thermal escalation → thermal runaway propagation.

1.3 Objectives and Research Scope

This research aims to establish a comprehensive multiphysics modeling and simulation framework for investigating mechanically induced thermal runaway in EV battery systems, followed by systematic passive safety optimization. The specific objectives are:

To develop and validate a high-fidelity electro-thermal coupled model of a prismatic lithium-ion cell using the Randles equivalent circuit approach with dynamic parameter optimization
To formulate an advanced thermal runaway criterion incorporating strain, temperature, and voltage as synergistic determination factors alongside multi-source heat generation (Joule heating and chemical reaction enthalpy)
To construct and analyze a full-scale EV battery system model subjected to bottom mechanical impact, examining mechanical-electrical-thermal coupled field responses
To propose and evaluate buffer and thermal insulation optimization strategies for passive safety enhancement of EV battery systems

Chapter 2: Theoretical Foundations of Lithium-ion Battery Simulation

2.1 Finite Element Method for Coupled Field Analysis

The finite element method (FEM) constitutes the fundamental numerical framework for this research. Originating from the aerospace structural analysis requirements of the 1940s, FEM has evolved through decades of development into the preeminent computational tool for engineering simulations. Within the displacement-based formulation central to structural analysis, the displacement field $ u(x,y,z) $ is discretized through shape functions $ N $ for nodal displacement vectors $ u_e $:

$$u(x, y, z) = N(x, y, z) \, u_e \tag{2.1}$$

$$u_x = \sum_{i=1}^{4} N_i u_{xi}, \quad u_y = \sum_{i=1}^{4} N_i u_{yi} \tag{2.2}$$

$$\varepsilon = B^T u_e \tag{2.3}$$

$$B = \begin{bmatrix} \frac{\partial u_x}{\partial x} & \frac{\partial u_y}{\partial y} & \frac{\partial u_x}{\partial y} + \frac{\partial u_y}{\partial x} \end{bmatrix} \tag{2.4}$$

$$\sigma = D \varepsilon \tag{2.5}$$

$$D = \frac{E}{1 – \nu^2} \begin{bmatrix} 1 & \nu & 0\\ \nu & 1 & 0\\ 0 & 0 & \frac{1-\nu}{2} \end{bmatrix} \tag{2.6}$$

For EV battery collision problems characterized by high-velocity transients, large deformations, high strain rates, and geometric nonlinearities, the explicit dynamics approach is preferred over implicit methods. The semi-discrete governing equation takes the form:

$$M\ddot{u} + C\dot{u} + Ku = F(t) \tag{2.7}$$

where $ M $, $ C $, and $ K $ represent the mass, damping, and stiffness matrices respectively, and $ F(t) $ is the external force vector. Explicit time integration advances the solution via:

$$\dot{u}_{t+\Delta t/2} = \dot{u}_{t-\Delta t/2} + \Delta t \ddot{u}_t \tag{2.8}$$

$$u_{t+\Delta t} = u_t + \Delta t \dot{u}_{t+\Delta t/2} \tag{2.9}$$

where nodal accelerations are computed directly from $\ddot{u} = M^{-1}(F – I)$, with $ I $ representing the internal force vector. The time step $ \Delta t $ must respect the Courant stability limit, making explicit analysis computationally efficient for high-speed impact scenarios.

2.2 Electrochemical Model Selection for EV Battery Analysis

The numerical simulation of lithium-ion batteries requires careful consideration of electrochemical model complexity. LS-DYNA provides three principal electrochemical modeling approaches, each characterized by distinct advantages and limitations:

Model Type	Theoretical Basis	Applications	Limitations
Newman P2D	Five coupled PDEs Doyle-Fuller-Newman theory	High-precision microscopic analysis High-rate cycling studies	Fine mesh requirements Extensive computational resources
Thermal Model	Macroscopic heat generation and transport	System-level temperature prediction	No electrochemical detail Limited failure simulation
Multiphysics Model	Integrated electro-thermal-mechanical coupling	Crash-induced thermal runaway analysis	Broad parameter requirements

For this research investigating thermal runaway triggered by mechanical impact in EV battery systems, the Multiphysics model was selected as it enables full-chain capture of the damage → short circuit → thermal runaway sequence. This approach balances the need for electrochemical fidelity against the computational demands of system-level simulation.

2.3 Contact and Collision Algorithms

EV battery systems contain numerous components whose mechanical interactions must be modeled with appropriate contact algorithms. LS-DYNA provides three principal contact formulations, each with distinct characteristics:

Kinematic constraint method — the earliest approach that enforces contact constraints by eliminating penetration of slave nodes through master surfaces, effectively modeling welded or bolted small-deformation connections.
Distributed parameter method — assigns half the slave surface mass to the contacted master surface, distributing internal stresses as surface pressures, primarily used for sliding interactions in large-deformation configurations.
Penalty function method — the most widely implemented approach, introducing interface spring forces proportional to penetration depth and contact stiffness. This algorithm conserves momentum, suppresses hourglassing, and may be tuned for stability, making it the standard choice for crash and impact analysis involving large deformations and sliding.

In our modeling strategy, the penalty function method was deployed for dynamic impact surfaces, while the kinematic constraint approach handled bolted or bonded joints in the battery pack architecture.

Chapter 3: Single-Cell Coupled-Field Simulation of EV Batteries

3.1 EM Solver and Joule Heating Formulation

The coupled-field simulation leveraged the Electromagnetic (EM) solver in LS-DYNA, which performs multi-field coupling between structural, thermal, and electromagnetic solvers. For electric-current-driven heating problems characteristic of thermal runaway scenarios, the EM solver operates on quasi-static electromagnetic equations:

The electrical conduction within battery conductors follows Ohm’s law under the quasi-static electric field assumption:

$$\mathbf{J} = \sigma \mathbf{E} = -\sigma \nabla \phi \tag{3.1}$$

where $ \mathbf{J} $ is the current density (A/m²), $ \sigma $ is the electrical conductivity (S/m), $ \mathbf{E} $ is the electric field intensity (V/m), and $ \phi $ is the electric potential (V).

The volumetric Joule heating source $ Q_{\text{joule}} $ is given by:

$$Q_{\text{joule}} = \frac{|\mathbf{J}|^2}{\sigma} \tag{3.2}$$

The resulting temperature distribution within the battery system is governed by the transient heat conduction equation:

$$\rho c_p \frac{\partial T}{\partial t} = \nabla \cdot (k\nabla T) + Q_{\text{joule}} \tag{3.3}$$

where $ \rho $ denotes density (kg/m³), $ c_p $ is the specific heat capacity (J/(kg·K)), $ k $ is the thermal conductivity (W/(m·K)), and $ T $ is the temperature (°C). Notably, the EM solver receives updated temperature fields from the thermal solver, enabling temperature-dependent characterization of electrical conductivity and Randles circuit parameters, thereby closing the electro-thermal feedback loop.

3.2 Distributed Randles Equivalent Circuit

The electrochemical behavior of lithium-ion cells was represented through a distributed Randles equivalent circuit model. This approach, schematically shown in the literature, connects positive and negative current collectors through a circuit governed by an ideal voltage source, internal resistance $ r_0 $, and n parallel RC branches. A first-order circuit (n = 1) with one RC branch provides sufficient accuracy for reproducing measured voltage damping characteristics:

$$\frac{dv_c}{dt} = \frac{i}{c_{10}} – \frac{v_c}{r_{10} c_{10}} \tag{3.4}$$

where $ v_c $ is the dynamic overpotential, $ i $ is the circuit current, and $ r_{10} $, $ c_{10} $ represent the charge-transfer resistance and double-layer capacitance of the RC branch, respectively. Key characteristics of the circuit parameterization include:

Open circuit voltage (OCV): Expressed as a function of state of charge (SOC) $ U = f(\text{SOC}) $
Internal resistance $ r_0 $: Captures the instantaneous Ohm’s law voltage jump during current steps
Charge-transfer resistance $ r_{10} $: Describes delayed, relaxation-type voltage response
Double-layer capacitance $ c_{10} $: Accounts for charge accumulation at electrode-electrolyte interface

The relationship between charge $ q $, capacity $ Q $, and SOC is given by:

$$q = Q \times 3600 \times \frac{\text{SOC}}{100} \tag{3.5}$$

$$q = \frac{Q}{C_Q} \times \text{SOC}, \quad C_Q = \frac{1}{36} \tag{3.6}$$

3.3 Dynamic Parameter Optimization of the Equivalent Circuit

Conventional Randles equivalent circuit characterization typically assumes fixed resistance and capacitance values. However, lithium-ion batteries operating in extreme environmental conditions—particularly those experienced during vehicle collisions or thermal excursions—exhibit substantial parameter drift. To capture these variations, we implemented a dynamic parameterization strategy based on HPPC (Hybrid Pulse Power Characterization) testing methodology.

The HPPC test procedure applies discrete charge/discharge current pulses at alternating SOC levels (typically 10% to 90% with 10% increments) across multiple temperature environments (25°C, 40°C, 50°C). From measured voltage responses, the parameters $ r_0 $, $ r_{10} $, and $ c_{10} $ are extracted as functions of SOC, temperature, and current direction. For example, the internal resistance $ r_0 $ was observed to vary with SOC and temperature as shown in the corresponding optimization curves.

Under transient thermal conditions, the circuit voltage response follows:

$$V = U + i \cdot r_0 + U_1 \tag{3.7}$$

$$U_1 = i \cdot r_{10}\left(1 – e^{-\frac{t}{r_{10}c_{10}}}\right) \tag{3.8}$$

This dynamic parameterization enabled accurate reproduction of cell behavior under rapidly varying thermal and electrical states, substantially improving simulation fidelity.

3.4 Finite Element Model of the Prismatic Cell

The single-cell reference geometry measured 173 mm × 21 mm × 85 mm, with 3.7 V nominal voltage and 100 Ah nominal capacity—specifications representative of cells used in modern EV battery traction packs.

Cell-level modeling considered the trade-off between fidelity and computational efficiency inherent to crashworthiness simulation. While detailed layered models precisely represent the internal architecture (current collectors, anode/cathode materials, separators), the computational overhead is prohibitive for system-level EV battery analysis. We therefore adopted a homogenized cell representation, which treats the cell core as an equivalent continuous medium while retaining the essential mechanical and electrochemical response characteristics.

The finite element discretization comprises the central cell domain modeled with 5 mm × 5 mm × 5 mm solid hexahedral elements (2,400 elements total), bounded by electrode terminal tabs. Mechanical behavior of the cell core was represented using a crushable foam constitutive model, which accurately captures the compaction behavior of layered jellyroll structures under compressive loading. The terminal tabs were represented through a linear elastic material model. Key material parameters for the cell core included: density of 2,223 kg/m³, elastic modulus of 1 × 10⁹ Pa, Poisson’s ratio of 0.05, tensile cutoff stress of 3 × 10⁹ Pa, and damping coefficient of 0.1.

The finite element model is further constrained through appropriate boundary conditions. The top and lateral faces of the cell are fully constrained in all translational and rotational degrees of freedom. The positive and negative terminal surfaces serve as equipotential boundaries for the electrical solution, with each cell node incorporating a first-order distributed Randles equivalent circuit. Electrical conductivities for the terminals and cell core materials were defined individually, thereby achieving comprehensive connectivity for the electrochemical simulation.

3.5 Mechanical Loading Scenarios

Three mutually orthogonal compressive loading configurations were investigated using a rigid cylindrical indenter (diameter = 20 mm, length = 120 mm). The indenter was driven at a constant velocity of 1 m/s along the X, Y, or Z directions over a compression depth of 13 mm. The contact formulation used a penalty-based approach with static and dynamic friction coefficients of 0.2.

These orthogonal loading orientations were selected to probe the anisotropic response of the homogenized cell, allowing systematic identification of the most critical loading direction for EV battery safety analysis.

3.6 Short-Circuit Criterion and Implementation

Lithium-ion cells in aluminum prismatic enclosures—the dominant architecture in EV battery applications—undergo catastrophic failure during mechanical abuse through a chain involving: (i) aluminum shell fracture and electrolyte leakage, (ii) separator rupture, (iii) direct internal short circuit between cathode and anode materials, and (iv) subsequent thermal runaway.

The internal short-circuit criterion in our model is based on the effective strain threshold derived from the mechanical deformation state. Aluminum shell elongation at break is approximately 0.08; when the local effective strain exceeds this critical threshold ($ \varepsilon > 0.08 $ ), the Randles equivalent circuit at the affected node is replaced by a pure resistance of 0.1 Ω:

$$\text{Short Circuit} = \begin{cases} \text{True}, & \varepsilon > 0.08 \\ \text{False}, & \varepsilon \leq 0.08 \end{cases} \tag{3.9}$$

This mechanical-electrical coupling was implemented through short-circuit criterion keywords and user-defined functional evaluation within the EM solver. The criterion evaluates the maximum principal strain at each material point, enabling spatial tracking of the short-circuit region.

3.7 Single-Cell Simulation Results and Analysis

3.7.1 Deformation and Strain Response

Simulation results revealed strongly direction-dependent mechanical response of the EV battery cell under compressive loading. The effective strain distributions (assessed through contour plots) displayed gradient patterns localized around the indentation zone. Peak strain values obtained for X, Y and Z directional compression equaled 0.32, 0.33, and 0.49, respectively. Importantly, the Z-direction (through-thickness) compression produced substantially higher peak strains than the X-direction or Y-direction cases, attributable to the anisotropic stiffness characteristics of the homogenized cell structure.

3.7.2 Electrical Response and Internal Short Circuit

Voltage distribution analysis showed highly localized potential drops corresponding to the mechanical compression zone. Under all three loading orientations, the minimum voltage reached approximately 0 V at the compression sites, while regions outside the deformation zone maintained potentials around 0.36–0.41 V. Since the standardized open-circuit voltage of a healthy cell was 4 V, the dramatic reduction to ~0.4 V following compression clearly evidenced internal short-circuit activation, consistent with the strain-based short-circuit criterion. Furthermore, the current density distribution varied spatially for compressed versus uncompressed zones. In particular, three material regions with distinct thermoelectric responses were observed in the Z-direction case:

Region A (strain 0.08–0.48): The highly compressed zone exceeding the critical strain threshold, thereby activating internal short-circuit conditions.
Region B (strain 0–0.08): Compressed but below the short-circuit threshold; cells here remained operational but experienced strong temperature rise through heat conduction.
Region C: Undeformed far-field region showing negligible current density or temperature increase.

3.7.3 Thermal Response

The temperature distributions directly correlated with the local current density fields, consistent with Joule heating as the dominant heat source during the short-circuit phase. For Z-direction compression, cell temperature increased from the initial 25°C to 311°C after 20 seconds, with a maximum heating rate of 81°C/s observed during the active compression phase. The extensive internal short circuit and electrolyte decomposition pathways at high strain levels accounted for the significantly higher temperature peak under Z-direction loading.

The transient electrical and thermal behavior satisfied the thermal runaway criteria specified in GB 38031-2020 for EV traction batteries:

Criterion Parameter	GB 38031-2020 Standard	Simulation Result	Status
Temperature rise rate	>1°C/s	81°C/s	✓ Exceeded
Voltage drop	Approaching 0 V	4 V → 0.4 V	✓ Confirmed
Current spike	Extremely high peak	2,500 A	✓ Confirmed
SOC depletion rate	Rapid decline	100% → 80% in 20 s	✓ Confirmed

3.8 Heat Generation Mechanisms in Lithium-ion Cells

The total heat generation rate of a lithium-ion EV battery during thermal abuse can be decomposed into four principal components:

$$Q_{\text{total}} = Q_p + Q_j + Q_r + Q_s \tag{3.10}$$

where:

$ Q_p = I^2 R_p $—polarization heat from concentration gradients
$ Q_j = I^2 R_j $—Joule heat (Ohmic heating)
$ Q_r $—reaction heat from reversible electrochemical processes
$ Q_s $—side-reaction heat from decomposition chemistries

During internal short circuit, $ Q_j $ dominates initially; however, once the cell temperature exceeds critical thresholds, exothermic decomposition reactions activate and become the principal heat source. As established by Ouyang and co-workers based on extensive adiabatic thermal runaway and scanning calorimetry testing, chemical reaction heat constitutes the dominant contribution to total energy release during the runaway cascade.

3.9 Thermal Runaway Chemistry Chain

Under mechanical abuse, the following chain of exothermic decomposition reactions drives the thermal runaway process:

Stage 1: SEI decomposition (80°C–120°C)—the initiation step of the runaway chain:

$$(\text{CH}_2\text{OCO}_2\text{Li})_2 \rightarrow \text{Li}_2\text{CO}_3 + \text{C}_2\text{H}_4 + \tfrac{1}{2}\text{O}_2 \tag{3.11}$$

$$2\text{Li} + (\text{CH}_2\text{OCO}_2\text{Li})_2 \rightarrow 2\text{Li}_2\text{CO}_3 + \text{C}_2\text{H}_4 \tag{3.12}$$

Stage 2: Electrolyte and lithiated anode reaction (~130°C):

$$\text{C}_3\text{H}_4\text{O}_3(\text{EC}) + 2\text{Li} \rightarrow \text{Li}_2\text{CO}_3 + \text{C}_2\text{H}_4 \tag{3.13}$$

$$\text{C}_5\text{H}_{10}\text{O}_3(\text{DEC}) + 2\text{Li} \rightarrow \text{Li}_2\text{CO}_3 + \text{C}_2\text{H}_4 + \text{C}_2\text{H}_6 \tag{3.14}$$

$$\text{C}_4\text{H}_6\text{O}_3(\text{PC}) + 2\text{Li} \rightarrow \text{Li}_2\text{CO}_3 + \text{C}_3\text{H}_6 \tag{3.15}$$

Stage 3: Cathode material thermal decomposition (~200°C)—releasing oxygen and accelerating electrolyte oxidation:

$$2\text{Li}(\text{NiCoMn})_x\text{O}_2 \rightarrow 2\text{Li}(\text{NiCoMn})_x\text{O}_{2-y} + \frac{y}{2}\text{O}_2 \tag{3.16}$$

Stage 4: Electrolyte thermal decomposition and combustion (>200°C)—the final and most exothermic stage:

$$\text{C}_3\text{H}_4\text{O}_3(\text{EC}) + \frac{5}{2}\text{O}_2 \rightarrow 3\text{CO}_2 + 2\text{H}_2\text{O} \tag{3.17}$$

$$\text{C}_4\text{H}_6\text{O}_3(\text{PC}) + 4\text{O}_2 \rightarrow 4\text{CO}_2 + 3\text{H}_2\text{O} \tag{3.18}$$

$$\text{C}_3\text{H}_6\text{O}_3(\text{DMC}) + 3\text{O}_2 \rightarrow 3\text{CO}_2 + 3\text{H}_2\text{O} \tag{3.19}$$

These reactions generate large quantities of flammable gases within the cell enclosure, increasing internal pressure until venting occurs—often with explosive ignition when the gas jet encounters ambient oxygen.

3.10 Enhanced Thermal Runaway Criterion for EV Battery Cells

While the initial modeling framework only captured electro-mechanical short-circuit-induced Joule heating, it neglected the critical chemical reaction enthalpy terms. In EV battery thermal runaway, chemical exothermic reactions release substantially greater energy than the internal short circuit itself, since Joule heating is merely the triggering mechanism. To rectify this deficiency, we formulated an enhanced thermal runaway criterion integrating multiple physical indicators.

The refined criterion is activated when two simultaneous conditions are satisfied:

$$\text{Chemical reaction source} = \begin{cases} \text{Active}, & T \geq T_{\text{SEI}} \text{ and } V < V_{\text{crit}} \\ \text{Inactive}, & \text{otherwise} \end{cases} \tag{3.20}$$

where $ T_{\text{SEI}} = 100\,^\circ\text{C} $ is the SEI decomposition threshold temperature, and $ V_{\text{crit}} = 0.1\,\text{V} $ is the voltage threshold identifying a shorted cell. The coupled criterion algorithm proceeds as:

$$\text{Thermal Runaway Condition} = \begin{cases} \text{Occurs}, & \varepsilon > 0.08 \text{ and } T \geq 100\,^\circ\text{C} \text{ and } V < 0.1\,\text{V} \\ \text{Does not occur}, & \text{otherwise} \end{cases} \tag{3.21}$$

This three-factor criterion (strain, temperature, voltage) provides a more sophisticated and accurate manifestation of the short circuit → thermal escalation → chemical decomposition pathway. Consequently, the coupled model incorporated heat contributions from both Joule heating $ Q_j = r_0 I^2 $ and chemical reaction sources, as illustrated by the optimized flow diagram.

Temperature contour comparisons between the standard and enhanced models revealed significant differences: the unmodified short-circuit-only model predicted temperatures of 206°C (5 s) and 284°C (15 s), whereas the full thermal-runaway model yielded 307°C (5 s) and 633°C (15 s), with a maximum asymptotic temperature around 750°C.

3.11 Model Validation

Experimental data from the literature on 100% SOC prismatic lithium-ion cell thermal runaway tests provided an independent basis for model validation. The physical experiments recorded maximum cell temperatures reaching 720°C—and the corresponding multiphysics simulation gave a peak temperature of 750°C, corresponding to a deviation within 5% of the experimental range.

Furthermore, the experimental measurement of cell voltage dropping rapidly to 0 V in conjunction with the temperature spike aligned well with our simulated voltage drop (from 4 V to 0.4 V) that accompanies thermal runaway. For comparison, additional simulations using the same coupled field framework were conducted at SOC levels below 70%, showing good agreement with experimental benchmarks that exhibit reduced peak runaway temperatures and delayed triggering. In summary,

Peak runaway temperature predicted within 40°C of experimental measurements
Voltage transient behavior consistent with experimental short-circuit dynamics
Temperature-rise-rate thresholds consistent with thermal runaway identification standards

Chapter 4: Battery System Coupled-Field Simulation for EV Applications

4.1 Heat Transfer Fundamentals in Battery Packs

Battery systems in EV applications comprise numerous cells in a compact layout with high thermal coupling between adjacent cells. Heat transfer from a cell undergoing thermal runaway to its neighbors involves three fundamental mechanisms:

Heat conduction (governed by Fourier’s law):

$$Q_{cd} = -kA_{cd}\frac{\Delta T}{\Delta x} \tag{4.1}$$

Heat convection (governed by Newton’s law of cooling):

$$Q_{cv} = hA_{cv}(T_s – T_f) \tag{4.2}$$

Heat radiation (governed by the Stefan-Boltzmann law):

$$Q_{rd} = \varepsilon \sigma A_{rd}(T_1^4 – T_2^4) \tag{4.3}$$

In dense EV battery pack configurations, heat conduction through the structural contact points is the dominant pathway for thermal propagation during the early stages, while radiation becomes increasingly significant at elevated temperatures above the runaway threshold.

4.2 Thermal Contact Modeling for Cell-to-Cell Heat Transfer

To validate the thermal coupling implementation, we first constructed a reduced three-cell model consisting of identical prismatic cells arranged in a line, with the central cell (Cell 2) subjected to a rigid spherical impact and adjacent cells (Cell 1 and Cell 3) serving as thermal receptors. The thermal contact between adjacent cells was defined through a surface-to-surface thermal contact formulation, where the heat-transfer coefficient is expressed as:

$$H_{cod} = \frac{c_f}{L} \tag{4.4}$$

$$H = \begin{cases} h_{tc}, & 0 \leq L \leq L_{\text{min}} \\ H_{cod}, & L_{\text{min}} < L \leq L_{\text{max}} \\ 0, & L > L_{\text{max}} \end{cases} \tag{4.5}$$

This contact formulation accounts for the radiation exchange between surfaces through the radiation factor $ f = \varepsilon \sigma $, and it properly models conduction through interstitial fluid layers whose gap thickness $ L $ may evolve dynamically during impact deformation.

Simulation results confirmed that when the central cell underwent internal short circuit and thermal runaway (reaching 716°C), significant heat transferred to both adjacent cells through the thermal contact surfaces. The inner faces of cells 1 and 3 showed marked temperature increases, while their outer surfaces remained near ambient—clearly indicating successful implementation of the conductive heat-transfer pathway between cells.

4.3 Full-Scale EV Battery System Model Development

The complete EV battery pack model measured 2062 mm × 1590 mm × 185 mm, containing 384 prismatic cells arranged in 8 groups. Each module group comprises 4 rows of 12 cells, separated by the longitudinal cross-members that also provide structural support and load distribution. Due to the complexity arising from combining mechanical, electrical, and thermal solvers, we performed a systematic model reduction:

Retained components: cover plate, bottom plate, housing shell, cross-members, cold plate, and all cell elements
Simplified components: bolts, busbars, thermal interface materials, wire harnesses, and small connectors—omitted or represented through equivalent contact conditions
Mesh strategy: 5 mm elements in the impact zone refined for high-resolution analysis, 10 mm elsewhere, resulting in 135,014 total elements

Table 4.1 — Material parameters for EV battery system components

Component	Material	Density (g/cm³)	Poisson’s Ratio	Elastic Modulus (MPa)	Yield Strength (MPa)	Tensile Strength (MPa)
Top cover	PCM	1.8	0.3	1,750	134.5	135.0
Housing	AL6061	2.7	0.3	70,000	260.0	313.3
Bottom plate	DP780	2.7	0.3	200,000	420.0	980.0
Cross-members	AL6061	2.7	0.3	70,000	260.0	313.3
Cold plate	AL3003	2.7	0.3	70,000	51.46	148.9

The boundary conditions fixed the outer frame in all degrees of freedom. Component interactions were treated through single-surface automatic contact, thermal contact between adjacent cells and other components, and offset-tied contact for permanent joints (e.g., cold plate to bottom plate, cells to cold plate) emulating bolted/bonded connections.

4.4 Bottom Impact Loading Configuration

Representative of a ground obstacle collision while driving, a rigid sphere with 100 mm diameter was positioned directly beneath the battery pack bottom plate and assigned an initial vertical velocity of 10 m/s. A deformable-to-rigid automatic transformation was defined after the initial deformation phase, allowing the model to transition from smaller-time-step electromechanical analyses during impact to longer battery thermal evolution afterward.

Due to computational resource constraints inherent to the scale of the EV battery model (hundreds of cells with mechanical-electrical-thermal coupling), only the impact zone cells contained the distributed Randles equivalent circuits, while the remaining cells were treated as passive thermal elements.

4.5 Battery System Simulation Results

4.5.1 Displacement and Strain Response

Upon application of the bottom impact to the EV battery system, only a localized region (identified as Zone A, directly beneath and surrounding the impact point) showed substantial displacement, reaching a maximum of 0.064 m. The remaining portion of the pack, designated Zone B, experienced negligible displacement. Zone A bottom plates attained a peak strain of 0.56, and more importantly the cells within this zone showed compressive strain values up to 0.4—dramatically exceeding the critical short-circuit threshold of 0.08.

4.5.2 Electrical and Thermal Response of the Battery System

Cells within Zone A, particularly those in the central affected module group, satisfied the strain-based short-circuit criterion. Their voltages collapsed from 4 V to essentially 0 V, SOC dropped from 100% toward 0%, and the currents exhibited large impulsive peaks as stored energy was dissipated through the short-circuit resistance.

The simulated temperature evolution of the affected cells within the EV battery system described a severe thermal cascade. The figure shows that the temperature spikes sharply during the early stage, with peak temperature reaching 845°C and a maximum heating rate of 164°C/s. The battery system thus exhibited unequivocal thermal runaway under bottom impact conditions. Furthermore, the results confirmed that the heat produced by the thermal runaway of one group of cells propagates toward adjacent groups through conduction and radiation pathways—posing a continued fire risk even after the initially impacted module has fully reacted.

Chapter 5: Passive Safety Optimization of the EV Battery System

5.1 Optimization Strategies for Mechanical, Thermal, and Electrical Protection

Given the physical chain of thermal runaway in EV battery systems—mechanical deformation → separator failure → internal short circuit → Joule heating → chemical decomposition cascade—passive safety design can accordingly be addressed along three parallel routes:

Mechanical safeguard: the addition of energy-absorbing layers within the structural layout to reduce mechanical strain on cells under impact loading
Electrical safeguard: high-voltage circuit protection via insulation and rapid circuit interruption (outside the scope of this study)
Thermal safeguard: incorporation of thermal insulation barriers to slow or arrest the propagation of thermal runaway between neighboring cells or module assemblies

5.2 Materials Selection: Buffer Layer

Various classes of structural cushioning materials were evaluated for their suitability as an EV battery system buffer. Important characteristics include crush resistance, thermal tolerance, density, fatigue resistance and manufacturing cost. After eliminating inadequate candidates:

Material Category	Representative Materials	Key Advantages	Critical Disadvantages
Polyurethane foam (PU)	PU foam	Excellent vibration damping	Poor high-temperature resistance; limited compressive strength
Polyethylene foam (PE)	PE foam	Low cost and lightweight	Insufficient thermal stability
Polypropylene foam (MPP)	Expanded polypropylene	Good mechanical strength, high temperature tolerance, corrosion resistance	Moderate cost
Thermoplastic polyurethane (TPU)	TPU elastomer	High elasticity and flexural performance	Density higher than foam; thermal limits
Fiber-reinforced composites	GFRP, CFRP	Very high stiffness-to-weight	Brittle; expensive for CFRP versions

Based on this evaluation, MPP foam was chosen as the buffer layer material. It combines moderate density with robust energy absorption, stable load-bearing characteristics, and appropriate thermal stability up to temperatures well beyond normal operating range—making it especially suitable for bottom-impact protection in EV battery packs.

5.3 Placement and Thickness Parametric Investigation

The buffer layer was located between the bottom plate and the cold plate of the battery system—the primary load path for bottom impacts. The passive protection configuration is conceived so that upon ground collision events, the bottom plate transmits force into the buffer layer, which crushes inelastically to absorb kinetic energy and disperse the loading footprint before it reaches the cells. Buffers of 10 mm and 20 mm thickness were compared under identical loading conditions (same boundary conditions, material properties, solver configuration and impact velocity).

5.4 Simulation Results of the Bottom Buffer Optimization

Application of the MPP buffer layer led to measurable attenuation of impact severity indicators:

Bottom plate peak strain reduction: 10% (10 mm) and 17% (20 mm) relative to the non-buffered configuration
Cell maximum effective strain reduction: 15% (10 mm) and 25% (20 mm) relative to baseline values of 0.4 (no buffer),
Cell temperature reduction due to shortened and weakened short-circuit zones
Active temperature response, plotted in Fig. 5.7, showed a reduced thermal peak.

Table 5.1 — Cell module performance comparison under varying cushioning configurations

Optimization Case	Module I Peak Temp (°C)	Module II Peak Temp (°C)	Module III Peak Temp (°C)	Module II Heating Rate (°C/s)
No buffer	39.0	845	43.1	164
10 mm buffer	36.9	811	40.4	157
20 mm buffer	32.8	791	37.4	153

Compared to the no-buffer case, the use of 10 mm and 20 mm thickness buffer layers reduced peak temperatures at the central thermal runaway module by 34°C and 54°C, respectively. Temperature values for the adjacent modules likewise diminished, reflecting lower conductive heat flux from the damaged module. Since the strain values (0.34 and 0.30) still exceeded the short-circuit threshold, the central cells continued to undergo thermal runaway—but with curtailed intensity. It may thus be concluded that within an EV battery context, the buffer effectively reduces—but does not fully suppress—thermal runaway severity, because the mechanical deformation required to trigger separator failure is too great to be fully accommodated by MPP foam thicknesses of practical relevance in the 10 mm–20 mm range.

5.5 Material Selection: Thermal Insulation Layer

To mitigate the inter-module propagation of thermal runaway, thermal insulation materials were evaluated according to the following requirements: low thermal conductivity, high temperature capability, low density, mechanical durability during crash, chemical compatibility with battery electrolyte, and manufacturing cost. A summary comparison of materials considered:

Material Category	Representative Materials	Thermal Advantages	Critical Disadvantages
Aerogel composites	SiO₂ aerogel felt, carbon aerogel, ceramic aerogel	Very low conductivity, stable at high temperatures	Production cost high for some variants; limited mechanical strength
Ceramic fibers	Aluminosilicate, polycrystalline mullite fibers	Excellent thermal stability up to ~1400°C	Fragility, health hazards during handling, cost
Mica sheets	Natural or synthetic mica	Good electrical and thermal insulation	Brittleness and high density
Intumescent materials	Graphite based sheets	Large volume expansion on heating	Networking and fill constraints; toxic off-gassing in some types
Polymer-based insulators	FRP, polyimide foams, thermal pads	Lightweight, good dielectric properties	Melting or decomposition below ~300°C

Among the various insulating options considered for the EV battery system, a silica aerogel blanket stood out thanks to its optimum combination of very low thermal conductivity, fire resistance, high operating temperature class, low density, and design adaptability. It was selected as the internal thermal insulator for module-to-module interfaces.

5.6 Thermal Barrier Placement and Simulation Results

Silica aerogel sheets were modeled as being mounted between the thermally active central module (Module II in the impacted zone) and the adjacent outer module columns (Mod I and Mod III). For this analysis the buffer configuration (with 10 mm MPP) was retained in addition to the insulation to allow the appraisal of both mechanisms independently. The arrangement was implemented through contact tie constraints that allowed heat transfer through the insulation while preserving its structural position.

Table 5.2 — Comparison of thermal performance with and without thermal insulation

Optimization Case	Module I Peak Temp (°C)	Module I Max Heating Rate (°C/s)	Module III Peak Temp (°C)	Module III Max Heating Rate (°C/s)
No buffer and no insulation	39.0	0.42	43.1	1.06
With buffer (10 mm)	36.9	0.32	40.4	0.66
With buffer (10 mm) plus insulation	33.1	0.24	34.4	0.36

The thermal insulation brought an additional reduction of roughly 15% in the peak temperatures of adjacent Modules I and III compared to the buffered-only case, together with lower heating rates (by 25% and 45%, respectively). Most notably, the heating rate of Module III declined from 0.66°C/s to 0.36°C/s—below the threshold of active thermal acceleration, indicating a reduced probability that neighboring cells will run away.

The coupled buffer-and-insulation strategy thus proves successful in achieving a simultaneous mitigation of direct impact-induced cell damage and cell-to-cell heat propagation. In practical EV battery systems, these lightweight porous barrier materials simultaneously fulfill two key functions: (1) absorbing crash energy through progressive compaction, and (2) interrupting the chain of thermal propagation that would otherwise extend the danger to the entire battery pack and potentially the vehicle occupant compartment.

Additional comparative evaluations of the model configurations (electrical and chemical heat generation enabled, with and without the optimized passive safety layers) confirmed that the central module continues to run away, but the thermal gradients are reduced, temperatures are lowered, and the fire risk to adjacent modules becomes significantly less acute.

Chapter 6: Conclusions and Future Perspectives

6.1 Summary and Conclusions

This research has presented a comprehensive mechanically-induced thermal runaway simulation framework from single-cell to full EV battery system level, accompanied by optimization of the passive safety layout. Key findings and technical advances can be summarized as follows:

A robust single-cell multi-physics simulation model of EV battery cells was established. Through comparison of available electrochemical models for LS-DYNA, the distributed Randles circuit method combined with the homogenized BatMac macro model demonstrated the best trade-off among accuracy, model size, and computational cost. The model reproduces the fundamental short-circuit → temperature escalation → chemical decomposition response.
Dynamic parameterization of equivalent circuit elements was introduced. In contrast to earlier static Randles approaches, the resistance and capacitance terms were expressed as functions of SOC and temperature, allowing accurate transients under extreme operating conditions and strongly time-varying fields.
Enhanced thermal runaway criterion function was formulated. The strain criterion alone can only trigger Joule heating. We incorporated both $ T \geq 100^\circ\text{C} $ (SEI decomposition onset) and $ V < 0.1\,V $ (short-circuit voltage state) for prediction of chemical-reaction heat generation, establishing a three-variable criterion that yields accurate thermal evolution. The predicted maximum temperature of 750°C matched experimental prismatic cell data to within 40°C—a deviation of under 5%.
A system-level EV battery pack model was developed using 384 cells and an explicit dynamics finite element mesh of 135,014 elements. Under bottom impact by a 100 mm sphere at a velocity of 10 m/s, the local cell assemblies in the impact zone reached peak effective strains of 0.4, causing a severe rise in temperature to 845°C with peak heating rates of up to 164°C/s. In contrast, strain levels in remote regions remained below the short-circuit threshold, confirming the localized yet potentially propagative character of thermal runaway.
Passive safety optimization studies were conducted for EV battery systems. A 20 mm MPP buffer layer located below the cells is shown to reduce maximum effective strain by 25%, lower the thermal runaway cell temperature from 845°C to 791°C, and simultaneously reduce bottom-plate strain by roughly 17%. SiO₂ aerogel thermal barriers placed between module groups effectively suppressed heat propagation, lowering the peak temperature of adjacent non-failed modules by 10%–20% and reducing their maximum heating rate coefficients below critical runaway threshold values.

6.2 Limitations and Future Work

Despite the progress made in this study, several limitations and open questions present opportunities for future investigations:

Direct experimental validation: A direct comparison with dedicated bottom-impact and crush tests of the specific cell architecture analyzed would further increase confidence. Due to facility limitations, the model was validated against available published thermal runaway test data rather than original testing. A future effort should combine abuse testing with simulation.
Expanded criterion function: The present three-factor (strain-temperature-voltage) thermal runaway criterion captures primary thermal pathway mechanisms. In real-world thermal runaway conditions in EV battery systems, additional parameters—including time at elevated temperature, SOC dependence in a larger range, heating rates, and pressure—may further increase fidelity. Multi-criteria, probabilistic and data-driven thermal runaway models may be combined to more closely reproduce thermal runaway arrest and propagation trends observed experimentally.
Advanced materials and multi-layer optimization: Further development of novel buffer/insulation materials (e.g., phase-change thermal management integrated with aerogel barriers) can potentially improve both energy absorption and thermal pathway interruption simultaneously.
Coupling with other safety systems: Structural protection in a real EV should be co-designed together with active monitoring systems (BMS), venting or fire-suppression mechanisms, so that conventional and passive safety complement one another.

With the rapid development of high-energy-density EV batteries, this integrated mechanical-electrical-thermal approach is becoming mandatory to guarantee safety in next-generation EVs. The combined use of accurate multi-physics simulation and simultaneous optimization of mechanical, electrical, and thermal layers will continue to serve as an essential engineering pathway for the development of inherently safe EV battery systems.