Thermal-Electrical-Mechanical Coupled Finite Element Models For Battery Electric Vehicle

1. Introduction

With the increasing number of battery electric vehicles (BEV) on public roadways, there is a need to develop a finite element simulation method to fully understand the crash safety of BEVs and their lithium-ion batteries. Due to the relative low energy density of lithium-ion battery, BEVs are typically 125% heavier than their Internal Combustion Engine (ICE) equivalent [1,2]. This additional mass has negative effect on safety performance and fleet compatibility [3]. ICE vehicles and BEVs differ significantly in terms of mass distribution and structural arrangement, particularly concerning vehicle’s frontal, side, floor, and roof structures. These distinctions have a substantial impact on crashworthiness performance during frontal, side, and rollover collisions. Most importantly, the use of lithium-ion batteries in BEV is associated with a potential fire hazard [4]. Battery fire accident can be caused by various factors, including mechanical damage [5] , electric malfunction [6,7], expose to high heat or external fire [8,9], involvement in a crash [10], post-crash fire [11,12,13] or other unknown reasons [6].

Given that battery fires are a significant failure mode in electric vehicles, a simulation method capable of modeling battery thermal runaway is essential for the crash safety design and analysis of electric vehicles and their batteries. This method should address critical questions: whether the battery catches fire in specific crash scenarios, the root cause of any battery fires, and the vulnerable parts in the battery and vehicle design during thermal runaway events. Moreover, the method must be computationally efficient for practical industrial application.

A. Lithium-ion Battery Thermal Runaway Process

The lithium-ion battery catches fire in a mechanism called thermal-runaway [14]. Thermal runaway occurs in situations where an increase in battery temperature changes the battery conditions in a way that leads to further increase in temperature, often leading to fire or explosion. To fully understand the physics that leads to thermal runaway, the sequence of how the battery fire started should be understood first. The triggering event of the thermal runaway can be excessive current, overcharging, subjected to elevated temperatures, manufacturing defect, mechanical damage, or their combined effect (cf. Figure 1). High environmental temperature can directly heat up the battery. Overcharging and excessive current can increase the battery temperature through Joule heating. Mechanical damage or manufacturing defects can create internal short-circuits, further raising the temperature through localized Joule heating. Once the temperature reaches 80 degrees Celsius, the heat can cause permanent damage to the anode material. At 110 degrees Celsius, the organic electrode begins to break down, releasing heat and flammable gas. The polyethylene separator film melts at 135 degrees Celsius. After the separator fails, the direct contact between the cathode and anode material forms an internal short-circuit, releasing a large amount of heat and flammable gas, which often leads to fire or explosion [15]. By the time temperature reaches 200 degrees Celsius, the cathode will breakdown and release flammable gas. At this stage, a fire and explosion is likely to occur [16]. Battery thermal runaway typically starts in an individual cell and then propagates to neighboring cells and eventually the entire battery pack. Separator failure is a critical event in thermal runaway. If the separator remains intact, the temperature can decrease after the external circuit is cutoff or colling is applied. However, once the separator fails and the internal short-circuit forms, a large amount of heat is released, causing a larger area of the separator to melt. This, in turn, leads to more internal short-circuits and the release of more heat. This self-sustained chain reaction, known as the thermal runaway, will leads to smoke, fire, and explosion of the battery [17]. The thermal runaway point refers to the starting point of uncontrolled temperature escalation, thus it should be considered the failure point for battery safety analysis [15]. Separator failure most likely triggers a thermal runaway. Therefore, the point of separator failure should be recognized as the critical failure point for battery safety analysis. The purpose of battery simulation is to capture the phenomenon involving the coupling of multiple physics that leads to the starting point of battery thermal runaway. This allows us to predict whether a BEV will catch fire under specific crash scenarios. Phenomena occurring after the initiation of thermal runaway, such as the ignition and combustion processes, are beyond the scope of this study.

B. Selection of Battery Simulation Techniques

1) Finte Element Simulation Methods used in battery simulation are summarized here: The classical method for battery simulation uses mechanical-only finite element analysis, [18,19]; While computationally efficient, this method cannot predict thermal-runaway and battery fire, which are primary concerns for battery safety. It is noted that when a battery cell is mechanically damaged, internal short-circuit and thermal runaway could happens before the cell loses its mechanical strength [20]. In simpler terms, the battery may catch fire before it breaks. Mechanical-only simulation can only predict the structure failure point, leading to the impression that the battery is safer than its actual safety level. This could be dangerous because it may lead to the implementation of inadequate safety measures or the use of the battery in situations where it is not actually safe. Another battery simulation method involves electrical-chemical-thermal-mechanical finite element analysis. Solving electrochemical equations with heat generation significantly impacts the computational efficiency [21], making it impractical to model the entire BEV together with its battery cells in full vehicle crash simulations. Another alternative is electromagnetic-thermal-mechanicalfinite element analysis [22], which involves solving Maxwell’s equation. This method is computationally costly due to the complexity of Maxwell’s equations and their interactions with thermal and mechanical phenomena. Since the widely adopted Randles Circuit Model for battery equivalent circuit [23,24,25] uses only capacitors and resistors; the involvement of electromagnetism is unnecessary. The resistant heating solver is much more efficient for battery simulation. Thermal Electrical Mechanical Coupled Finite Analysis is previously used in lithium-ion battery simulation [26,27,28,29,30,31,32]. In this study, thermal-mechanical-electricalfinite element analysis features the use of thermal solver, electric resistant heating solver, and mechanical solver in a fully coupled LS-DYNA simulation. A user-defined voltage source was utilized to model the capacitor, so the Randles circuit can be modeled with the resistant heating solver [26]. This method to model Randle circuit is previously used to predict the thermal runaway of a single lithium-ion battery cell under normal and abusive conditions [33] [26]. This study demonstrates that the thermal-mechanical-electrical coupled simulation is suitable for modeling the crash safety of lithium-ion batteries integrated with BEVs.

2) Thermal Electrical Mechanical Coupling Formulation in LS-DYNA: The thermal-electrical-mechanical coupling formulation in LS-DYNA is summarized in literature [34,35,36,37]. The Multiphysics solvers are required to concurrently solve Maxwell’s equations, Joule’s Law, and Newton’s Law.

Electrical: the Maxwell equations are approximated using the Eddy-current approximation for conductors with low frequency varying field. It is assumed that $\epsilon \frac{\partial \mathbf{E}}{\partial t}\ll \sigma \mathbf{E}$, where $\mathbf{E}$ is the electric field. $\epsilon$ is the permittivity of free space. $\sigma$ is the electric conductivity. This approximation implies a divergence free current density and no free charge accumulation. The Maxwell equations become:

Faraday Law:

$$\nabla \times \mathbf{E}=-\frac{\partial \mathit{B}}{\partial t}$$

(1)

Amper Law:

$$\nabla \times \frac{\mathbf{B}}{\mu }=\mathbf{j}+\epsilon \frac{\partial \mathit{E}}{\partial t}$$

(2)

Gauss Law:

$$\nabla \bullet \mathbf{B}=0$$

(3)

$$\nabla \bullet \epsilon \mathbf{E}=0$$

(4)

$$\nabla \bullet \mathit{j}=0$$

(5)

Ohm Law:

$$\mathbf{j}=\sigma \mathbf{E}+{\mathbf{j}}_s$$

(6)

where $\mathbf{B}$ is Magnetic flux density, $\mathbf{j}$ is Current density, ${\mathbf{j}}_s$ is a divergence free source current density，$\mu$ is permeability, $\epsilon$ is permittivity. The divergence condition (3) allows $\mathbf{B}$ writing as:

$$\mathbf{B}=\nabla \times \mathbf{A}$$

(7)

where $\mathbf{A}$ is the magnetic vector potential [38]. According to Equation (1), the electric field is:

$$\mathbf{E}=-\nabla \varphi -\frac{\partial \mathbf{A}}{\partial t}$$

(8)

where $\varphi$ is the electric scalar potential. The Gauge condition is used as:

$$\nabla \bullet \sigma \mathbf{A}=0$$

(9)

which allows a separate the vector potential from the scalar potential in the equations, implying Equations (5), (6), (8), and (9)

$$\nabla \bullet \sigma \nabla \varphi =0$$

(10)

which implies Equations (2), (6), (7), and (8)

$$\sigma \frac{\partial \mathbf{A}}{\partial t}+\nabla \times \frac{1}{\mu }\nabla \times \mathbf{A}+\sigma \nabla \varphi ={\mathbf{j}}_s$$

(11)

Equation (10), (11) is Maxwell equations in term of the 2 potentials. Let $\mathsf{\Omega }$ be a set of multiply connected conducting regions. The surrounding insulator exterior regions is ${\mathsf{\Omega }}_e$. The boundary between $\mathsf{\Omega }$ and ${\mathsf{\Omega }}_e$ is $\Gamma$, and the artificial boundary on $\mathsf{\Omega }$ at the end of the meshing region, where the conductors are connected to an external circuit, is called ${\Gamma }_c$. In the following, let $\mathbf{n}$ be the outward normal to surfaces $\Gamma$ or ${\Gamma }_c$. In ${\mathsf{\Omega }}_e$,we have $\sigma =0$ and $\mu ={\mu }_0$.

$$\mathbf{n}\bullet \nabla \varphi =0on\Gamma$$

(12)

$$\varphi ={\varphi }_{c}on{\Gamma }_c$$

(13)

$$\mathbf{n}\times \nabla \times \mathbf{A}=A_{e}on\Gamma$$

(14)

$$\mathbf{n}\times \mathbf{A}={\mathbf{A}}_{c}on{\Gamma }_c$$

(15)

Equation (13) and (15) allows the connection of the conductors to a voltage source and a current source, respectively. Once the potentials are computed, the electromagnetic fields are given by Equations and (6), (7), and (8).

Mechanical：After the EM fields are computed by solving Maxwell equations, the Lorentz force is evaluated at the nodes and added to the mechanical equation.

$$\rho \frac{D\mathit{u}}{Dt}=\nabla \bullet \mathbf{s}+\mathbf{f}+\mathbf{j}\times \mathbf{B}$$

(16)

where $\mathit{u}$ is the displacement. $\mathbf{s}$ is stress. $\mathbf{f}$ is the external force. $\rho$ is mass density. Notice the extra Lorentz force term $\mathbf{j}\times \mathbf{B}$ represents the electric mechanical coupling.

Thermal: The thermal equation

$$\frac{De}{Dt}=\mathbf{s}:{\dot{\mathsf{\epsilon }}}^{pl}+\frac{Dq}{Dt}+\frac{j^2}{\sigma {\rho }^e}$$

(17)

where $e$ is the internal energy. $\mathbf{s}$ is Stress tensor. ${\dot{\mathsf{\epsilon }}}^{\mathit{p}\mathit{l}}$ is Plastic strain rate tensor. $q$ is the external heat source. ${\rho }^e$ is electrical resistivity. The Joule heating term $j^2/\sigma {\rho }^e$is added to the thermal solver allowing to update the temperature.

Equation of State (EOS): The material’s mechanical properties are evaluated at Gauss point temperature. The conductivity also may vary with temperature change.

$$\sigma =\sigma (T,\rho )$$

(18)

where $T$ is temperature. $\rho$ is the mass density.

3) Geometry Modeling of Battery Cell. The lithium-ion battery consists of numerous stacked layers of very thin membranes [39]. For porch cells, mechanical behavior in the thickness direction is particularly important for determining the strength of the cell in hemispherical punch test, radial bending test, and side crushing test [39,40]. Utilizing solid elements is crucial for realistically representing the strength of the battery cell in the through thickness direction. However, the small membrane thickness (0.01mm~0.06mm) [41,42] requires a very small element size and timestep in simulation, while the presence of a large number of stacked layers (>50 layers) leads to a very large model [43,44]. These factors make the simulation computationally impractical if the cell is modeled with realistic geometry. In this study, the battery cell is modeled with a representative battery model: Each cell component, including cathode current collector, cathode, separator, anode, anode current collector, is modeled separately with a representation geometry. Since the lithium-ion battery features a jelly roll structure of many thin layers. The thickness of the layers is proportionally scaled up, so a handful of thick layers can be used to represent many thin layers. Realistic material properties of each cell component, including density, Young’s module, Poisson’s ratio, stress-strain relation, porosity, thermal conductivity, specific heat capacity, electric resistivity, etc., were used in the FE model. The alternative method is the homogenized battery model: the entire battery cell is treated as one piece of material. The phenomenon-based material model is fitted by matching cell-level mechanical tests [19]. This material model represents the overall deformation pattern of the battery cell, but it does not represent any realistic material. The representative model is closer to true physics than the homogenized model because it incorporates the realistic mechanical, electrical, and thermal properties of each battery component. The battery cell comprises a loosely packed jelly roll structure with materials such as metal, polymer, graphite, each having very different material properties such as Young’s modulus, stress-strain curve, thermal conductivity, resistivity, melting point, and more. It is not feasible to model the mechanical, electrical, and thermal behavior of each cell component accurately if they are bundled together and treated as a single material. For instance, the melting point for each cell component is different. The polymer separator fails around 135 °C, some ceramic coated separators can withstand 200°C. The aluminum current collector sheet fails at 660 °C. Modeling these components together as a single material in the homogenized approach would overlook such critical distinction. Moreover, the phenomenon-based material model in the homogenized model lacks physical representation. The material model fitted to match one test may not match another test conducted on the same cell under different testing conditions. The homogenized model also struggles with electric modeling because it cannot accurately represent the internal battery circuit and resistivity variations among different cell parts. In contrast, the representative model can realistically predict internal resistance, voltage distribution, current flow, Ohmic heating power, and temperature distribution for each layer of the battery. Previous studies using a multi-physics single-cell representative model have shown its capability to match mechanical and electrical test results [26]. Following a similar methodology, the model employed in this study is calibrated using external short circuit tests and spherical punch tests. It can simulate the coupled physical phenomena that lead to battery thermal runaway.

C. Full Vehicle Finite Element Models for Crash Analysis

Finite element analysis has been widely used to study the crashworthiness of automobile and roadside hardware. [45,46,47,48]. The National Highway Traffic Safety Administration (NHTSA) has sponsored the development of public domain finite element models for more than two decades to promote automobile safety and support the rule making effort [49]. Center for Collision Safety and Analysis (CCSA) at George Mason university has developed or managed over 20 public domain vehicle crash analysis models based on the reverse engineering modeling protocol [50]. (cf. Figure 2) These finite element models are created by reverse engineering the vehicle on the market, and they are also calibrated with the full vehicle crash tests to ensure physical representation. These simulation models provide a powerful tool for research in vehicle crash safety, occupant safety, pedestrian safety, restraint system design, biomechanics, fleet compatibility, roadside hardware, accident reconstruction, and many other field. These models also serve as a modeling guideline for the state of the art of vehicle crash analysis. The typical procedure for developing a reverse-engineered crash analysis model is summarized as follows: The vehicle is tested for its weight, center of mass, and moment of inertia before disassembly. At the time of disassembly, the geometry and position of each part is laser scanned, the thickness and mass of each part are also measured. The material test coupons are cut from parts to determine the stress-strain relation for the material model. Finite element mesh is generated based on the part geometry, thickness, and mass following a strict quality control criterion. The mesh of different parts is then connected with spotwelds, constraint rigid body, beam element, rotatable joints and tied contact. The suspension, added mass, airbags and restraint system are defined next. The contact is then defined with penetration checked. The finite element model is verified to match mass, center of gravity, and moment of inertia of the vehicle. The developed finite element model is validated against multiple crash tests, including NCAP Frontal Impact, IIHS Offset Deformable Barrier Test, Pole Impact Test, Roof Crash Test, etc., The modeling is seen to provide satisfactory result [47,51,52]. The models feature accurate geometry, realistic material property, a soft interior, deployable airbags, and full restraint system, making it a close representation of the actual vehicle.

The classical crash safety simulation models are built in mechanical-only analysis for the ICE vehicle. The reverse engineering FE model for battery electric vehicles requires a proof-of-concept study, especially for the method to model battery. There are two decisions to be made regarding the technological path, with computational efficiency being the primary concern: 1. Electrical-Thermal-Mechanical Analysis vs. Mechanical-Only Analysis: it is seen that thermal-electrical-mechanical coupled analysis is necessary for capturing the battery thermal runaway. However, it was uncertain whether it was practical to simulate the entire BEV with all active connected battery cells in a multi-physics analysis. If this analysis takes too long to run on a typical cluster computer, then the mechanical-only model would be the only option. 2. Battery Homogenized model vs. Battery Representation Model: The battery homogenized model has the entire cell modeled with one chunk of material. The representation model has each cell parts modeled separately, so the physical material property can be used. The representation model is more accurate but can be more expensive in computational time. It was uncertain whether the representation battery model runs fast enough for a full EV simulation on a typical cluster computer. If the representation model takes too long to run. Then the homogenized model would be the only option.

This proof-of-concept study has demonstrated that the thermal-electrical-mechanical coupled analysis with representation battery model is able to simulate the crash safety for the entire EV along its battery cells. A typical multi-physics EV crash simulation can be completed within 24 hours using a cluster computer with 24 CPUs.

2. Methodology and Result

The thermal-electrical-mechanical coupled crash analysis model is built in three stages: 1. A single battery cell; 2. A multi-cell battery module; 3. A generic battery electric vehicle with 24 active cells. The exterior and interior features of the generic BEV model represent the first-generation Google Waymo BEV. The body in white referenced the Toyota Yaris crash analysis model [53]. The battery cell and module represent the system used in first-generation Nissan Leaf [54]. The battery cell modeling methodology is developed in previous study [26].

A.

Single Cell Model

The geometry of the single cell model represents the battery cell used by Nissan Leaf [54]. A repurposed Nissan leaf battery module (first generation) with 4 cells has been disassembled to the cell level. The model closely represents the exterior geometry of the cell. The cell dimension is $261\mathrm{m}\mathrm{m}\times 216\mathrm{m}\mathrm{m}\times 79\mathrm{m}\mathrm{m}$excludes the plastic wrapper, which has a thickness of $0.14mm$. The parts dimension inside battery cell is proportionally scaled from the cells in previous studies [26,43]. Five plates are used to represent the multi-layer sandwich structure, each representing one battery part (cf. Figure 3). In simulation, the cathode current collector has thickness of 0.857mm, the cathode layer has thickness of 1.745mm, the gap between cathode and separator is 0.266mm, the separator has a thickness of 0.932mm, the gap between separator and anode is 0.266mm, the anode layer has a thickness of 3.054mm, the anode current collector has a thickness of 0.779mm [43]. The thermal, electrical, and mechanical properties of each cell materials are summarized in previous study [26]. Two layers of solid elements are used on each plate to capture the strength of the cell in the thickness direction. The cathode and cathode current collector plate are connected by sharing nodes; and the same also applies to the anode and anode current collector. A gap is necessary to represent the mechanical strength of the loosely packed batteries layers [26]. The plate representing separator is modeled with some “wrinkles” to reflect the loosely packed sandwich structure (cf. Figure 3). Only certain parts of the separator touch the anode and cathode through frictional contact, without sharing nodes. This configuration creates a gap between the separator and the two electrodes while still maintaining structural integrity under gravity

The electric potential difference is defined between Node Set 1 (red) and Node Set 2 (blue) as shown in Figure 3. Only two node sets are defined for each battery cell. Node Set 1 contains all nodes with positive voltage, while Node Set 2 contains all nodes with zero voltage. Since the voltage difference on the same battery component is usually very small, it is more efficient to define only one Randles circuit for each battery cell [23].

The electric current starts from Node Set 1 on the side of the cathode layer closer to the separator. The current flows through the cathode, cathode current collector, and cathode current collector tab to the external circuit. The circuit returns from the anode current collector tab, through the anode current collector and anode, and finally back to Node Set 2, which is on the other side of the anode layer closer to the separator. The capacitor in the Randles circuit is analytically modeled by a generalized voltage source (Figure 4), enabling the use of a resistance heating solver to simulate the Randles circuit. The detailed methodology and source code for the generalized voltage source are presented in previous work [26].

The external short-circuit test (cf. Figure 4) is simulated using an FE battery cell model in LS-DYNA. The cell undergo normal discharging until an external short circuit is formed at $5\times {10}^{-4}s$. The predicted current from the LS-DYNA simulation is compared to the result from CircuitLab, an electric circuit simulation software [55]. It is seen that the FE battery cell model accurately captures the dynamic circuit behavior and predicts the temperature rise of the separator. The FE simulation also precisely captures the peak current and the steady state current.

The mechanical property of the battery cell model is verified with spherical punch test. [56] A 1-inch diameter sphere indenter is used in the punch test. The LS-DYNA simulation result is compared with the test data in Figure 5. The results show that the battery cell FE model matches the mechanical property in general. It is noted that the porosity of the lithium material and graphite coating will influence the stiffness of the cell under deformation [26]. The Young’s material and plastic strain vs stress curve is scaled to 20% of the standard book value to reflect the porosity of the coated layers. Cell disassembly and detailed reverse engineering of the component part is necessary to further improve accuracy of the cell model.

Thermal contact within the battery cell is also incorporated into the model, as it is a cortical factor in predicting physical phenomenon that can lead to thermal runaway. The nonlinear transient thermal solver in LS-DYNA is used. A hybrid thermal contact model is used, considering the heat conduction, convection, and radiation. [57] The heat conductance between contact surfaces is calculated as

$$h=\left\{\begin{array}{c}{h}_0, 0\le l_{gap}\le l_{min}\\ h_{cond}+h_{rad}, l_{min}<l_{gap}\le l_{max}\\ 0, l_{gap}>l_{max}\end{array}\right.,$$

$$h_{cond}=\frac{K}{l_{gap}}$$

19

$$h_{rad}=f_{rad}(T_1+T_2)(T_1^2+T_2^2)$$

(21)

$$f_{rad}=\frac{\sigma }{\frac{1}{{\epsilon }_1}+\frac{1}{{\epsilon }_2}-1}$$

$l_{gap}$ is the gap between the two contact surfaces. $h_0$ is direct contact heat transfer conductance between two surfaces with a gap smaller than $l_{min}$. $h_{cond}$ is the heat conductance of the fluid filling the gap fluid between the contact surfaces $K$ is fluid’s thermal conductivity. $K=0.58W/m$ and $K=0.024W/m$ were adopted for thermal conductivity for the liquid electrolyte and air, respectively [58,59]. $f_{rad}$ is the radiant heat transfer factor between contact surfaces. Stefan-Boltzmann constant $\sigma =5.6703\times {10}^{-8}(W/m^2{K}^4)$. ${\epsilon }_1,{\epsilon }_2$ are the emissivity of two surfaces. The emissivity of black body ${\epsilon }_1={\epsilon }_2=1$ is assumed [57]. $T_1,T_2$ are the absolute temperature of two contact surfaces, respectively. Inside the wrapper of the battery cell, the thermal conductivity of water ($K=0.58W/m$) is used to represent the liquid electrolyte between the battery layers. Outside the battery cell, the thermal conductivity of air $K=0.024W/m$ is used for all other situations. Two different contact cards are set up separately to define the thermal contact within the battery cells and outside the battery cells. The same contact card also defines the mechanical contact between surfaces (cf. Figure 6,Figure 7) The surface-to-surface contact needs to be used for the combined thermal and mechanical contact. The dynamic friction coefficient, static fiction coefficient, and exponential decay coefficient are set to be $0.005,0.004$, and $0.001$, respectively [57].

B.

Battery Module Models

A Battery module consists of multiple individual battery cells assembled into a single unit. These cells are connected in serial or parallel to achieve the desired voltage and capacity. The module case offers structure integrity, housing the electric connection and protects the cells against external impact. Battery module is often the minimal replaceable part in an EV. The FE model of the battery module represents the system used in Nissan Leaf [54]. The battery module contains four lithium-ion battery cells in a 2s2p configuration, where two pairs of parallel-connected cells are connected in serial. The module exterior dimension is 303$\times$223$\times$55 mm, and the weight is 3.8kg [60]. In the simulation, the element size of the battery layers is 2mm in plain and 0.34~1mm in thickness. There are a total of 275,577 elements in this model. The cathode, cathode current collector, separator, anode, anode current collector, the module case and case cover are modeled with solid element (cf. Figure 8). The cushion and plastic wrap are modeled with shell element. The dimensions of the single cell model from the previous section has been slightly adjusted to fit into the case geometry of Nissan Leaf battery module. The material parameters for the battery module model are the same as those used in the single cell model [26]. Simulation of an event for 0.2 seconds in the battery module model takes approximately 3 hours. The electrical functions in this model operate as expected.”.

1)

Electric Verification

It is necessary to verify whether the use of user-defined voltage source to model the Randles circuit can be extended to multiple cells in electrical-thermal-mechanical simulation. Verification cases include two different battery cells connected in serial and parallel (Figure 9, Figure 10). Noted that the voltage and the capacitance for two battery cells are different. The result demonstrates that the electric simulation functions correctly for the two cells in live circuit connections. The method to model the Randles circuit with a generalized voltage source is described in the appendix of previous work [26]. It is important to ensure that each cell has a unique function name and ID for the generalized voltage source that are connected in circuit. Relative potentials must also be used in the user-defined voltage source, with only one ground defined with an absolute potential value for the entire circuit. The thermotical values from CircuitLab [61] are used to verify the electrical simulation in LS-DYNA. The results show that the dynamic effect of the circuit can be accurately captured in LS-DYNA’s multi-physics simulation. Therefore, the electric modeling method used for single cells can be extended to model multiple cells connected in circuit.

The simulation of an external short-circuit is used to verify the electrical behavior of the Leaf Battery module, which consists of 4 cells as shown in Figure 11. The result aligns well between the CircuitLab and LS-DYNA simulations. Notice the internal resistance of the cells (cf. R1, R2, R3, R4 in Figure 11) is determined through finite element analysis result and serves as input in the CircuitLab simulation. The LS-DYNA simulation calculates the internal resistant for each cell as 2.927$m\mathrm{\Omega }$. Table 1 compares the internal resistant value from the finite element simulation with the value from battery module testing and literature sources. The battery module is tested for its DC internal resistant and AC impedance at 1kHz using an electrochemical workstation CHI660E. Additionally, AC impedance is tested with battery impedance tester YR1035+. The finite element simulation result of 2.927$m\mathrm{\Omega }$ closely match the literature value ranging from 1.1 $m\mathrm{\Omega }$ to 6.0 $m\mathrm{\Omega }$ in literature [62]. The measured AC impedance and DC internal resistance are higher than the literature value, possibly due to the fact that the test subject is repurposed electric vehicle battery module retired from EVs.

The current density contour shows the expected higher current accumulation near the tab area, indicating that the electrical stimulation works effectively in the battery module model with 4 cells connected in a live circuit.

2)

Mechanical Verification

We simulated the spherical indentation test conducted on the Leaf module by Oak Ridge National Laboratory, using a 2-inch diameter sphere as the punch head at a punching speed of 0.1 inch/min [20,56,63]. The result shows a general match in the force vs. displacement response. It should be noted that the module case serves as the primary load -bearing structure in the battery module design. A detailed reverse engineering process and material modeling of the battery module case and cushion pads are necessary to further increase the simulation accuracy.

Figure 13. Experiment setup [20,63](Top Left), Simulation Model(top Middle) and Force vs. Displacement(Bottom) in Spherical indentation.

Figure 13. Experiment setup [20,63](Top Left), Simulation Model(top Middle) and Force vs. Displacement(Bottom) in Spherical indentation.

3)

Thermal Verification

The resistive heat generation, heat conduction and thermal contact of the battery module model are verified in an external short-circuit simulation. The heat generated by Joule heating gradually propagates throughout the entire module (cf. Figure 14). The primary heating source is observed in the cathode layer. The area near the current collector tab has higher temperature due to the higher current passing through that location. The heat transfers from the cathode to the metal cathode current collector, then through the plastic wrapper, the cell cushion pad and finally to the aluminum module cover.

4)

Thermal-Mechanical- Electrical Verification

The thermal-mechanical-electrical coupling is verified in a side impact simulation. In the Nissan Leaf, battery modules are installed horizontally [64]. Besides Spherical indentation, side impacts are common collision directions. Battery module may deform similarity in frontal or side impact. The pole impact simulation is shown in Figure 15, with the pole having a radius of 50 mm and an impact speed of 1 m/s. The temperature increase in the impact area is primarily attributed to the mechanical work caused by plastic deformation.

5)

Coarse Battery Module Model

A coarse battery module model containing only one cell was also developed for use in full electric vehicle (EV) analyses. The goal is to reduce the CPU time to 24 hours for a complete EV crash simulation with active battery cells in a multi-physics analysis. (cf. Figure 16). This coarse model uses a single 7.4V cell to represent the battery module circuit shown in Figure 11.

A side pole impact was also conducted for the coarse module mode to verify the functionality of thermal-electrical-mechanical coupling. The temperature and Ohm heating power demonstrated that the battery module can accurately simulate heat generation during impact. (cf. Figure 17). Besides the heat generation though work hardening process. It was observed that cell deformation influence the internal resistance, leading to localized current increase, thereby locally raising Ohm heating power and temperature.

C.

Battery Electric Vehicle with 24 Active Cells.

A generic battery electric vehicle with 24 active battery modules has been developed. Body-in-White is based on Toyota Yaris, the exterior and interior features the first-generation Google Waymo battery electric vehicle. The battery cells and modules follow the same configuration as the Nissan Leaf. This generic battery electric vehicle model served as a proof-of-concept study for future EV crashworthiness models. The battery pack consist of 24 active battery cells connected in serial. The battery modules are packed and integrated into vehicle floor structure (cf. Figure 18)

1)

BEV in Normal Usage & External Short-Circuit

The case of BEV with a normal discharging battery pack is simulated. An external resistor of $1\mathrm{\Omega }$ is connected to the battery pack to represent battery discharge under normal operating condition (cf. Figure 20). Recall that the BEV model consists of 24 battery modules in serial connection, with each module carrying a 7.4V Voltage. The simulation contour of electric potential confirms the serial connection of cells. The cell on upper left is connected to the ground and the cell on the lower right has the highest voltage of 177.6V (cf. Figure 19 left). All battery cells exhibit similar current density contours due to their serial connection, with higher densities observed near the current collector tab (cf. Figure 19 Right). Both the electric potential and current density contour suggest that the electrical simulation operates correctly. The temperature contour during normal discharged is shown in Figure 20.

For comparison, the case of external short-circuit is also simulated for the battery pack installed in the BEV. Due to the short-circuit, the voltage contour of the battery is almost zero (cf. Figure 21, Top Left). The current is over 30 times higher than the normal levels (cf. Figure 21, Top Right). These voltage and current results confirm the proper functionality of the electrical stimulation. The simulation results show a rapid temperature increase on the cathode, cathode current collector, and separator (cf. Figure 21, middle). The cathode is the primary heat source due to the high resistivity of cathode material. The heat generated at the cathode transfers to the adjacent cathode current collector and separator. The simulation shows the separators reach the melting point of 150$\mathrm{℃}$ and failed at t=0.2s, suggesting that an external short-circuit of the battery pack could lead to nearly instantaneous battery thermal runaway (cf. Figure 22).

2)

BEV in FMVSS Pole Impact (20mph impact speed):

The side pole impact could potentially lead to intrusion on the battery compartment. The full vehicle pole impact simulation model (cf. Figure 23) was developed in accordance with Federal Motive Vehicle Safety Standard (FMVSS) No. 201 and 214 [65]. The rigid pole has a radius of 127.5 mm and an initial velocity of 20 mph upon impact. The pole impact simulation required 24 hours and 40 mins to simulate the event within 0.2s. The battery system undergoes normal discharge, representing a typical driving condition. The thermal and electric simulation performed as expected in Figure 24. It is seen that the heat generated by internal resistance is concentrated near the current collector tab on the cathode layer of the cells. The battery has a 20°C temperature rise due to normal discharge. (cf. Figure 26). In normal usage of BEV, the heat will be carried away by the battery active cooling system. However, the active colling system will most likely being shut down by battery management system during crash. For simplicity, this simple conceptual model does not model the active cooling system for the BEV. It is observed that the cells near the impact area exhibit a 4°C higher temperatures compared to the cell farther away (cf. Figure 26), which are largely uninfluenced by the impact and can be considered as the benchmark of normal discharge. Notice the time for cell plastic deformation is from 0.01s to 0.11s as shown in Figure 26 (Top). the abnormal temperature rise on the cathode layer is also from 0.01s to 0.11s as show in Figure 26(Bottom) Therefore, the temperature rise is primarily caused by plastic deformation. Notice the 4°C additional temperature rise is minimal (cf. Figure 24). Therefore, the simulation does not predict an immediate risk of thermal runaway. The cathode layer shows a 1.5% of plastic deformation (cf. Figure 26). The non-uniformity caused by permanent change of battery geometry is relatively small, so it is estimated the cell internal resistance will not be significantly influenced. The post-crash cell is still safe in future charging and discharging cycles. The plastic porch bag is seen to have a 70% plastic deformation (cf. Figure 28). If we assume the plastic cell wrapper used to make the porch bag fails at 100% of strain. The cell wrapper will remain intact after impact. The post-crash battery also has no thermal runaway risk caused by cell leakage.

2)

BEV in Pole Impact with 25mph impact speed:

Higher Pole impact speeds of 25 mph were examined to assess the thermal runaway risk. A moderate deformation is seen for the battery module closer to the impact area (cf. Figure 28 Middle. The maximum temperature contours for the 25mph case are presented in Figure 25. It is seen the cells near the impact area exhibit a 12°C temperatures raise compared to the cell far away from impact (cf. Figure 26). Notice the time for cell plastic deformation is from 0.04s to 0.08s. However, the significant temperature rise on the cathode layer is from 0.08s to 0.1s as show in Figure 26(Left). This temperature rise is primarily caused by the additional joule heating caused by cell internal resistance change, which is in turn caused by cell plastic deformation. Plastic work is the secondary effect that causes the cell temperature raise. Notie the cell far away from the impact area has a maximum temperature of 313K, which is identical to the 20mph case. (cf. Figure 26) This verified the assumption that the cell far away from impact area is largely uninfluenced, therefore can be treated as normal discharging cells for comparison purposes. This temperature rise of 12°C is not enough to cause a thermal runaway (cf. Figure 24). Therefore, the simulation does not predict an immediate risk of thermal runaway for pole impact at 25mph speed. However, the cathode layer shows a 15% of plastic deformation (cf. Figure 28). If this damaged cell is still used in normal charging and discharging cycles, the permanently deformed cathode layer will cause uneven change of local resistivity, which may lead to the accumulated heat and a risk of thermal runaway. The plastic cell wrapper is seen to have a 70% plastic deformation (cf. Figure 29). If we assume the plastic cell wrapper fails at 100% of strain. The cell wrapper will remain intact. The cell will not leak.

Figure 25. Temperature Distribution at maximum intrusion moment in Pole Impact at 25mph impact speed.

Figure 25. Temperature Distribution at maximum intrusion moment in Pole Impact at 25mph impact speed.

3)

BEV Pole Impact with 30mph impact speed:

The battery module undergoes severe deformation at this impact speed. The current collector tab is bent, the plastic porch has multiple failure points with plastic strain greater than 150% at the time of 0.04s (cf. Figure 28, Figure 29). The broken porch bag will cause the disintegration of the battery structure and exposure the flammable materials to the open air. The damaged porch bag is identified as the primary failure mode that may lead to thermal runaway under this speed. The severely damaged cell will also distort the finite element mesh to cause a stability problem in simulation. Immediate thermal runaway is predicted for BEV Pole impact with 30mph impact speed. Table 2 summarizes the simulation thermal runaway prediction for all 3 impact speeds.

5. Conclusions

The thermal-electrical-mechanical coupling FE Model with representative geometry is developed to model the battery cell, battery module, battery pack and battery electric vehicle. The result proves the feasibility of modeling the EV and battery cells in a multi-physics simulation for thermal runaway prediction. The multi-physics simulation captures various mechanical, electrical, thermal, phenomena and their coupling effects, so the thermal runaway and battery fire can be predicted. The simulation methodology is efficient, expecting a 24-hour turnaround for a typical EV crash simulation involving active cells. Future studies should emphasize utilizing more realistic material models in the simulation. The material model should come from material test of the battery component material, including electrical, mechanical, and thermal tests. Reverse engineering with laser-scanning should be used in the modeling of battery cells, battery modules, battery packs, and BEV.