Calibration of Crushable Foam Models for the Jellyroll of Cylindrical Lithium-Ion Batteries

Abstract

Crushable foam plasticity models are employed to simulate material response under essentially monotonic loading. For the plastic part of the behavior, the default crushable foam model in Abaqus/Explicit is the volumetric hardening model, where the yield surface evolves by the volumetric compacting plastic strain, and the other available model is the isotropic hardening model, where the yield curve is centrally located at the origin in the pressure—the Mises stress plane. In this study, the characteristic of two models was examined by applying them to a simple 18650 lithium-ion cylindrical cell. The computation cell model consists of the shell casing and the homogenized jelly roll which represents the electrode assembly. Both crushable foam models were calibrated to represent the homogenized mechanical properties of the jellyroll, and the load–displacement relations were compared with the experimental results. Then, we examined the deformation characteristic of jellyroll for each crushable foam model.

1. Introduction

Nowadays the growing demand for an electric vehicle battery in the automotive industry requires a high energy density. As the energy density of lithium-ion batteries (LIBs) is continuously increasing, abusive environments might cause the battery to suffer from an exothermic phenomenon from a sudden internal short circuit [1,2,3,4], and it might lead to a thermal runaway such as fire or explosion [5,6,7,8]. Even though the possibility of its failure in actual use has a rare chance to occur, one of the critical challenges in LIBs is the safety issue. To address a range of possible abuses of LIBs, a number of standards and testing protocols have been developed by Underwriters Laboratories (UL). UL regulation [9] is currently used to assess specific safety risks associated with designed short-circuit abuse environments categorized as external short, electric, mechanical, and thermal abuses. Since the integrity of cell components or battery package systems might be deteriorated by events such as electric vehicle crashes [10,11], the mechanical abuse protocols of them are designed to damage LIBs intentionally by crushing [12], impacting [13,14], or penetrating a nail [15,16,17]. After the test proceeds, the safety level of LIBs is evaluated. To develop a cell or a pack satisfying the safety level, we need to understand the internal short-circuit mechanism. However, a cell or a pack is smoked and/or burned out completely due to the electrochemical reactions. Therefore, it is hard to predict the progressive sequence of local damage. To solve this problem, computational modeling analyses have been used, which require both a detailed structure modeling and the mechanical properties of each component [18]. However, a commonly used jellyroll in LIBs is manufactured through a winding process that involves a roll of layers consisting of anode/separator/cathode/separator. Replicating the complex structure and material layers in a simulation requires not only expensive preprocessing but also a high computational cost. As the scale of the analysis model increases, this inefficiency increases since the actual number of battery cells mounted on a vehicle is not just one, but hundreds.

As another simulation approach for analysis, a roll of layers in LIBs is considered as a homogenized material. A homogenized model needs fewer material parameters than a detailed model does, and modeling is simple. Owing to these advantages, numerous studies are conducted based on the homogenized model [19,20,21,22,23]. This approach introduces the concept of equivalent mechanical properties. Equivalent properties are derived by assuming a material yield model and then using the experimental load–displacement curve to correct the property parameters so that the experimental and analysis results match. Several authors have conducted investigations to find the material response of the cell under mechanical abuse conditions. For safety assessments, Greve and Fehrenbach [24] performed a quasi-static mechanical abuse test on cylindrical LIBs. They found that LIBs are compression-dependent characteristic. A macro-mechanical finite element crash simulation model was proposed considering both a pressure-dependent yield criterion and a pressure-dependent non-associated flow rule. E. Sahraei et al. [25] performed extensive mechanical testing on the most commercially available 18650 LIBs with a nominal capacity of 2600 mAh, including a lateral indentation of the cell by a rigid rod, crush test by a hemispherical punch, compression of the cell between two plates, and three-point bending test, and they developed an FE model to predict short-circuit detection. The stress–strain relation for the jellyroll was calculated by E. Sahraei et al. [19] for the case where the cell was crushed between two flat plates, estimating the load transfer mechanism inside the cell, and the constitutive behavior of the jellyroll in compression showed the parabolic function of the stress. Then, they developed a finite element model of the cell by using the isotropic crushable foam material in LS Dyna [25], in which the stress is a nonlinear function of the strain.

The constitutive modeling of the crushable foam was originally developed for polymeric and metallic foams showing a large compressive strength and a small tensile strength, since their Poisson ratio is small due to their pore microstructures [26,27,28,29,30,31]. Therefore, it is appropriate to simulate material response under monotonic loading. The standard crushable foam model in Abaqus/Explicit is the volumetric hardening model, where the yield surface evolves by the volumetric compacting plastic strain, and the other available model is the isotropic hardening model, where the yield curve is centrally located at the origin in the pressure—the Mises stress plane [32]. Both crushable foam models require proper material parameters to simulate the elastic/plastic behavior of an LIB cell under abusive environments. In this study, the calibration method is suggested to find the material parameters for both models. The proposed method applies to a simple 18650 cylindrical cell, which consists of the shell casing and the homogenized jellyroll which represents the electrode assembly. Both crushable foam models were calibrated to represent the homogenized mechanical properties of the jellyroll. For verification, two simulations were conducted, such as cell indentation between two rigid plates and cell indentation by a rigid cylinder rod. By comparing their results with the load–displacement curve of the experimental results, the calibration procedure was verified. Although a simple 18650 cell was used in this research, the proposed method is general, so we expect that it can apply to any type of cylindrical cells or prismatic cells. Further, we can expect that this method can be extended to model a detailed cell consisting of a roll of layers consisting of anode/cathode/separators to track the plastic behavior of each layer.

2. Numerical Modeling of LIB in Quasi-Static Conditions

2.1. Finite Element Model in Abaqus/Explicit

Two indentation problems are considered to illustrate the calibration procedures of crushable foam materials. One case is a deformable 18650 cylindrical LIB indented by a rigid square plate, and the other is the same kind of battery indented by a rigid, cylindrical punch. The 18650 cylindrical LIB consists of an outer shell casing and a jellyroll. The cell will deform symmetrically, so it is enough to analyze only half of it to save computational cost as shown in Figure 1.

The outer shell casing was modeled with the homogeneous shell elements without considering an endcap for simplification. The outer shell casing was made from deep drawn steel, and it had a height of 65 mm, a diameter 18 mm, and a wall thickness of 1.8 mm. It was modeled with S4R 4-node curved general-purpose reduced integration shell elements with hourglass control in Abaqus/Explicit. Contact conditions were defined between the shell casing and the inside of the jellyroll. The jellyroll assumed to be a crushable foam material was contained within the shell casing. The height of the jellyroll was 58 mm, and its diameter was 18 mm. It was modeled with C3D8R 8-node liner brick and reduced integration solid elements with hourglass control in Abaqus/Explicit.

Both the rigid floor and the indenter were modeled with R3D4 4-node 3-D bilinear rigid quadrilateral elements in Abaqus/Explicit, and the rigid cylindrical indenter was modeled as a 3D analytic rigid surface. Several different sizes of meshes and element types were tried in an iterative manner to match the localized deformation amount, and the section properties and mesh information are listed in

Table 1. Section properties and mesh information.

Parts	Element Type	Mesh Size (mm)	Element Number
Casing (without cap)	Shell	1	3360
Jellyroll	Solid	1	17,864
Indenter/rigid floor	Discrete rigid	5	400
Cylindrical indenter	Analytic rigid	-	-

.

The contact between the outer surface of the shell casing and the rigid indenter was modeled with a contact pair. Both the square and the cylindrical rigid indenters were modeled as analytical rigid surfaces using a surface definition in conjunction with a rigid body constraint. The contact between the inner surface of the shell casing and the jellyroll was modeled with a contact pair. Coulomb friction was modeled between parts with a friction coefficient of 0.3.

The indenter was fully constrained except in the vertical direction, and a reference point was set up to measure the reaction force when its motion was displaced quasi-statically when the cell was modeled without rate dependence. An amplitude curve with smoothing was used to specify the displacement of the indenter and to promote a quasi-static solution. The indenter’s loading rate was 0.35 m/s, in which the kinetic energy for the entire structure was a negligible fraction of the internal energy. In these simulations, an HP-Z840 workstation with 6 CPUs was used.

2.2. Materials

The material parameters for the shell casing were $E=207 \mathrm{GPa}$ (Young’s modulus), $\nu =0.3$ (Poisson’s ratio), and $\rho =7.85 {\mathrm{g}/\mathrm{cm}}^3$ (density), and the plastic property of the shell casing defined by a function of the equivalent plastic strain ${\overline{\epsilon }}^p$ is as follows [24]:

$$R_{\mathrm{housing}}\left({\overline{\epsilon }}^p\right)=700{(0.00801+{\overline{\epsilon }}^p)}^{0.1385} \mathrm{MPa}.$$

(1)

The crushable foam model was used for the jellyroll of the 18650 LIB. In this model, the elastic part of the response can only be specified as linear isotropic elasticity. For the plastic part of the behavior, the yield curve is an ellipse in the hydrostatic ($p$) and Mises ($q$) stress plane, and two hardening models in Abaqus/Explicit are available. For the volumetric hardening model, in which a hydrostatic tension loading is fixed, the yield surface evolves by the volumetric compacting plastic strain so that the ellipse grows only in the positive-pressure direction, and for the isotropic hardening model, the yield surface evolves symmetrically about the center of the ellipse, as shown in Figure 2. Deshpande and Fleck [14] developed the isotropic model for metallic foams.

The initial yield surface of the volumetric hardening model is defined by the ratio ($k={\sigma }_c^0/p_c^0$) of the yield stress in uniaxial compression (${\sigma }_c^0$) to the yield stress in initial hydrostatic compression ($p_c^0$), and by the ratio ($k_t=p_t/p_c^0$) of the yield stress in hydrostatic tension ($p_t$) to the yield stress in hydrostatic compression ($p_c^0$). The yield surface of the volumetric hardening is defined as follows:

$${\left[\alpha {\left(\frac{p_t-p_c}{2}+p\right)}^2+q^2\right]}^{1/2}=\alpha \frac{p_c+p_t}{2},$$

(2)

and

$$\alpha =\frac{3k}{\sqrt{(3k_t+k)(3-k)}} \mathrm{with} k=\frac{{\sigma }_c^0}{p_c^0} \mathrm{and} k_t=\frac{p_t}{p_c^0},$$

(3)

where $q$ is the Mises equivalent stress, $p$ is the hydrostatic pressure, and $\alpha$ is the shape factor of the yield ellipse. The flow potential is chosen as follows:

$$h=\sqrt{\frac{9}{2}{p}^2+q^2}.$$

(4)

The hardening curve can be obtained from uniaxial compression test data as a function of the value ($-{\epsilon }_{\mathrm{vol}}^{pl}$) of volumetric compacting plastic strain (PEEQ):

$$p_c\left({\epsilon }_{\mathrm{vol}}^{pl}\right)=\frac{{\sigma }_c\left({\epsilon }_{\mathrm{axial}}^{pl}\right)\left[{\sigma }_c\left({\epsilon }_{\mathrm{axial}}^{pl}\right)\left(\frac{1}{{\alpha }^2}+\frac{1}{9}\right)+\frac{p_t}{3}\right]}{p_t+\frac{{\sigma }_c\left({\epsilon }_{\mathrm{axial}}^{pl}\right)}{3}},$$

(5)

with the fact that ${\epsilon }_{\mathrm{axial}}^{pl}={\epsilon }_{\mathrm{vol}}^{pl}$ in uniaxial compression for the volumetric hardening model.

For the isotropic hardening model, the initial yield surface is only defined by the ratio ($k={\sigma }_c^0/p_c^0$), but the plastic Poisson’s ratio ($v^p$), which is the ratio of the transverse to the longitudinal plastic strains under uniaxial compression, needs to be additionally provided. The yield surface for the isotropic hardening model is defined as follows:

$$\sqrt{q^2+{\alpha }^2{p}^2}={\sigma }_c\sqrt{1+{\left(\frac{\alpha }{3}\right)}^2},$$

(6)

where

$$\alpha =\frac{3k}{\sqrt{9-k^2}} \mathrm{with} k=\frac{{\sigma }_c^0}{p_c^0}.$$

(7)

The flow potential for the isotropic hardening model is chosen as follows:

$$G=\sqrt{q^2+{\beta }^2{p}^2} \mathrm{with} \beta =\frac{3}{\sqrt{2}}\sqrt{\frac{1-2{\nu }_p}{1+v_p}}.$$

(8)

The hardening law defines the yield stress in uniaxial compression as a function of the absolute values of the axial plastic strain [32].

For both hardening models in this study, uniaxial compression test data are input as the following function [19]:

$${\sigma }_c=0.8+848{\left({\overline{\epsilon }}_{\mathrm{axial}}^p\right)}^{2.7} \mathrm{MPa} .$$

(9)

This function is obtained from the following plastic flow rule [24]:

$$R({\overline{\epsilon }}^p)=\left[{\sigma }_{\mathrm{plateau}}-({\sigma }_{\mathrm{plateau}}-{\sigma }_{\mathrm{yield}})\times \mathrm{Exp}\left[-\frac{{\overline{\epsilon }}^p}{{\overline{\epsilon }}_{\mathrm{ref}}}\right]\right] \times \left[1+s{({\overline{\epsilon }}^p)}^m\right] \mathrm{GPa},$$

(10)

where the calibrated parameters are ${\sigma }_{\mathrm{yield}}=0.0003 \mathrm{GPa},$ ${\sigma }_{\mathrm{plateau}}=0.0006 \mathrm{GPa},$ ${\overline{\epsilon }}_{\mathrm{ref}}=0.009,$ $s=1000$, and $m=2.7$.

2.3. Calibration of Crushable Foam Material Parameters

For the volumetric hardening model, the rate-independent calibration was performed as follows. (1) The elastic properties: the Young’s modulus $(E)$ obtained from the calibration along with the plastic strain hardening curve was used, and the elastic Poisson’s ratio $(\nu )$ was assumed to be almost zero, since the major part of the porous microstructure of a coating allows for substantial compressibility. (2) The crushable foam properties: the uniaxial compression data in Equation (9) were used for the crushable foam hardening, and the compression yield strength ratio ($k={\sigma }_c^0/p_c^0$) was increased to a value within $0<k<3$ until the simulated result of the force–displacement curve for the crush test between two plates agreed with experimental one. In addition, the tension yield strength ratio ($k_t\ge 0$) is usually assumed to be the magnitude of the strength of the foam in hydrostatic tension, since jellyroll is hardly tested in tension. In this study, the tensile strength ($p_t$) was assumed to be within five to ten percent of the initial yield stress in hydrostatic compression ($p_c^0$), which is always positive. Then, the deformed geometries for the rod indentation were checked to find the proper minimum value of the ratio ($k_t$). After using compression and tension yield strength ratios, the force–displacement for the rod indentation simulation was compared with the experimental one. If the results did not fit each other, the analyses ran with the increased compression yield strength ratio until both numerical results corresponded closely with the experimental results.

For the isotropic hardening model, the elastic properties, the compression test data for the hardening, and the compression yield strength ratio are required, so the same values as the volumetric hardening model’s ones were used. The plastic Poisson’s ratio $(v_p)$ for the isotropic hardening model is required instead of the hydrostatic compression. For the parameter calibration, the compression yield strength ratio which is found by the volumetric hardening model was used, and then the analyses were run with different plastic Poisson’s ratios until both numerical results corresponded with the experimental results.

3. Results and Discussion

The true stress and logarithmic strain values in Figure 3a are used to represent the uniaxial compression curve for the three-dimensional isotropic and volumetric cylindrical cell models. For the indentation simulation with the two identical plates, the jellyroll’s material parameters for both models were calibrated by using the procedures described in the previous Section 2.3, and they are listed in

Table 2. The calibrated parameters for the jellyroll.

Parameter	$\rho {(\mathrm{g}/\mathrm{cm}}^3)$	$E \left(\mathrm{GPa}\right)$	$\nu$	$k$	$k_t$	${\nu }_p$
Value	2.27	1.5	0	2.5	0.1	0.25

. For the volumetric model, the tensile strength was approximated by setting the hydrostatic tension equal to 10% of the initial yield stress in hydrostatic compression, since it should not have a strong effect on the simulation, and then the yield strength ratio within the possible range was found iteratively. For the isotropic model, the plastic Poisson’s ratio was calibrated by fixing the same yield strength ratio as the one used in the volumetric model. For verification purposes, the load–deflection responses were compared with the experimental results given by Sahraei et al. [19]. The numerical results plotted in Figure 3b show that the simulated load–deflection responses for both models were in good agreement with the experimental results.

The meshes in the three-dimensional model were created using the reduced integration with hourglass control, and they were refined until the load–displacement curve approached a certain limit. At the end of the loading step, the rigid indenter was displaced by about 6.5 mm. The deformed configurations of the cylindrical cell with the volumetric and isotropic hardening models are shown in Figure 4a,b, respectively. For both cases, the maximum stresses occurred at the outer shell case, and the permanent plastic deformation appeared along the region in contact with the indenter. The crushing strength of the casing did not have a strong influence on the load–displacement response even though the thickness increased from 0.18 mm up to 0.3 mm, so that the jellyroll had a major effect on the compressive response compared to the outer shell casing. The stress distributions of the jellyroll at the end of the loading step are shown in Figure 4c,d, respectively. For the volumetric model, the uniform maximum stress occurred along the center of the jellyroll in the longitudinal direction, while the maximum stress occurred at both ends of the core zone of the jellyroll for the isotropic model.

Figure 5a,b show the contour plots of the volumetric compacting plastic strain (PEEQ) for the volumetric hardening model and the equivalent plastic strain (PEEQ) for the isotropic hardening model, respectively. The figures show that the plastic strain magnitude uniformly approaches 35% along the core of the jellyroll for the volumetric model, while the maximum equivalent plastic strain approaches 36% at both ends of the core of the jellyroll for the isotropic model. The equivalent plastic strains for both cases reach the permanent deformation, so that the short circuit initiates from the location with the maximum plastic strain. The simulation results are in good agreement with the experimentally observed short-circuit locations [24].

For further verification, both models were applied to indentation with the cylindrical rod, as shown in Figure 1b, with the above calibrated material parameters, and the simulated load–displacement results for both cases were compared with the experimental one by Sahraei et al. [19]. The numerical results plotted in Figure 6a show good agreement with the experimental results.

In this simulation, the rigid indenter was displaced by about 6.2 mm at the end of the loading step. The deformed geometries for indentation with the cylindrical rod are shown in Figure 7a,b, respectively. The maximum stresses for both models occur at the edge of the outer shell contacting the rod indenter, and the magnitude of both models shows similar values, while PEEQ (20%) for the volumetric model is 2% higher than PEEQ (18.6%) for the isotropic model. The stress distributions of the jellyroll at the end of the loading step for both models are shown in Figure 7c,d, respectively. The magnitude of maximum stress is similar for both models, and the local stress concentration occurs right below the rod indenter.

Figure 8 shows the contour plots of the volumetric compacting plastic strain (PEEQ) for the volumetric hardening model and the equivalent plastic strain (PEEQ) for the isotropic hardening model, respectively. The figures show that the compacting plastic strain in the vicinity of the punch to the bottom of the jellyroll uniformly approaches 36%, while the equivalent plastic strain approaches 37% in the vicinity of the punch and the bottom. Therefore, for the rod’s compressive loading simulation, the cell fracture will start right below the rod indenter or at the bottom of the cell. The predicted result shows a good agreement with the experimental infrared measurement [24].

In addition, a case study was carried out to investigate the characteristics of the crushable foam hardening model. The yield strength in hydrostatic tension was set as 5% of the initial yield stress in hydrostatic compression, and then the strength ratio k was calibrated in an iterative manner in order to find the same load–displacement curve for the indentation experiment between two rigid plates. By using the calibrated values of k and k_t for the volumetric model, the numerical result showed a good agreement with the experimental result as shown in Figure 6b. Although two different material parameters were used for the volumetric model, the permanent plastic strain and the deformed configuration in Figure 9a are the same as the ones in Figure 9b for the indentation simulation. However, for the cylindrical rod indentation, the deformed configuration showed an unreasonable, large deformation in the longitudinal direction as shown in Figure 9d. Five percent of the initial yield stress in hydrostatic compression led to less resistance than ten percent did for the radial direction. Therefore, we can conclude that the calibration procedure for the parameters of the crushable foam model is required to check whether the force–displacement responses for all loading cases are consistent with the experimental results. Similarly, the isotropic model also showed the same characteristics for the different plastic strain (v_p) values. Therefore, the calibration procedure for all loading cases is also needed.

4. Conclusions

The crushable foam plasticity model was developed for the analysis of energy absorption structures or to simulate material responses under a compressive loading. In order to represent the mechanical responses on a cylindrical cell under abusive environments, the crushable foam model has recently been used. To use this model efficiently, we need to know characteristic of the plastic part of the behavior. In this study, calibration procedures for the plastic behavior of the cell were investigated with two different crushable foam models available in Abaqus/Explicit. To calibrate the material parameters for both crushable models, the experimental load–displacement curve of the mechanical crush test was used. By performing the simulation of a cylindrical rod indentation and comparing it with the experimental result, the calibration procedure was validated. The deformed configurations of the volumetric and isotropic hardening foam models showed similar results obtained from the experimental infrared images. Therefore, the calibration procedure for both crushable foam models could be used to predict the mechanical responses of the cylindrical cell under compressive loading conditions.

The proposed calibration procedure of the crushable foam models encourages us to develop a detailed modeling of a cell consisting of a set of anode, cathode, and separators. The simulation of the detailed modeling could predict the weakness of each layer, which would be helpful to design a cell.