Modeling of deformation of battery cells using thick shell element formulation

Available online at www.sciencedirect.com

ScienceDirect Comput. Methods Appl. Mech. Engrg. 362 (2020) 112840 www.elsevier.com/locate/cma

Modeling of deformation of battery cells using thick shell element formulation Srdjan Simunovica ,∗, Lee P. Bindemanb , Abhishek Kumarc b

a Oak Ridge National Laboratory, P.O. Box 2008-6164, Oak Ridge, TN 37831-6164, United States of America Livermore Software Technology, an ANSYS company, 7374 Las Positas Road, Livermore, CA 94551, United States of America c Northeastern University, 360 Huntington Avenue, Boston, MA 02115, United States of America

Received 24 June 2019; received in revised form 3 December 2019; accepted 8 January 2020 Available online xxxx

Abstract We describe a new approach for modeling nonlinear deformation and stress distribution of battery cells using a new thick shell finite element formulation with a through-thickness calculation of stresses and strains that satisfy equilibrium conditions. Battery cells are transversely layered materials that contain numerous thin layers in a repeating sequence. The layers are made of materials with significantly different mechanical responses. Explicit discretization of the layers is computationally impractical except for very small domains, while homogenized material models cannot account for stress and strain variations and partition through the thickness. The new formulation allows for calculation of the transverse and interlayer stresses based on the material properties of the individual cell layers. The application to problems where the transverse stresses have strong influence is also possible. Published by Elsevier B.V. Keywords: Thick shell finite element formulation; Mechanics of batteries; Transverse deformation of layered materials

1. Introduction The new automotive Li-ion batteries store large amounts of energy and power and need to be protected from mechanical impact similarly to the gas tanks in the internal combustion vehicles. The battery’s primary protection components are made out of structural materials for which sufficiently accurate mechanical models exist. However, once those barriers are breached, the battery’s internal electro-chemical assembly, called jellyroll, can be exposed to the mechanical loads. The mechanical response and models of battery jellyrolls are still under development. The scope of our research is modeling of the mechanical response of the jellyroll under external or internal mechanical loads. The main practical concern with jellyroll deformation is an onset of an internal electrical short, i.e. the mechanical and electrical contact of electrodes of opposite polarities [1], that can result in uncontrollable energy discharge, large electrical currents, heating, and runaway chemical reactions that can ignite the battery [2,3]. The configuration of the internal short will influence its severity, and is a result of the mechanical response of jellyroll during loading. The ∗ Corresponding author.

E-mail addresses: [email protected] (S. Simunovic), [email protected] (L.P. Bindeman), [email protected] (A. Kumar). https://doi.org/10.1016/j.cma.2020.112840 0045-7825/Published by Elsevier B.V.

2

S. Simunovic, L.P. Bindeman and A. Kumar / Computer Methods in Applied Mechanics and Engineering 362 (2020) 112840

Fig. 1. Battery cell (a) Pouch cell, and (b) electrode stack inside the pouch cell.

cell jellyrolls are assemblies of large number of repeated, thin layers of positive electrodes (cathodes), separators, and negative electrodes (anodes) that can be folded or wound in different patterns. Electrodes are assemblies of dissimilar materials, as well. They consist of a metal (e.g. Al, Cu) foil (current collector) coated (usually on both sides) with electro-chemically active material that occupy the largest volume in the cell. A schematic of a pouch cell Li-ion battery and its jellyroll stack are shown in Fig. 1. Thin black and gray layers depict current collectors of opposite polarities, single and cross hatched layers denote opposite polarity active materials, separated by polymer separators depicted by white layers. The jellyroll segment marked by the curly brace in Fig. 1(b) denotes the repeating sequence of layers, usually repeated up to 20 times in the automotive pouch cells. The characteristic thicknesses of the foils, separators, and active layers range from a few to tens of micrometers, whereas the in-plane dimensions can be tens of centimeters or more. Due to this length-scale disparity, modeling of jellyrolls by finite element discretizations that resolve the geometry of each of the layers leads to large models that exceed practical computing capabilities for all but very simple configurations. In order to avoid the computational limitations of explicitly discretized layers, a common modeling approach is to treat the entire jellyroll as a homogeneous material. A prototype constitutive model is selected and the model parameters are set to experiments [4–6]. The homogenized models are very effective for engineering design, but are not readily transferable between different jellyroll types and cannot model deformation and stress partitioning in the layers. In addition, under external loading, models give maximum stresses on the surface and cannot account for subsurface failures and strain partitioning in the cell layers. We are proposing an alternative approach to modeling mechanics of the repeated layered sequence of jellyrolls that is based on a combination of finite element formulation, stress and strain integration for explicit time incrementation, and constitutive material models. The approach takes into account the large variations in the mechanical response and properties between different jellyroll layers. The proposed method does not currently include models for failure of the jellyroll assembly although it can include damage in the layers based on constitutive models. Our approach up-scales deformation in the jellyroll using a finite element formulation and avoids modeling each layer of cell with separate finite elements and their respective materials. The main idea is to incorporate into the model formulation the layered arrangement of the jellyroll, its materials, and interfaces, and embed it into a new element formulation. Mechanics of the battery jellyroll is described in Section 2. Formulation for the new element is given in Section 3. In Section 4, we have compared the stresses and strains obtained from the new layered solid element to benchmark problems discretized by assemblies of solid elements. Apart from faster computational time, embedding multiple cell layers in one element allows for easier modeling of cooperative faults, and interfacial fault between layers. The proposed approach is also applicable to thick shell problems where through-thickness distribution of transverse axial and shear stresses is important. 2. Mechanics of battery jellyroll Battery cell jellyroll is a loosely bonded laminate composite, with alternating layers of electrodes and separators of different mechanical properties. Deformation and stresses of electrode and separator layers are influenced by their stacking sequence, loading and the packaging constraints. Further down, the active materials are porous, lightly bonded aggregates [7–9] saturated with electrolyte; the current collectors are thin metal foils; and separator is

S. Simunovic, L.P. Bindeman and A. Kumar / Computer Methods in Applied Mechanics and Engineering 362 (2020) 112840

3

Fig. 2. Thick shell element form 5 with n integration points.

usually a thin porous polymer membrane. The difference between the mechanical material properties of the layers can introduce significant deformation variations in the thickness direction [10,11] of the jellyroll. The bond between the active material particles and current collector foils is weaker than the bond between the active particles [12], so that the shear forces transmitted from layer to layer are assumed to be small [13]. Metal foils are very thin and can cooperatively buckle and progressively fold [4,14] under in-plane compression. At large, non-uniform, indentations in the cell’s through-thickness direction, metal foils experience membrane tensile forces that are limited by their inherent brittleness [15]. Under an out-of-plane indentation, which is the focus of this study, the jellyroll is primarily deforming by the through-thickness compression and shear in the active layers. The resulting bending creates in-plane tension in metal and polymer layers, and their in-plane compression is limited by the stability of the assembly and the shear strength of the active material [16]. The developed computational model of cell jellyroll is focused on the modeling of through thickness compression problems. The extension of the model to more complex loading situations is underway. 3. Layered thick shell element formulation Formulation for the proposed eight node layered thick shell is presented in this section. The element has been implemented for both implicit static and implicit dynamic solvers, and also for explicit dynamic finite element solutions. Highly nonlinear and large deformation solutions are possible with explicit dynamics and this element is designed to work well in an explicit dynamic code. Because explicit dynamic solutions are conditionally stable on a small solution time step, a computationally efficient element is needed for the element to be useful for large problems. The main idea of the proposed approach was to tailor a finite element formulation to the characteristics of the jellyroll deformation. We wanted to avoid an increase in the number of degrees of freedom that would be needed if the layers were explicitly modeled by nodes or special shape functions. As noted above, the jellyroll is a loosely bonded laminate so that the laminated composite theories can be employed. A review of the various composite laminate finite element formulations can be found in standard textbooks [17] and review papers [18]. The developed layered 8 node brick element shown in Fig. 2 is implemented as thick shell formulation 5 in LS-DYNA [19]. Four nodes define the bottom surface and another four nodes define the top element surface. For computational efficiency, each layer has one in-plane integration point. At least 2 layers are needed through the thickness, but there is no limit to the number of layers that may be defined. The through-thickness orientation of the integration points gives the element the ability to capture bending like a plate or shell element, but it can also capture the thickness stress like a standard brick element. Like the standard brick element (LS-DYNA formulation 1), the layered element is based on the trilinear shape functions. Before describing how the strain is (evaluated, ) (it is useful ) to introduce three coordinate systems, (a) spatial ˆ yˆ , zˆ ≡ xˆ1 , xˆ2 , xˆ3 , and (c) referential (ξ, η, ζ ) ≡ (ξ1 , ξ2 , ξ3 ). The map (x, y, z) ≡ (x1 , x2 , x3 ), (b) corotational x, to the corotational system is a simple rotation from the global spatial coordinate system. The corotational system has its origin at the element centroid, and is defined by vectors that pass through the centers of opposite faces. For skewed elements, the corotational system is orthogonalized and may not pass exactly through face centers. In the equations that follow, a hat on a variable or array indicates that it is evaluated using this corotational system. In the referential system, the element is a cube with edge lengths of 2. The origin of the referential system is at the center and the axes intersect the element faces at −1 and +1. Fig. 3 shows the coordinate system.

4

S. Simunovic, L.P. Bindeman and A. Kumar / Computer Methods in Applied Mechanics and Engineering 362 (2020) 112840

Fig. 3. Element geometry in global, corotational, and referential coordinate systems.

In the equations that follow, the uppercase subscript, I , refers to node numbers 1 to 8. Lower case subscripts refer to the spatial or corotational coordinates, 1 to 3. Greek subscripts indicate range 1 to 4. Repeated subscripts in expressions imply summation over the range of the subscripts unless otherwise noted. Lower case subscripts following comma denote partial derivatives with respect to the corotational coordinate. Carets are added over variables to indicate they have been rotated to the corotational coordinate system or that they have been evaluated from nodal coordinates or velocities that have been rotated to the corotational coordinate system. Boldface characters indicate arrays which are typically dimensioned 8 × 1 which contain one value for each element node. The well-known trilinear shape functions N I (ξ, η, ζ, ) shown in Eq. (1) can be written in terms of the referential coordinates, and the value of the referential coordinates at the nodes, either +1 or -1. 1 (1) N I (ξ, η, ζ ) = (1 + ξ I ξ ) (1 + η I η) (1 + ζ I ζ ) 8 In the referential system, we can evaluate a variable at any location (ξ, η, ζ ) within the element in terms of nodal values, u i I , using the trilinear shape functions from Eq. (1) as shown in Eq. (2). u i (ξ, η, ζ ) = N I (ξ, η, ζ ) u i I

(2)

The development of a strain field follows the development for the brick element with 1 integration point by Belytschko and Bindeman [20], where it was shown that the displacement or velocity field could be written in ˆ orthogonal hourglass shape vectors γˆ of Belytschko, terms of the spatial derivatives of the shape functions, b, et al. [21] and the nodal values of displacements or velocities vˆ , as: ) ( T vˆi = vˆ0i + xˆ j bˆ j + h α γˆ αT · vˆ i , (3) ∫ ∫ 1 ∂ N (ξ, η, ζ ) ∂ N I (ξ, η, ζ ) 1 dΩ , bˆi I = dΩ , (4) bˆ i = V Ωe ∂ xˆi V Ωe ∂ xˆi ( ) ] 1[ γˆ α = hα − hαT xˆ j bˆ j . (5) 8 In Eqs. (3) to (5), xˆ i and vˆ i denote vectors of the nodal coordinates and velocity directional components i: = (xi1 , xi2 , xi3 , xi4 , xi5 , xi6 , xi7 , xi8 ) = (vi1 , vi2 , vi3 , vi4 , vi5 , vi6 , vi7 , vi8 )

xˆ iT vˆ iT

, .

(6)

T

The bˆ i and γˆ αT terms in Eqs. (3) to (5) are transposed 8 × 1 arrays so that a dot product with the 8 × 1 array of nodal velocities vˆ i produces a scalar component i of the velocity. The symbol vˆ0i represents velocity component in xˆi direction of the element center which drops out of the element strain equations as strain is not a function of the rigid body motion. In Eq. (3), h 1 to h 4 are defined in terms of the referential coordinates as shown in Eq. (7). h 1 ≡ ηζ

,

h2 ≡ ζ ξ

,

h3 ≡ ξ η

,

h 4 ≡ ξ ηζ

(7)

Evaluating expressions for h 1 to h 4 at the 8 nodes of an element creates four hα , 8 × 1 arrays of +1 and -1: h1T h2T h3T h4T

≡ ≡ ≡ ≡

(+1, +1, −1, −1, −1, −1, +1, +1) (+1, −1, −1, +1, −1, +1, +1, −1) (+1, −1, +1, −1, +1, −1, +1, −1) (−1, +1, −1, +1, +1, −1, +1, −1)

, , , .

(8)

S. Simunovic, L.P. Bindeman and A. Kumar / Computer Methods in Applied Mechanics and Engineering 362 (2020) 112840

5

The four arrays in Eq. (8) are the base vectors in the referential system for the hourglass modes of deformation which are neglected by one-point quadrature in hexahedral element. In Eq. (4), V is the element volume, and the integral is over the element domain, Ωe . Closed form expressions of bˆ in terms of nodal coordinates are presented in Belytschko et al. [21]. Using Eq. (3) with nodal velocities, vˆ , and taking the derivative with respect to the corotational coordinates, we ˆ as get terms of the spatial gradient of velocity in the corotational system, L, ) ( T ∂ vˆi ∂ NI Lˆ i j = = vˆi I = bˆ j + h α, j γˆ αT · vˆ i . (9) ∂ xˆ j ∂ xˆ j ˆ which is the symmetric part The strain measure of the element will be calculated from the rate-of-deformation, D of Lˆ in Eq. (9). If the 6 unique terms of the rate-of-deformation are arranged in a 6 × 1 array, then the calculation of them can be written as ⎡ ⎫ ⎧ ˆ 11 ⎪ ⎪ D ⎪ ⎪ ⎪ ⎢ ⎪ ⎢ ⎪ ˆ 22 ⎪ ⎪ ⎪ D ⎢ ⎪ ⎪ ⎪ ⎪ ⎬ ⎢ ⎨ ˆ D33 ⎢ =⎢ ˆ ⎢ ⎪ ⎪ 2 D 12 ⎪ ⎪ ⎢ ⎪ ⎪ ⎪ ⎪ ˆ ⎢ ⎪ ⎪ 2 D ⎪ 13 ⎪ ⎪ ⎪ ⎣ ⎭ ⎩ ˆ 2 D23

a matrix multiplication as T bˆ 1 + h α,1 γˆ αT

0

0

⎤

0

T bˆ 2 + h α,2 γˆ αT

0

0

0

T bˆ 2 + h α,2 γˆ αT T bˆ 3 + h α,3 γˆ αT

T bˆ 1 + h α,1 γˆ αT

⎥ ⎥ ⎧ ⎫ ⎥ ⎥ ⎨ vˆ 1 ⎬ ⎥ vˆ 2 ⎥· ⎭ ⎥ ⎩ vˆ 3 ⎥ ⎥ ⎦

0

0 T bˆ 3 + h α,3 γˆ αT

T bˆ 3 + h α,3 γˆ αT

0 T bˆ 1 + h α,1 γˆ αT T bˆ 2 + h α,2 γˆ αT

.

(10)

In Eq. (10), the repeated α subscripts on h α,i and γˆ αT indicate a summation from 1 to 4. Index i after comma in t subscript h α,i denotes partial derivative with respect to corotational axis xˆi . Calling the matrix with the bˆ i and γˆ αT ˆ we can write Eq. (10) as a matrix–vector multiplication terms B, ˆ =B ˆ · vˆ D

.

(11)

Eq. (11) gives the rate of deformation that is consistent with the trilinear shape functions. If this was used as a strain ˆ measure for the element, it would exhibit volumetric locking and excessive shear stiffness. Modifications to the B matrix were made in Belytschko and Bindeman [20] to correct these problems in a 1-point integration hexahedral element based on the QBI (Quintessential Bending Incompressible) assumed strain field [22]. The thick shell also uses the QBI assumptions, not only for hourglass control, but also to modify the strain rate at the integration points to achieve element stress with the same beneficial properties. Compared to the 1-point hexahedron, the layered thick shell element can achieve much better accuracy in coarse mesh bending tests with nonlinear material since the stress is evaluated in each layer. To prevent locking in incompressible materials, the QBI strain field adds terms in the upper half of the matrix to approximate Poisson’s effects. It also removes terms from the bottom half to prevent excessive shear stiffness during bending. This correction to the shear stiffness is consistent with assuming the top and bottom surface of the element is curved during bending. To simplify equations, we introduce notation with up to four subscripts that indicate which terms of the ˆ Yˆ , or Zˆ are used to indicate the variable for the partial derivative. Examples summation are retained, and X, of this notation are shown in Eq. (12) Xˆ 1234 ≡ h α,1 γˆ αT

,

Yˆ 1234 ≡ h α,2 γˆ αT

,

Zˆ 1234 ≡ h α,3 γˆ αT

,

(12)

where summation over index α is assumed. Notation for fewer terms in the sums is Xˆ αβ ≡ h α,1 γˆ αT + h β,1 γˆ βT , Yˆ αβ ≡ h α,2 γˆ αT + h β,2 γˆ βT , Zˆ αβ ≡ h α,3 γˆ αT + h β,3 γˆ βT Xˆ α ≡ h α,1 γˆ αT , Yˆ α ≡ h α,2 γˆ αT , Zˆ α ≡ h α,3 γˆ αT where repeated indices do not imply summation.

(13)

6

S. Simunovic, L.P. Bindeman and A. Kumar / Computer Methods in Applied Mechanics and Engineering 362 (2020) 112840

ˆ matrix for the QBI assumed strain is shown in Eq. (14) Using Eqs. (12) and (13), the B ⎡ ⎤ T T T T T T bˆ 1 + Xˆ 1234 −ν Yˆ 3 − ν Yˆ 24 −ν Zˆ 2 − ν Zˆ 34 ⎢ ⎥ T T T T T ⎥ T ⎢ ⎢ −ν Xˆ 3 − ν Xˆ 14 bˆ 2 + Yˆ 1234 −ν Zˆ 1 − ν Zˆ 34 ⎥ ⎢ ⎥ T ⎢ ⎥ T T T T T ⎢ −ν Xˆ 2 − ν Xˆ 14 −ν Yˆ 1 − ν Yˆ 24 bˆ 3 + Zˆ 1234 ⎥ ⎥ ˆB = ⎢ ⎢ ⎥ ˆ T + Yˆ T ˆ T + Xˆ T ⎢ ⎥ b b 0 2 12 1 12 ⎢ ⎥ ⎢ ⎥ T T T T ˆ + Zˆ ˆ + Xˆ ⎢ ⎥ b 0 b 3 13 1 13 ⎣ ⎦ T T bˆ 3 + Zˆ 23

0

(14)

T T bˆ 2 + Yˆ 23

where ν Poisson’s ratio and νˆ = ν/ (1 − ν). The reasoning for this form of the strain operator is given in [20]. Because the corotational coordinate system and the referential coordinate system are aligned as they both rotate with the element, we can make the simplifying assumptions in Eqs. (15) ∂ξi ∂ξi ≈ 1 = ( 1 ) , = 0 for i ̸= j 1 ∂ xˆi ∂ xˆi ∂ xˆ j . (15) T Λi xˆ i 8 ∂ξi For skewed elements, Eqs. (15) are approximate, but Belytschko and Bindeman [20], concluded that the assumptions do not affect element behavior. In Eq. (15), the lowercase subscripts have values of 1, 2, or 3 and refer to components of the corotational coordinate system. The Λ terms, shown in Eq. (16), are 8 × 1 arrays of +1 and −1 that represent the stretching modes in the referential system. Λ1T ≡ (−1, +1, +1, −1, −1, +1, +1, −1) Λ2T ≡ (−1, −1, +1, +1, −1, −1, +1, +1) Λ3T ≡ (−1, −1, −1, −1, +1, +1, +1, +1)

(16)

With Eq. (15), and the fact that ξ = 0 and η = 0 for all integration points, we can consider the partial derivatives of h 1 to h 4 from Eq. (7). There are some terms that will be zero throughout the element domain because the partial derivatives are zero everywhere. Others will be a function of the ζ coordinate. A third type will be zero at the integration points, but nonzero elsewhere in the element domain. Examples of the 3 types are shown in Eqs. (17) to (19): ∂h 1 ∂η ∂ζ = ζ +η zero throughout the element domain, (17) ∂ xˆ ∂ xˆ ∂ xˆ ∂η ∂ζ ∂h 1 = ζ +η a function of through thickness axis ζ , only, (18) ∂ yˆ ∂ yˆ ∂ yˆ ∂h 1 ∂η ∂ζ = ζ +η a function of only in-plane axis η, only. (19) ∂ zˆ ∂ zˆ ∂ zˆ We will omit the first type shown in Eq. (17) from B because they will not contribute anything. The second type shown in Eq. (18) will be included in the calculation of strain. The third type shown in Eq. (19) will be used in the hourglass control. We can then split B into 2 parts, B in Eq. (21) which will generate strain at the integration points ˜ in Eq. (22) which will be used to construct hourglass control. and is passed to the stress update routine, and B ˆ = B+B ˜ B ⎡ T T bˆ + Xˆ 2 ⎢ 1 ⎢ ⎢ 0 ⎢ ⎢ T ⎢ −ν Xˆ 2 ⎢ B=⎢ T T ⎢ bˆ 2 + Yˆ 1 ⎢ ⎢ T ⎢ bˆ 3 ⎣ 0

(20) 0 T bˆ 2

T + Yˆ 1 T −ν Yˆ 1 T T bˆ 1 + Xˆ 2

0 T bˆ 3

0

⎤

⎥ ⎥ 0 ⎥ ⎥ T ⎥ bˆ 3 ⎥ ⎥ ⎥ 0 ⎥ ⎥ T ⎥ bˆ 1 ⎥ ⎦ ˆbT 2

(21)

S. Simunovic, L.P. Bindeman and A. Kumar / Computer Methods in Applied Mechanics and Engineering 362 (2020) 112840

⎡

T Xˆ 34

⎢ T T ⎢ ⎢ −ν Xˆ 3 − ν Xˆ 4 ⎢ T ⎢ −ν Xˆ 4 ˜ =⎢ B ⎢ ⎢ 0 ⎢ ⎢ T ⎢ Zˆ 1 ⎣ 0

T T −ν Yˆ 3 − ν Yˆ 4

T T −ν Zˆ 2 − ν Zˆ 4

⎤

⎥ T T ⎥ −ν Zˆ 1 − ν Zˆ 4 ⎥ ⎥ T ⎥ Zˆ 124 ⎥ ⎥ ⎥ 0 ⎥ ⎥ T ⎥ ˆ X3 ⎦

T Yˆ 34 T

−ν Yˆ 4 0 0 T Zˆ 2

7

(22)

T Yˆ 3

With B we can evaluate the rate-of-deformation at the integration points which can be used to incrementally update the stress to generate nodal forces, as shown in Eq. (23). ˆ = B · vˆ , D ( ) ˆ σ n = σ (n−1) + S C, D∆t, ω ∫ n T f = Ωe B · σ n dΩ .

,

(23)

In Eq. (23), the superscripts n − 1 and n indicate sequential solution time steps. The function S represents the stress ˆ the solution time step, update which is a function of constitutive properties, C, the current rate of deformation, D, ∆t, and possibly history variables, ω, such as plastic strain and current yield stress. The stress update evaluates an increment of stress and adds this to the stress from the previous solution cycle to obtain a new stress. Integration is over the element domain, Ωe . It is done by evaluating the terms within the integral at the integration points, and using appropriate numerical integration weighting factors. Eq. (23) is used to evaluate element nodal forces at each integration point in the element. For example, an element with 3 layers would evaluate forces at the three Gauss √ √ integration points (ξ, η, ζ ) = (0, 0, 0), (0, 0, −1/ 3), and (0, 0, 1/ 3). ˜ from Eq. (22) we can update the hourglass forces ˜f using a procedure that is similar to Eq. (23) as Using B follows: ˜ · vˆ , q=B n (n−1) Q = Q + C · q∆t , (24) ∫ T n ˜f n = B · Q dΩ . Ωe ˆ and σˆ as in Eq. (23), we have intermediate terms q and Q in Eq. (24) which we will call Instead of evaluating D the hourglass deformation rate and hourglass stress. Q is a function of material properties, C. For many materials, C is the elastic constitutive matrix, but for plasticity models, the average tangent modulus is used. The use of elastic constants is usually sufficiently accurate as the energy in the hourglass mode of the thick shell is typically very small compared to the strain energy, particularly because the modes for bending about the xˆ axis and yˆ axis are included in the stress update. Therefore, plate bending with thick shells does not require stacking elements as it does with bricks. The 2 × 2 × 2 numerical integration in Eq. (24) is done in closed form so no looping over integration points is needed. The element forces calculated by Eqs. (23) and (24) must be transformed from the corotational element coordinate system to the global coordinate system before adding them to a global nodal force array. The Dˆ 13 and Dˆ 23 terms of rate-of-deformation calculated by Eq. (23) are the same at all integration points. It is well known that for a rectangular cross section, the out-of-plane shear stress will have a parabolic distribution with maximum stress at the mid-plane and zero stress at the outer surfaces. To achieve this, Dˆ 13 and Dˆ 23 are scaled for each layer. The scale factors can be continuously recalculated for plasticity models such that proper shear distribution will be maintained as layers of the element reach the yield surface of the material. To this point, the discussion has assumed that the material is homogeneous. In other words, all layers use the same material model with the same properties. For modeling layered composites of fiber reinforced materials, the strong fiber direction may vary by layer. Other composite models may use different materials or material properties in different layers. When material properties vary by layer, force from stress due to in plane deformation, Dˆ 11 , Dˆ 22 , and Dˆ 12 can be obtained simply by using the appropriate constitutive matrix in Eq. (23). The out-of-plane shears, Dˆ 13 and Dˆ 23 and the through-thickness deformation, Dˆ 33 pose a bigger challenge, particularly for Dˆ 33 since

8

S. Simunovic, L.P. Bindeman and A. Kumar / Computer Methods in Applied Mechanics and Engineering 362 (2020) 112840

through-thickness stress is coupled with in-plane stresses due to Poisson’s effect. Scaling Dˆ 13 and Dˆ 23 is unlikely to work well because the relationship between stress and strain is no longer consistent between layers. Therefore, options exist in the current formulation to apply an additional correction such that the shear stress distribution is either is parabolic, constant, or user-defined. Similarly, the constant assumption for the Dˆ 33 term may not be accurate for composite models. Ideally, for rigidly bonded layers, the element formulation should have continuity in both displacements and stresses across the thickness direction. The two main approaches in laminate element formulations are equivalent single layer, and discrete layer theories. The discrete layer theories [23,24] can meet the above continuity requirements but their total number of degrees of freedom is proportional to the number of layers so that the computational performance resembles explicit discretization of layers with solid elements. The equivalent single layer theories assume the functional form of through-thickness deformation and thereby reduce the number of degrees of freedom considerably. The continuity of through-thickness stresses requires discontinuity of throughthickness strains. Therefore, smoothing the displacements by increasing the order of functions does not help. Using deficient displacement approximation is necessary, i.e. through-thickness displacements having selective discontinuities of the first derivative (zig-zag effect). Carrera [25] provides a review of zig-zag theory for modeling through-thickness stress and strain variations in composite materials. A recent example [26] uses zig-zag shape function [27] in the through-thickness direction for a problem of alternating soft-stiff layers. In this case, stiff layers were assumed to have almost identical rotations of their respective normals, whereas the soft layers had different rotations. The zig-zag shape functions add degrees of freedom which complicates the element and time integration formulations. An alternative method for a recovery of the through-thickness normal and shear stresses is the post-processing method [18], also known as two-step method. In the 1st step, the strain increment is predicted using the displacement field from the equivalent single layer theory and in-plane stresses are calculated from the equilibrium equations. In the 2nd step, the stresses are corrected by equilibrium equations of elasticity across thickness. The method does not add degrees of freedom and can be relatively simply added to the existing element theories. The two-step method was first used by Pryor and Barker [28]. In the standard implementation, the calculation of transverse shear stresses from the equilibrium calculations requires higher order approximations of in-plane stresses or displacements and their gradients. Some authors obtain these stress gradients by assuming kinematics (e.g. cylindrical bending) [29,30], Isogeometric Analysis methods [31], or using multi-element patches [32,33]. As noted before, the active materials, which constitute the largest volume fraction, have very low shear modulus, and mechanically behave as weakly bonded granular materials [16]. They deform inelastically almost from the onset of deformation so that the higher order quadrature or displacement functions through the thickness may not significantly improve the accuracy of the stress distribution [34]. In addition, active materials have weak adhesion to the current collector foils [9], so that the interfaces between the layers are not rigid. Complex theories to account for the interface displacement discontinuities [35] are not computationally practical. Experiments show that it is reasonable to assume that the interlaminar shear stresses are low and that the deformation in the through thickness direction is dominated by normal stresses. Given a low cost of the two-step method, and its reasonable accuracy [25], we have implemented it in an incremental form of correction for the through-thickness strains to approximate a constant through-thickness stress distribution. With this correction, soft layers will have larger strain increments than stiff layers. This correction may not only improve the solution in the thickness direction but may also improve the in-plane stress due to Poisson’s effects. These modifications for non-homogeneous materials enable the element formulation to capture some of the complexity in the behavior of layered composites. The use of reduced integration within the layers enables fast explicit solutions which is essential for battery impact simulations as a part of larger structural impact simulation such as vehicle crash. 4. Numerical examples The developed modeling formulation has been compared against the models based on solid element assemblies for different modes of loading. The solid element benchmark models use for discretization of layered structure at least one element through the thickness of each layer. The solid assembly benchmark has many more degrees of freedom and represents the best possible outcome of the new element formulation. We demonstrate the performance of the new element formulation on the problems of through-thickness compression, shear and bending.

S. Simunovic, L.P. Bindeman and A. Kumar / Computer Methods in Applied Mechanics and Engineering 362 (2020) 112840

9

Fig. 4. Transverse uniaxial compression (a) Solid element assembly, and (b) thick shell element. Unit of dimension is mm.

4.1. Transverse uniaxial compression We performed simulations of a through-thickness compression on an assembly of solid blocks with equal dimensions of 0.5 mm in x and y directions. The thickness in millimeters (z-direction) of each block is shown in Fig. 4(a). Each block is represented by one standard constant stress solid element with one integration point. Solid element formulation is selectively reduced integrated element (element formulation −1 in LS-DYNA [19]). Equivalent thick shell element is shown in Fig. 4(b). The layered solid element has 5 integration points (one integration point for each solid block) through the thickness (vertical) direction. For this example, we choose 2 different materials (Material 1 and Material 2) and assigned them to 5 layers as shown in Fig. 4(a). Material models for both these materials are isotropic elasto-plastic model with piecewise-linear isotropic hardening, i.e. MAT-24 in LS-DYNA [19]. Elastic modulus, E, for material 1 and material 2 is 2.0 GPa and 1.0 GPa, respectively. Poisson’s ratio, ν, is 0.2., yield stress, σ y , for material 1 is 50 MPa, while that of material 2 is 25 MPa. Tangent modulus, E T , for material 1 and material 2 are 2 and 1 MPa, respectively. Incremental vertical displacement is imposed on the top nodes of the two models at the rate of 0.03 m/s. Transverse stress and strain distributions through the thickness for the layered-solid element are very similar to that of the solid element assembly counterpart, as shown in Fig. 5. The number of the degrees of freedom in solid element model is 3 times larger. The differences between the two models normalized by the values in the layered thick shell model are given in Figs. 5(e) and 5(f). In the next example, we performed a uniaxial compression on a layered structure as shown in Fig. 6(a). The model represents a single repeating sequence of the layers in a battery jellyroll. The sequence of layers of the battery jellyroll is (Cu-A-S-C-Al-C-S-A), where Cu denotes copper foil current collector, A is anode, S is separator, C is cathode and Al is aluminum foil. Cu, Al and S layers are modeled using isotropic elasto-plastic material models with piecewise linear isotropic hardening (MAT-24 in LS-DYNA [19]). Anode and cathode active layers are modeled using crushable foam material model (MAT-63) [36,19,13]. In the stacked solid element model, Fig. 6(a), one element is used for each layer. The equivalent thick shell element representation is shown in Fig. 6(b). Horizontal dimensions of the cells were x = 0.5 mm and y = 0.5 mm. Vertical stacking dimensions are given in Fig. 6. Material model parameters used for each layer are given in Tables 1 and 2. Symbols E, ν, σ y , E T , σc , d in the tables denote Young’s modulus, Poisson’s ratio, yield stress, tangent modulus, cut-off stress, and damping coefficient, respectively. For active layer material models (crushable foam), the curves describing yield stress as a function of volumetric strain are given in Fig. 7. Mechanical properties for all these components have been taken from [37]. Top nodes of the model are compressed in z-direction with rate of 0.02 m/s. As shown in Fig. 8, stress–strain distribution along the thickness direction for the new thick shell element formulation is in excellent agreement with the stacked solid element model. Computational time for thick shell element is about 18 times shorter compared to its solid element counterpart. Maximum normalized differences between the models are 6.2% for Z strain and 3.9% for Z stress.

10

S. Simunovic, L.P. Bindeman and A. Kumar / Computer Methods in Applied Mechanics and Engineering 362 (2020) 112840

Fig. 5. Stress and strain values for solid element (a, c) and thick shell element (b, d) at different layers. Normalized differences for (e) stress and (f) strain between the models. Table 1 Material parameters for current collectors and separator. Component

E (GPa)

ν

σ y (MPa)

E T (MPa)

Cu foil Separator Al foil

110 0.5 70

0.35 0.20 0.36

210.0 10.0 180.0

1100.0 50.0 700.0

Table 2 Material parameters for active layers. Component

E (MPa)

ν

σc (MPa)

d

Anode Cathode

330.0 330.0

0.0 0.0

10.0 10.0

0.05 0.05

4.2. Shear loading Fig. 9 depicts shear loading configuration for comparison between the benchmark model made with solid elements, Fig. 9(a), and the equivalent thick shell model, Fig. 9(b). For this example, we choose two different materials (material 1 and material 2) and assigned them to five layers as shown in Fig. 9(a). Material models for both these materials are isotropic elasto-plastic models with isotropic piecewise linear hardening (MAT-24). Elastic moduli for material 1 and material 2 are 100.0 GPa and 0.5 GPa, respectively. Yield stresses for material 1 is 50 MPa,

S. Simunovic, L.P. Bindeman and A. Kumar / Computer Methods in Applied Mechanics and Engineering 362 (2020) 112840

11

Fig. 6. Transverse uniaxial compression of a single layer of battery cell (a) Solid element assembly, and (b) thick shell element. Unit of dimension is mm.

Fig. 7. Yield stress vs volumetric strain curve for active material taken from [37].

while that of material 2 is 1 MPa. Tangent moduli for Material 1 and Material 2 are 2 and 1 MPa, respectively. Displacements for the top nodes are imposed in x direction at rate of 0.03 m/s. As before, the new thick shell element formulation, Fig. 10(b, d), gives very similar stress and strain distributions as the solid element assembly, Fig. 10(a, c) even for a large property contrast between the layers. Maximum normalized differences between the models are 7.37% for ZX strain and 6.87% for ZX stress in the early loading stage and less than 1% afterwards. 4.3. Bending problem A bending of a long beam was modeled with thick shell elements and with solid elements. The displacement of the nodes on two extreme ends of a beam were fixed and displacement loading at the rate of 0.002 m/s in the vertical direction was applied at the mid-span of the beam, as shown in Fig. 11. The span of the beam was 20 mm (x), width 0.5 mm (y), and thickness 0.5 mm (z). In this simulation, we have compared the bending stress (σxx ) calculated using the solid element assembly model against the layered solid element model. Bending stress is calculated across the thickness for the 13th column of elements from the left side. In case of solid element assembly

12

S. Simunovic, L.P. Bindeman and A. Kumar / Computer Methods in Applied Mechanics and Engineering 362 (2020) 112840

Fig. 8. Transverse stresses and strains for solid element model (a, c) and thick shell element (b, d) at different points of a single battery cell.

Fig. 9. Transverse shear (a) Solid element assembly, (b) thick shell element. Unit of dimension is mm.

we have picked 1st, 2nd, 3rd, 4th and 5th element from the top direction in that column. The beam material is a linear elastic material (MAT-01, E = 1.0 GPa, ν = 0.2.). Stress response of the solid element assembly is very similar to that of the thick shell element model, Fig. 12. Maximum normalized difference of X-stress at the end of the simulation is less than 1%. 4.4. Indentation of a battery cell A schematic of the indentation of a cell pouch with spherical indenter is shown in Fig. 13(a). The cell length (x) and width (y) are 40 mm, and its thickness (z) is 6.405 mm. The spherical punch of diameter 25.4 mm was modeled as a geometrical analytical rigid body. The lower surface of the cell rested on a rigid wall. We prescribed rigid body motion for the sphere and controlled its vertical movement as a linear function of time such that the velocity of sphere in z direction is constant.

S. Simunovic, L.P. Bindeman and A. Kumar / Computer Methods in Applied Mechanics and Engineering 362 (2020) 112840

13

Fig. 10. Transverse shear stresses and strains for solid element model (a, c) and thick shell element (b, d) at different layers.

Fig. 11. Bending test using (a) solid elements, and using (b) thick shell elements.

The pouch cell consists of 17 jellyroll layers along the thickness (z) direction. Each jellyroll layer is composed of one positive and negative current collector and two cathode anode and separator layers as shown in Fig. 6(a). The material model used to represent each component along with thickness and material property of the component is shown in Tables 1 and 2. Mechanical properties for all these components have been taken from [37]. The pouch cover of thickness 0.1 mm at top and bottom of the cell is modeled using an elastic material (MAT-01) with Young’s modulus of 0.25 GPa, and Poisson’s ratio of 0.25. The solid finite element model of the cell indentation is shown in Fig. 13. One solid element was used for every jellyroll layer, as shown in Fig. 4(a). Pouch covers were added on the top and bottom of the jellyroll. The graded

14

S. Simunovic, L.P. Bindeman and A. Kumar / Computer Methods in Applied Mechanics and Engineering 362 (2020) 112840

Fig. 12. Stress distribution through the thickness for (a) solid elements (5 elements) and (b) thick shell element.

Fig. 13. (a) Solid element model of cell indentation, and (b) cell deformation after indentation.

finite element discretization was necessary to keep the element aspect ratios in the indentation region in a reasonable interval. The deformed shape of the cell using the solid element model is shown in Fig. 13(b). The corresponding thick shell element model is shown in Fig. 14. Integration points are stacked as in Fig. 4(b). The pouch covers were added as integration points to the top and bottom layer of elements, respectively. The total number of elements in the solid element model is 76,728, and in the thick shell model is 1700. The calculation time for the layered solid element is 55 times faster, not accounting for the failure of the solid model at the 2/3 of the loading. Fig. 15 shows the load displacement curve obtained from the indentation test simulation. The solid element model breaks down due to a negative element Jacobian. Such a problem does not arise in the thick shell element solution despite a significantly faster computational time. Coarsening of the thick shell model by increasing the number of layers per element for a given in-plane discretization gives essentially the same indentation force response. The deformed configurations of the models with 8 and 1 thick shell elements through thickness of the cell, are shown in Figs. 16(a) and 16(b), respectively. The corresponding displacement–force data is shown in Fig. 17(a). The effect of mesh refinement in plane of the cell is shown in Fig. 17(b). The difference between the two in-plane discretizations is attributed to a better accommodation of the sphere contact by the finer in-plane mesh. Computational times for different through thickness discretizations are shown in Table 3. In all cases, element processing accounts for 97% of the computational time, The speedup is primarily due to the increased stable time step with increasing element thickness since the number of integration points in all the three models is the same.

S. Simunovic, L.P. Bindeman and A. Kumar / Computer Methods in Applied Mechanics and Engineering 362 (2020) 112840

15

Fig. 14. (a) Thick shell element model of cell indentation, and (b) cell deformation after indentation.

Fig. 15. Indentation force obtained using (a) solid element, and (b) thick shell element model.

Fig. 16. Thick shell element model of cell indentation: (a) 800 elements, and (b) 100 elements.

5. Conclusions We have developed a new version of a thick shell element that accurately calculates through-thickness transverse stresses and strains. The element gives an excellent agreement with the standard solid element benchmarks in compression, shear and bending. The method is orders of magnitude faster and capable of modeling large

16

S. Simunovic, L.P. Bindeman and A. Kumar / Computer Methods in Applied Mechanics and Engineering 362 (2020) 112840

Fig. 17. Indentation force for element mesh refinement (a) through thickness, (b) in plane of the cell. Table 3 Processing time in seconds for different layered thick shell element discretizations. Elements

Total time

Element processing

10 × 10 × 1 10 × 10 × 8 10 × 10 × 17

65.8 529.9 1151.6

58.4 509.4 1108.9

deformations. In the present study, we have used different material models that can be applicable to model the components of a battery cell. The use of thick shell elements offers a possibility to model the response in different components of a battery cell under mechanical loading with acceptable CPU time, and the approach can be scaled to module and pack levels. We will also incorporate failure theories to model the onset of the formation of an internal short circuit in future work. Acknowledgments This work was supported by the Vehicle Technologies Office, Office of Energy Efficiency and Renewable Energy, U.S. Department of Energy under DOE Agreement DE-EE0007288. This manuscript has been authored by UTBattelle, LLC under Contract No. DE-AC05-00OR22725 with the U.S. Department of Energy. The United States Government retains and the publisher, by accepting the article for publication, acknowledges that the United States Government retains a non-exclusive, paid-up, irrevocable, world-wide license to publish or reproduce the published form of this manuscript, or allow others to do so, for United States Government purposes. The Department of Energy will provide public access to these results of federally sponsored research in accordance with the DOE Public Access Plan (http://energy.gov/downloads/doe-public-access-plan). Authors acknowledge discussions with Chulheung Bae, Jie Deng, and Theodore Miller from Ford Motor Company. References [1] P. Ramadass, W.F. Fang, Z.M. Zhang, Study of internal short in a li-ion cell i. Test method development using infra-red imaging technique, J. Power Sources 248 (2014) 769–776, http://dx.doi.org/10.1016/j.jpowsour.2013.09.145. [2] A. Thaler, D. Watzenig, Automotive Battery Technology, Springer, 2014. [3] Q.S. Wang, P. Ping, X.J. Zhao, G.Q. Chu, J.H. Sun, C.H. Chen, Thermal runaway caused fire and explosion of lithium ion battery, J. Power Sources 208 (2012) 210–224, http://dx.doi.org/10.1016/j.jpowsour.2012.02.038. [4] E. Sahraei, R. Hill, T. Wierzbicki, Calibration and finite element simulation of pouch lithium-ion batteries for mechanical integrity, J. Power Sources 201 (2012) 307–321, http://dx.doi.org/10.1016/j.jpowsour.2011.10.094. [5] S. Kim, Y.S. Lee, H.S. Lee, H.L. Jin, A study on the behavior of a cylindrical type Li-ion secondary battery under abnormal conditions, Mater.wiss. Werkst.tech. 41 (5) (2010) 378–385, http://dx.doi.org/10.1002/mawe.201000612. [6] I. Avdeev, M. Gilaki, Structural analysis and experimental characterization of cylindrical lithium-ion battery cells subject to lateral impact, J. Power Sources 271 (2014) 382–391, http://dx.doi.org/10.1016/j.jpowsour.2014.08.014.

S. Simunovic, L.P. Bindeman and A. Kumar / Computer Methods in Applied Mechanics and Engineering 362 (2020) 112840

17

[7] M. Ebner, F. Geldmacher, F. Marone, M. Stampanoni, V. Wood, X-ray.tomography.of. porous, X-ray tomography of porous transition metal oxide based lithium ion battery electrodes, Adv. Energy Mater. 3 (7) (2013) 845–850, http://dx.doi.org/10.1002/aenm.201200932. [8] A.J. Stershic, S. Simunovic, J. Nanda, Modeling the evolution of lithium-ion particle contact distributions using a fabric tensor approach, J. Power Sources 297 (2015) 540–550, http://dx.doi.org/10.1016/j.jpowsour.2015.07.088. [9] T. Woehrle, Lithium-Ion Batteries: Basics and Applications, Springer, Berlin, Heidelberg, 2018, http://dx.doi.org/10.1007/978-3-66253071-9_9. [10] M. Disciuva, Development of an anisotropic, multilayered, shear-deformable rectangular plate element, Comput. Struct. 21 (4) (1985) 789–796, http://dx.doi.org/10.1016/0045-7949(85)90155-5. [11] H. Murakami, Laminated composite plate-theory with improved inplane responses, Trans. ASME, J. Appl. Mech. 53 (3) (1986) 661–666. [12] J.C. Chen, J.Y. Liu, Y. Qi, T. Sun, X.D. Li, Unveiling the roles of binder in the mechanical integrity of electrodes for lithium-ion batteries, J. Electrochem. Soc. 160 (9) (2013) A1502–A1509, http://dx.doi.org/10.1149/2.088309jes. [13] T. Wierzbicki, E. Sahraei, Homogenized mechanical properties for the jellyroll of cylindrical lithium-ion cells, J. Power Sources 241 (2013) 467–476, http://dx.doi.org/10.1016/j.jpowsour.2013.04.135. [14] M.Y. Ali, W.J. Lai, J. Pan, Computational models for simulations of lithium-ion battery cells under constrained compression tests, J. Power Sources 242 (2013) 325–340, http://dx.doi.org/10.1016/j.jpowsour.2013.05.022. [15] H. Wang, T.R. Watkins, S. Simunovic, P.R. Bingham, S. Allu, J.A. Turner, Fragmentation of copper current collectors in li-ion batteries during spherical indentation, J. Power Sources 364 (2017) 432–436, http://dx.doi.org/10.1016/j.jpowsour.2017.08.068. [16] H. Wang, S. Simunovic, H. Maleki, J. Howard, J. Hallmark, Internal configuration of prismatic lithium-ion cells at the onset of mechanically induced short circuit, J. Power Sources (2016) 424–430. [17] J. Reddy, Mechanics of Laminated Composite Plates and Shells, CRC Press, 2004. [18] I. Kreja, A literature review on computational models for laminated composite and sandwich panels, Central Eur. J. Eng. 1 (1) (2011) 59–80. [19] J.O. Hallquist, LS-DYNA Keyword User’s Manual, Livermore Software Technology Corporation, Livermore, CA, 2014. [20] T. Belytschko, L.P. Bindeman, Assumed strain stabilization of the 8 node hexahedral element, Comput. Methods Appl. Mech. Engrg. 105 (2) (1993) 225–260. [21] T. Belytschko, J.S.J. Ong, W.K. Liu, J.M. Kennedy, Hourglass control in linear and nonlinear problems, Comput. Methods Appl. Mech. Engrg. 43 (3) (1984) 251–276, http://dx.doi.org/10.1016/0045-7825(84)90067-7. [22] T. Belytschko, W.E. Bachrach, Efficient implementation of quadrilaterals with high coarse-mesh accuracy, Comput. Methods Appl. Mech. Engrg. 54 (3) (1986) 279–301, http://dx.doi.org/10.1016/0045-7825(86)90107-6. [23] E.J. Barbero, J.N. Reddy, J. Teply, An accurate determination of stresses in thick laminates using a generalized plate-theory, Internat. J. Numer. Methods Engrg. 29 (1) (1990) 1–14, http://dx.doi.org/10.1002/nme.1620290103. [24] Y.M. Ghugal, R.P. Shimpi, A review of refined shear deformation theories of isotropic and anisotropic laminated plates, J. Reinf. Plast. Compos. 21 (9) (2002) 775–813, http://dx.doi.org/10.1106/073168402025748. [25] E. Carrera, A priori vs a posteriori evaluation of transverse stresses in multilayered orthotropic plates, Compos. Struct. 48 (4) (2000) 245–260, http://dx.doi.org/10.1016/s0263-8223(99)00112-9. [26] Y. Liang, B.A. Izzuddin, Nonlinear analysis of laminated shells with alternating stiff/soft lay-up, Compos. Struct. 133 (2015) 1220–1236, http://dx.doi.org/10.1016/j.compstruct.2015.08.043. [27] A. Toledano, H. Murakami, A composite plate-theory for arbitrary laminate configurations, Trans. ASME, J. Appl. Mech. 54 (1) (1987) 181–189. [28] C.W. Pryor, R.M. Barker, Finite-element analysis including transverse shear effects for applications to laminated plates, AIAA J. 9 (5) (1971) 912, http://dx.doi.org/10.2514/3.6295. [29] R. Rolfes, A.K. Noor, K. Rohwer, Efficient calculation of transverse stresses in composite plates, in: 1997 MSC Aerospace User’s Conference, Newport Beach, CA, 1997, pp. 1–18. [30] R. Rolfes, K. Rohwer, M. Ballerstaedt, Efficient linear transverse normal stress analysis of layered composite plates, Comput. Struct. 68 (6) (1998) 643–652, http://dx.doi.org/10.1016/s0045-7949(98)00097-2. [31] J.E. Dufour, P. Antolin, G. Sangalli, F. Auricchio, A. Reali, A cost-effective isogeometric approach for composite plates based on a stress recovery procedure, Composites B 138 (2018) 12–18, http://dx.doi.org/10.1016/j.compositesb.2017.11.026. [32] O.C. Zienkiewicz, J.Z. Zhu, The superconvergent patch recovery and a-posteriori error-estimates. 1. The recovery technique, Internat. J. Numer. Methods Engrg. 33 (7) (1992) 1331–1364, http://dx.doi.org/10.1002/nme.1620330702. [33] F. Daghia, S. de Miranda, F. Ubertini, Patch based recovery in finite element elastoplastic analysis, Comput. Mech. 52 (4) (2013) 827–836, http://dx.doi.org/10.1007/s00466-013-0847-6. [34] T. Belytschko, W.K. Liu, B. Moran, Nonlinear Finite Elements for Continua and Structures, Wiley, New York, 2000. [35] D. Liu, L. Xu, X. Lu, An interlaminar bonding theory for delamination and nonrigid interface analysis, J. Reinf. Plast. Compos. 12 (11) (1993) 1198–1211, http://dx.doi.org/10.1177/073168449301201105. [36] M. Shaw, T. Sata, The plastic behavior of cellular materials, Int. J. Mech. Sci. 8 (7) (1966) 469–478, http://dx.doi.org/10.1016/00207403(66)90019-1. [37] C. Zhang, S. Santhanagopalan, M.A. Sprague, A.A. Pesaran, A representative-sandwich model for simultaneously coupled mechanicalelectrical-thermal simulation of a lithium-ion cell under quasi-static indentation tests, J. Power Sources 298 (2015) 309–321, http://dx.doi.org/10.1016/j.jpowsour.2015.08.049.