A Simulation Framework for Battery Cell Impact Safety Modeling Using LS-DYNA

The mechanical model furnishes the solution to the momentum equation. The individual components, namely the positive current collector, negative current collector, positive active material coating, negative active material coating, and separator are modeled with a single layer of thick shell (TSHELL) elements. The adjacent layers share common nodes, since it is assumed the exterior pouch is completely sealed and exterior atmospheric pressure prevents sliding motion between the layers. Element formulation 5 is selected because it enables 3D stress update with strain reduced integration,¹⁶ giving the TSHELL element solid-like behavior and ability to handle compressive loading.

The separators are thin porous polymer sheets, modeled using Material Model 181: Simplified Rubber/Foam (MAT-181), which is a simplified "quasi"-hyperelastic rubber or foam model defined by a single uniaxial load curve or by a family of curves at discrete strain rates. This material model is selected for its ability to simulate soft polymeric materials. The detailed element description and complete formulae can be found in Ref. 17. The input load curve for the material model is procured from tensile and compressive experiments performed on separators extracted from cell teardowns.

The electrode sheets are thin metal current collector foils coated by porous coating materials on each side. The coatings consist of active material particles which host lithium upon intercalation, carbon conductive additives, and polymeric binder. In the finite element model, all three layers (coating-substrate-coating) are resolved for both the positive and negative electrodes. The metal substrate layers are modeled using Material Model 123: Modified Piecewise Linear Plasticity (MAT-123),¹⁷ a material model that is capable of simulating the non-linear yielding behavior of metal and is expandable to incorporate strain rate effects in future studies. The coating layers use the previously mentioned MAT-181. Literature data for pure copper¹⁸ and aluminum¹⁹ materials are used to populate the inputs of the substrate material models. While the strength of the substrate dominates the overall strength of the coated electrode sheets, the properties of the coating materials control its failure strain as discussed in more detail in the Results and discussion section. Therefore, the failure strain and ultimate stress of the coating materials are used as tuning parameters. The tuning objective is to make the peak force, and elongation at peak force, in simulation match that of the physical test.

The standard LS-DYNA resistive solver²⁰ is capable of determining the potential and current distributions within an inert conductor, such as a current collector, by solving the Laplacian equation subject to arbitrary boundary conditions. To empirically represent the electrochemistry occurring within the active materials and separator of a cell, an equivalent circuit model is coupled to a modified version of the Laplacian equation by dividing each current collector into numerous finite elements and then assuming equivalent circuits govern the relationship between transverse current and local voltage across each unit cell. The relationship between the in-plane and transverse governing equations is illustrated in Figure 2. The governing equations of the generic equivalent circuit are adopted from literature,²¹ where the rationale and implementation of these models is well-documented and therefore not expanded upon here.

Zoom In Zoom Out Reset image size

The electromagnetism resistive solver computes the Laplacian equation for the scalar potential, ϕ,

where σ is the electric conductivity. Equation 1 represents the physical equation, which after projection against the nodal basis functions, φ, using a standard finite element method, reads:

where S is the σ-dependent Laplace operator. Equation 2 is solved along with adequate boundary conditions discussed further in the Case study setup section. In Eq. 2, φ represents the vector of nodal values of the scalar potential, hence φ₁ = φ(N₁) is the value of the scalar potential at node N₁, for example. From the scalar potential, the physical in-plane current density is computed as:

The coupling to the local equivalent circuits creates a non-zero right-hand side in Eq. 2, which represents the potential drop between the nodes of opposing current collectors defined as

where φ₂ and φ₁ are the local potentials at the positive and negative current collector nodes N₂ and N₁ respectively, i is the local transverse current density, r₀ is the local instantaneous resistance, u is the local open-circuit voltage, and v_D is the local diffusion overpotential. All local variables are assigned to an opposing node pair on the positive and negative current collectors, so the total number of local variables of a given type is equal to the number of in-plane, surface nodes on either the positive or negative current collector (38400 in this instance). For numerical stability reasons, the r₀i term of Eq. 4 is included within the Laplacian operator instead of in the right-hand side. The modified system thus reads

where, like S(σ), D(r₀) is a matrix of size (n, n), where n is the number of nodes in the model. The matrix D(r₀) is defined as

for all pairs of nodes (N₂, N₁) associated to the equivalent circuits. Furthermore, b is defined as

Once the system defined by Eq. 5 is solved, the local current in each equivalent circuit can be computed by

which allows actualization of the equivalent circuit governing equations. There are two state variables associated with each local equivalent circuit. The first is the diffusion overpotential, defined as the voltage across the resistance-capacitor (RC) pair given by

where c_D is the equivalent capacitance of the RC pair, and r_D is its resistance. The state of charge (SOC) of the local equivalent circuit is computed by

where c_Q is a conversion factor and Q is the capacity delivered by a unit cell pair. During mechanical deformation, a short circuit may occur based on a threshold criterion selected for the simulation, as discussed further in the Case study setup section. When a short circuit is detected, the equivalent circuit governing equations at a particular node are replaced with pure resistive elements, which act as a constraint for the remaining equivalent circuits. The electrical model parameters r_o, r_D, c_D, are assumed to depend on temperature according to

where β, α, and E_a are fitting constants and R is the universal gas constant.

The thermal solver computes the finite element solution to the transient, three-dimensional heat diffusion equation

where T is the temperature, ρ is the mass density, c_p is the specific heat capacity, k is the thermal conductivity, and represents all volumetric heat source terms. The volumetric heat source terms considered within this model are

where represents the Joule heating that results from in-plane electron transport within the current collectors, i²r₀ is irreversible heat generation associated with ohmic and kinetic losses across the unit cell, and is the combined reversible heat generation associated with the electrode couple. These heat source terms are applied at the node existing at the interface between the active material coating and current collector, and diffuse into the remaining components accordingly. They are also corrected for the volume of the element, V_e, to maintain consistent units with Eq. 12, whereas the in-plane Joule heating is already consistent because it is calculated directly from the mesh. Note that the diffusion overpotential is not included as an irreversible heat generation term since heat of mixing effects are assumed to be negligible.²² A single, isotropic thermal constitutive model is implemented with parameters from literature²³ for battery components soaked with electrolyte. If a short circuit is initiated, the heat source terms are replaced within the short circuited region by

where r_S is the resistance of the local short circuit. At the time of the short, the total energy remaining in the corresponding equivalent circuit is added locally to the thermal solver. As presented, the model does not conserve energy since the heat produced by the diffusion overpotential is not included in the heat diffusion equation. This would produce erroneous results for situations involving normal cycling and has been corrected in future releases of the model. However, data is presented in the Results and discussion that indicates the heat generated by the diffusion overpotential is less than 1% of the magnitude of Eq. 14 for the cases presented in this paper.

The case studies analyzed here assume that rapid mechanical deformation is imposed on the cell on a time scale of approximately 0.5 ms by an indenter with fixed velocity. It is further assumed that all kinetic energy is converted into mechanical strain energy of the cell components within that time span, and the indenter instantaneously stops at 0.5 ms. Within the initial 0.5 ms where mechanical deformation occurs, explicit time integration is used to solve the structural equations and implicit time integration is used to solve the electrochemical and thermal equations. All sets of physics in this initial phase use a time step of 2.0 × 10⁻⁷ s. Following this initial phase, the mechanical deformation is assumed to be rigid and fixed for the remaining simulation time. Therefore, the time step of the implicit time integration is increased to 1 s without any stability issues, and no further calculations are performed for the structural equations.

In two separate cases, cylindrical and spherical indenters of 75 mm diameter are implemented in simulation as rigid bodies using ^*RIGIDWALL_GEOMETRIC. These cases create different failures within the simulated cell to investigate the parametric dependence of the total short circuit resistance that results from varying types of mechanical deformation. The failure criterion for the onset of short circuit is a compressive strain threshold applied at the unit cell level. When the distance between the current collectors is reduced by ten percent of its initial value the local equivalent circuit is replaced by a purely resistive constraint between the opposing current collectors. The resistivity value chosen in this study is 12.5 mΩm² which leads to a short circuit that discharges at a rate that is slow enough to observe the thermal and electrical dynamics. This failure condition is arbitrary but useful for demonstrating the effectiveness of the software tool, and it is expected that future users will incorporate their own comprehensive failure conditions that are representative of the system they are studying. For example, data representing contact resistance between various pairs of components could be leveraged.

The boundary conditions imposed on the cell encompass mechanical, thermal, and electrical constraints. It is assumed that no external electrical load is placed on the cell tabs during the crush, so that current exiting the cell tabs is zero. The negative terminal potential is also set to zero to consider only the potential difference between the two terminals. It is further assumed that the cell exchanges heat with its ambient environment by natural convection from its top surface, with a convection coefficient of 10 W/m²/K, and all other external surfaces are adiabatic. A rigid wall supports the cell along its bottom surface and prevents displacement in the x,y, and z-directions, as well as rotation about those nodes.

The stress-strain curves for the tensile test of each component are plotted in Figure 3. In the case of the separator, the material stress-strain data is used directly as a model input. The difference between stress-strain curves for different orientations is not significant compared to sample variation for this particular separator type, so the simulation input is based on the average values of data from all three orientations. In the case of the electrodes, the test results show that, on average, the positive electrode specimen ruptures under a stress of approximately 17 MPa at 0.6% elongation, and the negative electrode specimen ruptures under a load of 22 MPa at 1.6% elongation. The failure strain of the electrode specimen is significantly lower than that of the pure current collector metals. One reasonable explanation of this phenomenon is stress concentration on the metal foil caused by small fractures of the active material coating. During elongation test of a tensile specimen, the active material coating will fracture at a strain of one to two percent, because the coating is a mixture weakly held together by polymer binder. As soon as the fracture occurs, almost all further elongation of the specimen is concentrated on the exposed metal substrate which is less stiff than the unexposed portion, leading to rapid rupture of the specimen at the exposed location. Thus the overall elongation of the tensile specimen is much lower than pure metal. This reasoning is used to develop the mechanical models of the battery electrodes implemented in this work, but other explanations such as embrittlement due to the calendaring process are also possible.

Zoom In Zoom Out Reset image size

In order to confirm this hypothesis, finite element models of the tensile test were created, and the output of these simulations are also plotted for Figures 3b and 3c. From comparing the simulation and experimental curves for positive and negative electrode specimens, it is noted that although the peak force and elongation when force peaks in simulation match that of the physical test, the simulation force curves do not drop to zero immediately after the peak, as demonstrated by the test curves. In the physical test, after the active material coating ruptures, the stress is concentrated in the extremely small area of the exposed metal substrate and the entire specimen ruptures immediately. In the simulation, however, due to the finite size of a substrate element, stress is not concentrated on an infinitesimally small area after the coating ruptures. Therefore, the force drops from the peak as the coating ruptures to a lower value, and then to zero as the substrate ruptures after additional elongation. Figure 4 shows a CAE tensile specimen with coating element deletion. The red elements are exposed metal substrate which further elongates during the CAE run until final rupture.

Zoom In Zoom Out Reset image size

Finite element models were also constructed to simulate the compression tests. Initially, the stress-strain curves from the tests are directly used as the material load curves in the finite element model. Then, a scale factor is applied to the finite element material load curve so that the force-displacement curve in simulation matches that of the test. Figure 5 shows calibrated simulation curves versus test data.

Zoom In Zoom Out Reset image size

A linear least squares²³ approach is used to identify the circuit model parameters from the cell-level current and voltage response. As mentioned previously, the Arrhenius dependence of Eq. 11 is assumed for the equivalent circuit parameters to interpolate between, and extrapolate beyond, temperatures other than the measured values. The average resistance (or capacitance) across all state of charge values at a given temperature is first computed. Then, an Arrhenius plot is constructed and the activation energy is identified from the slope of the best fit line to the data at the tested temperatures. The reported temperature dependence of electrical conductivity for aluminum²⁵ and copper²⁶ are applied to the current collectors. The average of C/10 charge and discharge voltages is assumed to represent the open-circuit voltage as a function of state of charge. The identified electrical parameters for 25°C are plotted in Figure 6 and the Arrhenius characteristics shown in Table I.

Zoom In Zoom Out Reset image size

The simulation described in the Case study section is carried out using parallel processing on 48 CPUs. The total wall time is 7 hours, 36 minutes for the full cell model with 150687 elements. Analyzed further, the initial simulation of 0.5 ms physical time takes 1 hour, 48 minutes, and the longer time scale simulation encompassing a physical time of 250 s takes 5 hours, 48 minutes. Clearly these simulation times are dependent upon the time step chosen by the user, and in the case of the fully implicit solution over the long time scale, larger time steps are feasible to reduce the simulation time.

The shorted area time series data are plotted in Figure 7. This simulation output represents the number of elements which have exceeded the compressive strain failure criterion and therefore have been replaced by resistance constraints as a function of time. A constant resistivity as a function of applied pressure is assumed, so the shorted area is inversely proportional to the total short circuit resistance. Both indenters exhibit a period during the initial mechanical deformation where compression occurs without initiating a short circuit. The first failures occur around 0.3 ms for the spherical indenter, and 0.37 ms for the cylindrical indenter. The broader contact area of the cylindrical indenter leads to a more rapid increase in the short circuit contact area, as compared to the spherical indenter. By assumption, the mechanical deformation stops at 0.5 ms, and the electrochemical and thermal responses are driven by the steady state value for the remainder of the simulation. The steady state response is a shorted area of 0.026 m² for the cylindrical indenter and 0.0015 m² for the spherical indenter, or the equivalent of 1040 and 600 local circuit models replaced by shorts respectively. Combining this analysis with the fringe plots of Figure 10 reveals that all layers through the cell thickness short, with greater occurrence (width of the shorted region) in the layers closest to the indenter.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The force versus displacement curves for each indenter, along with nodal displacement fringe plots at the final state, are shown in Figure 8. The characteristic strain hardening observed at the component level is also present when compressing the complete cell. As expected, the strain is concentrated under the contact zone and decays toward zero at a sufficient distance from the failure area. Furthermore, the spherical indenter case, with its more localized contact area, demonstrates a much lower force required to achieve the 0.5 mm compression.

Zoom In Zoom Out Reset image size

The voltage response, as simulated at the cell terminals, is plotted in Figure 9. A nearly instantaneous drop in voltage occurs as the mechanical deformation creates the internal short circuit in less than 0.5 ms. Following the initial drop, dictated by the short circuit magnitude and the cell internal resistance, the equivalent circuit parameters and cell geometry define the voltage response. In the case of the spherical indenter, the short circuit is of a high resistance and the overall discharge rate is slightly less than 1C, which leads to a voltage response that is slightly but steadily decreasing. In the case of the cylindrical indenter, the discharge rate is approximately 10C. This aggressive discharge rate creates immediate and sustained temperature gradients throughout the cell, which causes the electrical resistance parameters to decrease dramatically according to their Arrhenius dependence. The rapid temperature rise overcomes the loss of stored energy that steadily decreases the open-circuit voltage, and the terminal voltage actually increases as a result.

Zoom In Zoom Out Reset image size

State of charge fringe plots are shown for selected time steps in Figure 10. The response is characterized by ohmic limitations of the current collector. For example, the regions nearest the failure area discharge through the short circuit more rapidly than those positioned further away. The cylindrical indenter case experiences a gradient of around 15% state of charge between the regions nearest the indenter and those toward the cell perimeter. The spherical results are not pictured because the overall change in state of charge is nearly negligible due to the slow discharge rate.

Current density within the outermost positive current collectors, viewed from the cell top and bottom surfaces, is shown in Figure 11. In the case of the cylindrical indenter, the failure area bisects the cell into upper and lower regions that are coupled by the resistive constraint created in the failure area. For the sphere, a centralized short circuit location leads to a radial distribution of the current density. In both cases, an interesting pattern is displayed where the top cell layers have a larger short circuit area than the layers further from the indenters. Therefore, the current density near the tabs shows that the layers in the cell interior, through to the bottom layer, actually supply discharge current to the more aggressive short circuit created in the top unit cell.

Zoom In Zoom Out Reset image size

Temperature fringe plots and the nodal temperature time series data for a node at the geometric center of the top surface of the cell are shown in Figure 12. The thermal response is dominated by the ohmic heating emanating outward from the short circuit location, which averages 3300 W for the cylindrical indenter and 300 W for the spherical indenter over the simulation time. By comparing the maximum element temperature predicted by the model with known onset temperatures of exothermic side reactions for a cell design, an engineering safety assessment can be made. Furthermore, since local state of charge variations are incorporated in the model, the dependence of side reactions on state of charge can also be accounted for when assessing safety. In both cases, the temperatures initially rise rapidly before the rate slows due to the free convection imposed on the top cell surface. Figure 13 summarizes the progression of the mechanical and thermal responses for the spherical case to emphasize the separate time scales involved.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

There are several areas of future work that are critical to validating the simulation capability and to move from simulation development toward virtual verification. First, the understanding of the configuration, area, and contact resistance of the internal short circuits that are formed is crucial. This improvement hinges upon both high fidelity mechanical models and experimental methods for creating and measuring internal short circuits that are representative of those that occur during blunt impact loading. The compressive strain criterion is a simple, easily implemented failure criterion that was chosen to demonstrate a method for one-directional coupling for this study, but experimental studies suggest the mechanics responsible for forming an internal short circuit are more complex. The simulation format used in this paper is easily extendable to short circuits defined by other means, such as element deletion due to tensile or compressive failure, and further combined modeling and experiments will be required to select a robust set of failure criteria. Though an assumed constant value for short circuit contact resistance was used in this work, the simulation can accept a tabular function of contact resistance versus contact pressure which could also be used to parameterize the resistance values. Second, a more detailed investigation into the element formulations and sizes that define the mechanical model is warranted, to maintain accuracy while reducing computation time. As previously mentioned, thick shell elements with strain reduced integration are deemed a reasonable first choice, as they are recommended for applications featuring a high aspect ratio and coarse through-thickness discretization. To confirm this element formulation choice, full cell testing will be performed with measurement of the indenter contact force and displacement during a variety of mechanical loading conditions. The elements used in the presented simulation exhibit aspect ratios as high as 250. These poor aspect ratios lead to hourglass modes at higher degrees of crush, so further reduction of the mesh size in the indenter contact zone may also prove useful for increasing the robustness of the simulation. However, increasing the number of elements will increase the computational cost per time step. Therefore, methods for creating transition meshes which employ elements with moderate aspect ratios in areas that experience intrusion and high aspect ratios further away from the indenter are needed to balance robustness with computational time. Given the relatively modest computation time of the presented model, a finer mesh could be used and still be tractable. These aspects will be explored in future work as the simulation tool is validated with experimental data for impact scenarios.