Homogenized modeling methodology for 18650 lithium-ion battery module under large deformation

Effective lithium-ion battery module modeling has become a bottleneck for full-size electric vehicle crash safety numerical simulation. Modeling every single cell in detail would be costly. However, computational accuracy could be lost if the module is modeled by using a simple bulk material or rigid body. To solve this critical engineering problem, a general method to establish a computational homogenized model for the cylindrical battery module is proposed. A single battery cell model is developed and validated through radial compression and bending experiments. To analyze the homogenized mechanical properties of the module, a representative unit cell (RUC) is extracted with the periodic boundary condition applied on it. An elastic–plastic constitutive model is established to describe the computational homogenized model for the module. Two typical packing modes, i.e., cubic dense packing and hexagonal packing for the homogenized equivalent battery module (EBM) model, are targeted for validation compression tests, as well as the models with detailed single cell description. Further, the homogenized EBM model is confirmed to agree reasonably well with the detailed battery module (DBM) model for different packing modes with a length scale of up to 15 × 15 cells and 12% deformation where the short circuit takes place. The suggested homogenized model for battery module makes way for battery module and pack safety evaluation for full-size electric vehicle crashworthiness analysis.


Introduction
With the strong support from the government [1] and major technology breakthrough for lithium-ion batteries (LIBs) [2,3], electric vehicles (EVs) have been witnessed to boom over the past recent years [4][5][6]. The major reason for LIBs to become the primary choice for EVs is due to the combination advantage of high energy/power density, lightweight, and safety [6,7]. Many research works were conducted on the optimization approach for electrified vehicles to increase cost competitiveness and reduce carbon emissions [6,[8][9][10]. Additionally, the battery management system plays an important role in maintaining battery lifetime without unduly sacrificing its performance. Some key technologies of battery management system that monitor the unmeasurable internal states of the battery have been extensively studied [11][12][13][14]. PLOS  Battery models such as electrochemical [13], thermal [15], and high-order physics-based model [16] are also used for evaluation of existing charging strategies, estimator and controller development, simulation, and optimization. However, none of these works have considered the mechanical performance of the battery. Because of the inevitable crash or impact for vehicles during traffic accidents, there is a high possibility of internal short-circuit [15], thermal runaway [17], and even explosion [18,19] for LIBs subjected to external mechanical loading. Thus, the crashworthiness issue of EVs by consideration of LIB safety remains a paramount concern in electric vehicle safety.
In the past years, many pioneering efforts were made to elucidate the mechanical safety behavior of LIBs on multiple length scales, ranging from component material scale to battery pack scale. On the component level, the mechanical properties of the case [20], separator [21], anode, and cathode foils [22] of the cell have been investigated. These results will help to understand the failure mechanism of internal short circuit of lithium-ion battery, but can not characterize the global mechanical behavior of the battery. Many experiments were conducted for the investigation of mechanical behaviors of the battery cells under various loading conditions [20,23,24]. The finite element (FE) models were used to understand the mechanical properties and predict extreme cases. Since coating of active material coating and separator are highly porous and soaked in the electrolyte, the detailed modeling including each component and interaction among them is very complicated. Therefore, it is reasonable and acceptable to take jellyroll as a homogeneous material [23,25,26]. Furthermore, the dynamic behavior [27] and SOC effect [28] of the battery have been studied and the results suggested that higher SOC leads to higher structure stiffness. In this paper, a computational model of a single cylindrical battery is established and validated based on homogeneous modeling technique.
For battery module and pack, the mechanical safety performance is closely related to sizes and packing modes of the module and modeling every single cell in detail would be costly. Therefore, equivalent homogenized model for the module is fully necessary for many applications such as vehicle level crashworthiness analyses and optimum design. Sahraei et al. [29] used a homogenized crushable foam core to simulate the interior containing battery cells of an idealized battery pack to model the drop test. The battery pack was taken as a linear elastic material in [30] for crash analysis of a conceptual electric vehicle. However, modeling the module by using a simple bulk material would result in sacrificing computational accuracy. Additionally, Lai et al. [31] adopted macro homogenized material models calibrated by the test data to simulate the punch test of a small-scale module specimen. Nevertheless, it is very difficult to directly measure the integrated mechanical properties of the battery module.
This paper responds to the challenge by extracting a representative unit cell (RUC) with the periodic boundary condition applied on it to analyze the homogenized mechanical properties of the module. An elastic-plastic constitutive model is established to describe the equivalent battery module (EBM) model. Further, a small-scale battery module is tested to compare the mechanical behavior with those obtained from the EBM model and detailed battery module (DBM) model. Upon the satisfactory comparison results, EBM model is further generalized for larger battery modules and different packing modes and the feasibility of the established model is discussed.

Computational model of a single cell LIB
Mechanical behavior of a single cell LIB surely dominates that of a battery module. The 18650 LixC6/LiCoO2 batteries are used in this study, provided by SONY. The main cell geometry and dimensions are provided in Fig 1A. The cylindrical LIB cell is composed of winding cathode, anode, and separator films and encapsulated into a steel can. So it is extremely complicated to simulate each component and interaction among them in detail. Thus, the entire LIB is taken as an orthotropic homogeneous material. For simplicity, the x1-x 2 plane that is perpendicular to the axis of the battery is considered as isotropic and a transversely isotropic elastic plastic model has been adopted to express the mechanical properties.
The stress-strain laws for transversely isotropic elasticity are shown as follows: where p and t stand for "in-plane" and "transverse," respectively. Thus, E p and E t are Young's modulus of radial and axial directions, E 1 = E 2 = Ep and E 3 = E t ; G p and G t are the shear modulus of the corresponding planes, G 12 = G p and G 13 = G 23 = G t ; ν p , ν tp and ν pt are the Poisson's ratio, ν 12 = ν p and ν 13 = ν 23 = ν pt ; and subscripts "1", "2" and "3" stand for the coordinate system. The anisotropic yield behavior is modeled through yield stress ratio, R ij , which is applied in Hill's potential function. The yield stress ratio is defined as follows: where each " s ij is the measured yield stress value; σ 0 is the user-defined reference yield stress set The plastic-hardening model for the compression is expressed as the form of the Eq. (3) proposed by Greve and Fehrenbach [23]: where A and B are the parameters to be determined, σ is the stress, ε is the plastic strain, and n is the hardening exponent.
The above developed mechanical model is calibrated through the radial compression and axial compression tests, as shown in Fig 1B and 1C. The compression tests were conducted on INSTRON 8801 platform with a loading rate of 2 mm/min, where the indenter and the bearing were both flats. The force and displacement sensors were set on the indenter. For radial compression, the nominal stress-strain curve can be converted from the load-displacement curve as follows: where F is the load and S c stands for the contact area; l c is the length of the cell and b c is the contact width; R is the radius of the cell and s is the displacement of the indenter. The radial compressive nominal stress-strain curve distinctly has three stages, as follows: linear region, plateau region, and densification region. The axial modulus is extracted from the second stage of the curve because the cap collapsed during the axial compression test. The parameter values for the mechanical model of a single cell are summarized in Table 1. To validate the model, the simulations of radial compression and bending were conducted based on the ABAQUS/Explicit platform on which the proposed constitutive model was used. The bending experiment data were cited from open literature and more detailed information could have been referred from another study [26]. The computational model of the 18650 lithium-ion battery cell was composed of eight-node linear brick, reduced integration solid element with the size of 0.8 mm. As illustrated in Fig 2A and 2B, the model can well predict the mechanical behavior of the cell under radial compression and bending loading conditions, where the coefficient of determination ðy i À À yÞ 2 ¼ 0:97 (where y i and f i are the experimental and simulated data, respectively, and " y is the average of the experimental data).

RUC extraction for the battery module
The packing method of cylindrical cells can be treated as a geometrical model, where packing sizes and modes of the module were uniformly defined with b, l, and θ as shown in Fig 3. b and l are the number of cells on each row and column and θ describes the relative position of rows,  which varies from 2π/3 to π. The packing density of the module can be expressed as follows: Since the battery module can be regarded as the periodic arrangement of RUCs, RUC of the packing module is extracted. As shown in Fig 2, RUC is related to θ, where the nominal density should be calculated as the following: Naturally, ρ RUC equals the packing density of a module when b and l are regarded as infinite. In this paper, the mechanical properties of RUC are assumed to be equivalent to the properties of the module. When θ = π and θ = 2π/3, i.e., cubic dense packing and hexagonal dense packing, the RUC is symmetrical in structure, and the structure of module is stable during compression deformation.

Computational model for RUC
For a battery module, the length of the cell is sufficiently long relative to its axial deformation, and every cross section of the module perpendicular to the axes of the cell is identical.  Therefore, it is reasonable to be simplified into a two-dimensional plane strain problem. In order to verify the simplification from 3D to 2D, constrained compression simulations using 3D model and 2D model, respectively, were performed, as shown in Fig 4. The results of 2D model and 3D model are almost the same with each other. Then, 2D model is chosen for the next works.
2D model was developed to represent the 3D RUC and the RUC FE model was composed of four-node bilinear quadrilateral plane strain, which reduced integration 2D solid elements with the size of 0.8 × 0.8 mm. A penalty based contact was set for the contacting parts of adjacent cells by reasonably assuming a friction coefficient of 0.1. To understand the compressive mechanical properties of RUC, a periodic boundary condition that has been investigated in Ref [32][33][34][35] was applied to the boundaries of RUC and the boundary equation, as follows: where u is the displacement of note, " ε ik are the average strains, index ''j+" means along the positive X j direction and ''j−" means along the negative X j direction and Dx j k are distances between pairs of opposite boundary surfaces. The top and bottom sides were forced and fixed separately to achieve the mechanical properties.

Homogenized method for battery module modeling
A homogenized method is proposed to replace the detailed modeling method for the battery module. The mechanical properties of the homogenized EBM model can be equivalent to the properties of the RUC model. Fig 5 shows the true stress-strain curves of RUCs for two typical packing modes during uniaxial compression deformation. The true stress-strain curve could be converted from the nominal stress-strain curve as follows: where ε nom and σ nom can be calculated from the load force and displacement divided by the contact area and the original height of the RUC. The curves first increase linearly to a plateau stage and then increase sharply after the plateau stage. Apparently, there is a difference between the two curves that the stress for hexagonal dense packing is higher than that of cubic dense packing because of its larger density. An elastic-plastic material model can be used to predict the homogenized EBM behavior. Young's modulus E (θ = π) and E (θ = 2π/3) are calculated from above true stress-strain curves as 55 and 90 MPa, respectively. The plastic stage of the curves are used as input file of the plastic part of the material model.

Experiment setups for EBM model verification
To validate the proposed homogenized EBM model, the constrained compression experiments in quasi-static loading condition that commonly occurred during traffic accidents were conducted on INSTRON 8801 platform with a loading rate of 2 mm/min as shown in Fig 6. The battery modules included nine fully discharged battery cells, which are in contact with one another for typical cubic dense packing and hexagonal dense packing. Connection tabs and cooling systems existing in a realistic battery module were neglected. The width and height of the module were within 65 mm in a small size. Moreover, the battery module was restrained by two steel plates with thickness of 20 mm. The distance between the two plates could be adjusted to the size of the module and the plates are fixed by bolts. For reduction of friction effect, the width of load transfer punch is slightly smaller than the distance between the plates. The battery module was compressed by 10 mm, which was about 20% of its height.

Results
Experiments, homogenized EBM model and corresponding DBM model were conducted on the case for the cubic dense packing mode, i.e. θ = π. The settings of boundary conditions for both computational models were exactly the same as in the experiment.  For the EBM computational model, four-node bilinear quadrilateral plane strain, which reduced the integration of 2D solid elements, was used; this was also verified via convergence study. Moreover, the proposed homogenized mechanical properties in Section 2.4 were adopted. For the DBM model, nine single battery cell models were included that developed in Section 2.1. The supporting platform and the indenter were set as discrete rigid. The indenter was designated a certain displacement with quasi-static compressive loading. A penalty-based contact was set for the contacting parts by reasonably assuming a friction coefficient of 0.1.
The stress evolution of the module during deformation is shown in Fig 7. The stress distribution of EBM model is mean, and the corners near the load are under high stress state. For the DBM model, the stress wave spreads from contact area along the load direction to the rest part of the cell, and the stress values of contacting areas of the neighboring cells are higher. In addition, the maximum stress value of the EBM model is smaller than that of DBM model under the same deformation because of larger contact area. It should be noted that, as the deformation increases, the stress growth of the EBM model slows down. As presented in the subsequent section, in the overall mechanical response, the deviation of the load between the EBM model and the DBM model will increase dramatically when the deformation exceeds a certain range. Therefore, the homogenized EBM model can only precisely predict the mechanical behavior of a battery module within a certain deformation range.
Comparison of the mechanical responses of the battery module for cubic dense packing is shown in Fig 8A. The load-displacement curves are close to one another (where r 2 = 0.92-0.96), thereby indicating that the EBM model can well predict the integral mechanical behavior of the battery module under constrained compression loading conditions on a small size. It is noted that the EBM model could not distinctively simulate the local mechanical behavior such as contacts between neighboring cells.
Additionally, as illustrated in Fig 8, the DBM model can well reflect the mechanical behavior and deformation of the module for cubic dense packing (where r 2 = 0.96) and hexagonal dense packing (where r 2 = 0.94), so it could be used as a method to validate the EBM model if it is difficult to conduct an experiment.

Applicability for the battery module under large size
In general, there are more than nine cells in a realistic battery module. Therefore, the applicability of the homogenized modeling method for the module should be validated under large size. The EBM model and DBM model under different sizes were conducted in constrained compression loading condition for two typical packing modes, i.e., cubic dense packing and hexagonal dense packing. Because the maximum force generated during experiments for the module including nine cells already reached almost 100 KN load limit of the INSTRON 8801, the validation under larger size will be conducted by the DBM model.
For the cubic dense packing mode, the packing density of a module is equivalent to that of corresponding RUC under different sizes. Fig 9 shows the comparisons of mechanical responses between the EBM model and DBM model under different sizes. For observation convenience, the mechanical responses are expressed as the form of load-strain curves and the subscript "b × l" denotes that a number of b and l cells are placed in each row and column, respectively. It is found that there is a good agreement between the EBM model and DBM model under each size (where r 2 = 0.94-0.99), which indicates that the mechanical behavior can be predicted well by this model with length scale up to 15 × 15 in the range of 15% deformation where the structure of cubic dense packing is stable. When the size reaches a greatly extent, the slight dislocation of individual cells in the module will occur due to deformation of cells, which leads to a load decline.
For the hexagonal dense packing mode, the packing density of a module increases as the increase of its size b and l from Eq. (8). Fig 10 shows the comparisons of mechanical responses between the EBM model and DBM model as well as the packing density of the module under different sizes. It is observed that the difference between the EBM model and DBM model becomes small gradually as the size increases (where r 2 increases from 0.30 to 0.99). Since the effect of size could be eliminated according to the aforementioned result of cubic dense packing, it is reasonable that the EBM model can more precisely describe the module when the density of the module approaches that of its RUC.
When the battery module is taken as a homogenized composite where cells and air are seen as fibers and matrix, respectively, its elastic modulus depends on the volume fraction and elastic property of components according to the composite rule of mixture. Therefore, the elastic

Applicability for different packing modes
Besides aforementioned two typical packing modes where the structure of the module is stable during compression deformation, the cylindrical battery cells can also be arranged for other packing modes, i.e., 2π/3 < θ < π. Fig 12A shows the mechanical behavior of modules with the same size ("10 × 10") for different packing modes. The load plateau can be observed, which is negligible in the whole compression process. After removing the load plateau, the remainder of the curve is parallel to the curve that describes the module for hexagonal dense packing (the dotted lines represent the remainder of the curves). Evidently, the reason for this phenomenon is that the individual cells in the module overcome the friction to conduct relative motion and tend to form hexagonal dense packing. Thus, the verified EBM model for the hexagonal dense packing mode can be used to predict the mechanical behavior for the asymmetric packing modes, only removing the load plateau stage.
As demonstrated in Fig 12B, when the module fully converts to the hexagonal dense packing mode, assuming the cell was not deformed, the theoretical plateau strain ε Ã p can be calculated through the geometrical relationship as follows: The equation shows that the length of the load plateau increases with the augment of θ, and that is in accord with what is presented in Fig 12A. In fact, the plateau length is always smaller than the theoretical value because of the limit of size and boundary condition. Fig 12C shows the mechanical responses of the module for θ = 5π/6 with variables b and l. It is suggested that the plateau length hardly depends on the size while closely related to packing density. An empirical formula predicting the length of plateau can be proposed as follows: As shown in Fig 12D, since the range of packing density of a module relies on its packing mode, g(ρ) can be fitted as the following: where ρ RUC is the density of RUC calculated from Eq. (9) and the ratio ρ/ρ RUC is chosen to eliminate the effect of packing mode. To validate this model for asymmetric packing modes, two packing forms were chosen with different sizes, modes, and densities, namely, θ = 13π/18, where b = 8 and l = 6, and θ = 8π/9, where b = 7 and l = 8. Fig 13 shows the mechanical responses neglecting the load plateau of the module, and it is observed that there is a satisfactory agreement with the EBM model for hexagonal packing mode under the same size (where r 2 = 0.93-0.98), suggesting the accuracy and capability of the EBM model in predicting the mechanical behavior of the battery module.

Concluding remarks
Correctly predicting the mechanical behavior of the module subjected to mechanical loading is critical to understanding the mechanical integrity of the battery pack and ensuring full vehicle crash safety. In this paper, a general method to establish the computational homogenized EBM model is proposed. For the convenience of study, we uniformly defined the packing size and mode of the module. Moreover, the RUC was extracted from the module and the periodic boundary condition was applied on it to acquire the homogenized mechanical properties. The experiments and corresponding DBM model were performed to validate the proposed homogenized EBM model under constrained compression loading conditions. For the cubic dense packing mode, the EBM model can well predict the mechanical behavior with length scale up to 15 × 15 and 15% deformation. For the hexagonal dense packing mode, the revised EBM model based on the packing density can represent well the mechanical properties of the module in the range of 12% deformation under different sizes. For asymmetrical packing modes where the structure of the module is unstable, after removing the load plateau stage that was caused by the relative motion of individual cells, the mechanical behavior of the module can be predicted well by the EBM model for hexagonal dense packing with the same size. Therefore, the homogenized modeling method is widely applicable for different packing modes.