3D Isogeometric Analysis of intercalation-induced stresses in Li-ion battery electrode particles

Comput. Methods Appl. Mech. Engrg. 268 (2014) 225–244

Contents lists available at ScienceDirect

Comput. Methods Appl. Mech. Engrg. journal homepage: www.elsevier.com/locate/cma

3D Isogeometric Analysis of intercalation-induced stresses in Li-ion battery electrode particles P. Stein ⇑, B. Xu Mechanics of Functional Materials Division, Institute of Material Science, Technische Universität Darmstadt, Petersenstrasse 32, 64287 Darmstadt, Germany

a r t i c l e

i n f o

Article history: Received 23 July 2013 Received in revised form 4 September 2013 Accepted 22 September 2013 Available online 3 October 2013 Keywords: Isogeometric Analysis Lithium-ion batteries Intercalation Coupled diffusion Drifting

a b s t r a c t This paper is concerned with a three-dimensional coupled chemo-mechanical model for intercalation-induced stresses. Thereby, a diffusion model is enhanced by a drifting term involving the gradient of the hydrostatic stress state. A straightforward Finite Element discretization of the coupled diffusion would require C 1 -continuous shape functions. The implementation is thus based on the concept of Isogeometric Analysis, permitting a discretization purely in terms of the displacements and the concentration field. Furthermore, it allows to set up an operator matrix for the computation of the hydrostatic stress gradient in terms of the primal variables. The model is verified, in a simplified form, using analytical results from literature. It is subsequently employed for the study of the mechanical behavior of spherical and ellipsoidal particles under galvanostatic boundary conditions. The simulation data show a stress relaxation effect that becomes significant both in stiff electrode materials and for high charge rates. During charging, a characteristic tensile core-compressive shell structure develops. Once the stresses have reached a material-dependent threshold, they enhance the diffusion and drive ions away from regions of high compressive hydrostatic stress, that is, towards a particle’s core. Over time, this leads to reduced concentration gradients which, in turn, reduces the stresses in the particle already during the charging process. Models that neglect the stress effect consequently predict stress levels that are not only of higher magnitude, but that also do not decay. A comprehensive study on the influence of geometry, material constants, and charge rate has been performed. It sheds light on the distributions of stresses in ellipsoidal electrode particles. These have been reported to exhibit lower stress levels than spherical particles. Using a transformation of Cartesian stresses into a prolate spheroidal coordinate system, it can be shown that these particles equilibrate stresses by deformations along the semimajor axis. At the same time, a region of high stresses develops as a ‘‘belt’’ around the particles’ equator. Their intensity depends on the shape of the particle but is below that observed in spherical and oblate spheroidal particles. This offers an explanation for the longevity of these particles under cyclic charging processes. Ó 2013 Elsevier B.V. All rights reserved.

1. Introduction Lithium-ion batteries have become a key component of modern lifestyle: they provide a compact, lightweight energy source for all kinds of mobile electronic devices, for instance notebooks or smart phones. The batteries offer a high energy

⇑ Corresponding author. Tel.: +49 61511675681. E-mail addresses: [email protected] (P. Stein), [email protected] (B. Xu). 0045-7825/$ - see front matter Ó 2013 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.cma.2013.09.011

226

P. Stein, B. Xu / Comput. Methods Appl. Mech. Engrg. 268 (2014) 225–244

density, show low self-discharge, allow for high operating voltages, and do not show a memory effect [1]. Their cells exhibit a simple structure that is illustrated in Fig. 1. Every cell consists of two electrodes, each of which is connected to a current collector. The cathode is commonly made of a Lithium compound such as LiFePO4, LiMnO2, or LiCoO2. The anode, on the other hand, consists mostly of graphite, silicon, or tin. Depending on their fabrication, diverse morphologies can be achieved, e.g., powder-structured electrodes, thin films, or brushes of nanowires. Inactive binders and conductive additives enhance both mechanical bonding and electric conductivity. The electrodes are separated by either a solid or a liquid, non-aqueous electrolyte. During charge, Lithium ions are extracted from the cathode and diffuse through the electrolyte towards the anode. Here, they are intercalated into the material, that is, they are inserted into the crystal lattice of the host material, in this case graphite or silicon. When the battery is discharged, the respective electrodes change roles (the former cathode becomes the anode and vice versa) and the electrochemical process of extraction and intercalation is reversed. An overview of the different thermodynamic and electrochemical aspects of the intercalation processes can be gained from [2–4]. Experiments show that, over several hundred charge cycles, the electrodes experience delamination and contact loss, material fatigue, and literal pulverization of the electrode particles. This suggests that high stresses develop in the electrodes. This degradation occurs in both cathode and anode, and it can also be observed for small deformations of the electrode [5]. Being detached from the conductive network, the particle’s fragments become inactive and no longer participate in the intercalation process. As more and more particles turn into passive material the capacity of the cell fades, and hence the overall capacity of the battery. In recent years, there has been a renewed interest in these batteries as they possess a high potential for application in hybrid-electric vehicles [6]. This puts the batteries under extreme operating conditions, namely short charge cycles with a high current density. Consequently, intensive research is necessary in order to increase the batteries’ lifetime, attainable charge rates, and their capacity. Theoretical studies of intercalation batteries focus on two scales of this problem, the first of which being the behavior of a complete battery. A framework for the two-dimensional simulation of the diffusion in the electrolyte and the electrodes has been put forward in [7,8], which highlights the role of the cell microstructure for optimal charge processes. Similar but onedimensional models have been proposed in [9–11], gauging the effects of cell-scale parameters on properties such as capacity fade. The second group of studies is directed towards the mechanical behavior of single particles. An important contribution in this regard is the work of Christensen and Newman [12] who have set up an analytical model for the coupled chemomechanical problem in a spherical particle. Their description of principal effects, for instance the emergence of a core–shell structure, motivated a series of studies. Investigated aspects comprise a particle’s behavior under potentiostatic and galvanostatic boundary conditions and the influence of surface mechanical effects on the behavior of a particle [13,14]. Other works considered varying elastic moduli and phase transformations [15–17], plasticity [18], and fracture [5,19,20]. These studies assume smooth, regularly shaped particles in order to reduce the model equations to a one-dimensional problem. Hence, they cannot account for realistic particle and electrode shapes. In particular, they tend to underestimate the stress

Fig. 1. Schematic composition of a Lithium-ion battery.

P. Stein, B. Xu / Comput. Methods Appl. Mech. Engrg. 268 (2014) 225–244

227

concentrations produced by surface irregularities by an order of magnitude [21]. Furthermore, they neglect electrode morphologies such as nanowires or thin films (that show a mechanical behavior distinct from spherical particles [22–24]). A number of groups dealt with the numerical simulation of electrode particles and provided three-dimensional models for diffusion-induced stresses. These simulations have been used to study the stress distribution in spherical, ellipsoidal, and cylindrical particles [1,25,26], as well as the effect of temperature [27], stress generation due to phase transitions [28], and fracture [29–31]. In order to allow verification of results, above studies focused on spherical and ellipsoidal particles. An exception is the work of Seo et al. [21] who obtained realistic particle shapes from atomic force microscopy. Unfortunately, the equations describing diffusion-induced stresses do not allow for a direct implementation in terms of Finite Elements. In particular their convergence is not guaranteed for this problem. Intercalation is modeled as a coupled chemo-mechanical problem in which the stress field contributes a drifting term to the classical Fick diffusion. This term, derived from the chemical potential, depends on the gradient of the hydrostatic stress. Accordingly, the weak form of the coupled diffusion problem possesses a variational index of 2, requiring C 1 -continuous shape functions. Implementations based on C 0 -continuous Lagrange shape functions can consequently show erratic convergence behavior that can be attributed to the violation of the compatibility requirement [32]. In order to overcome these problems, one can employ a mixed-variational formulation such as in a recent work by Gao et al. [33]. There, quadratic shape functions are used for the interpolation of the displacement field, and the hydrostatic stress field has been introduced as a redundant variable, interpolated by linear basis functions. The disadvantage of mixed formulations is that they are not straightforward, as only certain combinations of interpolation schemes lead to stable computational schemes. This paper employs the concept of Isogeometric Analysis (IGA), introduced by Hughes et al. [34], as an alternative approach. It allows to employ a wide range of smooth, higher-order basis functions, e.g., Non-uniform rational B-Splines (NURBS) or T-Splines [35], whose continuity can be adapted both in geometric modeling and during discretization. This proves to be beneficial for a wide range of applications, from structural vibrations [36,37], to shell analysis [38,39], up to fluid flow and turbulence [40–42]. The compact representation of geometries as provided by NURBS and T-Splines allows for efficient model representations [43] and accommodates efficient optimization algorithms [44,45]. Of particular interest for the field of material science is the isogeometric treatment of phase-field models [46–48]. In the context of this work, a single set of shape functions of higher continuity allows for an implementation of the governing equations purely in terms of the displacement and concentration fields. Following this introduction, the mathematical model for coupled diffusion in a bulk electrode particle is described in Section 2. Its implementation using the concept of IGA is given in Section 3. The beginning of Section 4 shows the verification of the model using an analytical model from literature. Thereby, a simplified model is employed, using only a partial coupling of the solution fields. The section subsequently deals with the simulation of spherical and ellipsoidal electrode particles under galvanostatic boundary conditions. It demonstrates the effect of the aforementioned drifting term and contains a comprehensive study of the effect of geometry, material parameters, and charge rate on the resulting stress levels. The paper concludes with a summary and discussion of the results in Section 5.

2. Problem description An electrode particle is described as a three-dimensional domain B with the boundary @ B ¼ S. The evolution of mechanical stresses is given as a coupled chemo-mechanical problem of two fields, namely the displacements and the concentration field. Thereby, the mechanical problem is treated as quasi-static since the characteristic time for atomic migration is much slower than mechanical deformation [49]. Hence, the equilibrium of internal and external forces in the particle is described by

rij;j þ bi ¼ 0 8xi 2 B;

ð1Þ

with rij being the Cauchy stress tensor and bi the vector of body forces. As common, a partial derivative with respect to spatial coordinate xi is denoted by ðÞ;i . Furthermore, the summation convention is employed, i.e., paired indices imply a summation over their range. Following Prussin [50] and Li [51], the influence of the concentration field on the stresses is modeled in analogy to thermal effects:





rij ¼ Cijkl ekl  c  cref ekl

 1 ¼ uk;l þ ul;k ; 2

 X 3

 dkl ; ð2Þ

where Cijkl denotes the components of the fourth-order stiffness tensor, ui the displacement vector, ekl the symmetric, second-order infinitesimal strain tensor, c the current concentration, cref the concentration in the stress-free initial configura the partial molar volume of the intercalated ions in the host material, and dkl the second-order unit tensor. The latter tion, X represents an isotropic expansion/shrinking of the host lattice which, strictly speaking, does not hold for materials with a

228

P. Stein, B. Xu / Comput. Methods Appl. Mech. Engrg. 268 (2014) 225–244

layered structure such as LiCoO2 or graphite [6]. Similarly, the elastic moduli are taken to be independent of the current concentration. The evolution of the Li ion concentration is driven by the mass conservation

c_ þ J i;i ¼ 0 8xi 2 B;

ð3Þ

where c_ denotes the derivative of the concentration field with respect to time, and with J i being the ion flux within a particle. The latter is given by

  X J i ¼ D c;i  c rh;i ; RT

ð4Þ

with the diffusion coefficient D, the gas constant R, the absolute temperature T, and the hydrostatic stress rh , in turn defined by

1 3

rh ¼ rij dij :

ð5Þ

The derivation of (4) from the chemical potential can be found, e.g., in Zhang et al. [25]. Thereby, the classical Fick diffusion is enhanced by a drifting term based on the gradient of the hydrostatic stress field. This causes ions to be driven away not only from regions of high concentration towards lower concentration, but also away from regions of high compressive stresses, corresponding to a compression of the host lattice. The assumption of isotropic diffusion may conflict with experimental observations [52] but it gets its justification from a lack of available material data [6,15,53]. Combination of (3) and (4) yields the governing equation for the coupled concentration field problem:

  X 0 ¼ c_  D c;i  8 xi 2 B c rh;i RT ;i   DX DX ¼ c_  D c;ii þ c;i rh;i þ c rh;ii : RT RT

ð6Þ

The body under consideration is restrained against rigid-body movements by suitable displacement and traction boundary conditions

^i ui ¼ u

8xi 2 S u  S;

rij nj ¼ ^ti

ð7aÞ

8xi 2 S t  S;

ð7bÞ

where S ¼ S u [ S t ; S u \ S t ¼ £ must be fulfilled. Furthermore, the following boundary conditions are prescribed for the concentration field and for the ion flux across the boundary,

c ¼ ^c 8xi 2 S c  S; in J i ni ¼ ^f ¼ F ¼ Dc;i ni þ

ð7cÞ

8xi 2 S J  S;  DX c rh;i ni ; RT

ð7dÞ

in which in refers to the applied electric current density, and with F as Faraday’s constant. As before, the conditions S ¼ S c [ S J and S c \ S J ¼ £ must hold. Boundaries of the different fields may overlap, though. 2.1. Dimensionless form A straightforward implementation of above governing equations resulted, in early implementations, in equation systems with condition numbers on the order of 1012, which hindered the convergence of the nonlinear iteration scheme. This can be attributed to the varying orders of magnitudes of the physical parameters of the chemical and the mechanical problem, given in Table 1. Therefore, the field equations have been normalized via the ansatz

c k2 xi ¼ k xi ; ui ¼ k ui ; c ¼  ; t ¼ t c t  ¼ t ; D X where k represents the characteristic length scale of the considered problem, xi dimensionless unit coordinates, ui the dimensionless displacement field, c the dimensionless concentration field, t  a ‘‘unit time’’, and tc the characteristic time scale of the problem. Furthermore, the material parameters Cijkl are scaled such that the stress field becomes dimensionless:



rij ¼ Cijkl ekl 

i  1h  c  cref dkl ; 3

ð8Þ 

 cref is the dimensionless reference concentration, C ¼ X Cijkl are the scaled components of the constitutive where cref ¼ X ijkl RT  X tensor, and rij ¼ RT rij are dimensionless stresses. This ansatz turns the governing equations into

229

P. Stein, B. Xu / Comput. Methods Appl. Mech. Engrg. 268 (2014) 225–244 Table 1 Material data for the governing equations’ coefficients, taken from [25]. Parameter

Value

 Partial molar volume X Gas constant R

8:3144621 J mol

Absolute temperature T Diffusion coefficient D

283 K

Young’s modulus E

1010 Pa 0:3

Poisson’s ratio

m

Electric current density in Faraday’s constant F

where

3:497  106 m3 mol 1

1

K1

7:08  1015 m2 s1

2 A m2 96:485336521 C mol

1

rij;j þ bi ¼ 0;

ð9aÞ

@ c  c;ii þ c;i rh;i þ c rh;ii ¼ 0; @ t

ð9bÞ

rh is the dimensionless hydrostatic stress, and bi are the scaled body forces, 

bi ¼

 kX bi : RT

Above equations are subject to the boundary conditions

^i u ^ i ; ¼u k  X rij nj ¼ ^ti ¼ ^ti ; RT  ^c ¼ ^c ; c ¼ X ui ¼

ð10aÞ ð10bÞ ð10cÞ

 c;i ni þ c rh;i ni ¼

 k in X D F

¼ ^f  :

ð10dÞ

As a consequence of this transformation, the equation systems’ condition number reduced to an order of 104 and the convergence behavior of the nonlinear equation solver improved significantly. Remark: In subsequent sections, only the dimensionless form of the governing equations will be employed. The superscript stars are hence omitted for notational simplicity. 2.2. Transformation into the weak form Following standard arguments as given in, e.g., [32], the weak form of the mechanical problem results from (9a) as

0¼

Z

deij rij dV þ

Z

B

dui t i dA þ

St

Z

dui bi dV:

ð11Þ

B

The weak form of the coupled diffusion equation is obtained in a similar fashion. Multiplication of (9b) with the virtual concentration field dc yields

Z dc c;i rh;i dV þ dc c rh;ii dV B B B B Z Z Z Z Z ¼ dc c_ dV þ dc;i c;i dV  dc c;i ni dA þ dc c rh;i ni dA  dc;i c rh;i dV:

0¼

Z

dc c_ dV 

B

Z

dc c;ii dV þ

B

Z

SJ

SJ

B

The two integrals over S J can be combined using (10d), giving the weak form of the concentration field problem:

Z Z Z 0¼ dc c_ dV þ dc;i c;i dV  dc;i c rh;i dV B B Z B  þ dc f dA:

ð12Þ

SJ

3. Discretization Eqs. (11) and (12) are discretized by the isoparametric ansatz

ui ¼ NI uIi ;

c ¼ N I cI ;

ð13Þ

230

P. Stein, B. Xu / Comput. Methods Appl. Mech. Engrg. 268 (2014) 225–244

with N I ¼ N I ðn; g; fÞ; uIi , and cI being the NURBS shape functions, the displacements, and the concentration associated with node/control point I, respectively. Details of the computation of NURBS functions and their properties are omitted here. Instead, we refer to the host of publications on NURBS and IGA, a small subset being given by [34,36,37,54,55]. In vector form–or Voigt notation–the constitutive law (8) reads



r¼C e

c  cref  / ; 3

ð14Þ

with

r ¼ ð rx ry rz rxy ryz rxz ÞT ; e ¼ ð ex ey ez 2exy 2eyz 2exz ÞT ; / ¼ ð 1 1 1 0 0 0 ÞT : The strains are computed from the nodal displacements via

e ¼ BuI uI ; using the strain–displacement matrix

2

NI;x

0

6 0 6 6 6 0 u BI ¼ 6 6N 6 I;y 6 4 0 NI;z

0

NI;y 0 NI;x NI;z 0

3

0 7 7 7 NI;z 7 7 0 7 7 7 NI;y 5 NI;x

and

 uI ¼ uI1

uI2

uI3

T

:

Similarly, the concentration gradient follows from

rc ¼ BcI cI ; where

 c T BI ¼ ½ NI;x

NI;y

NI;z :

By virtue of the shape functions’ higher continuity, the gradient of the hydrostatic stress can be expressed purely in terms of the displacements and the concentration field,

rrh ¼

 1 DI uI  f7 BcI cI ; 3

ð15Þ

where DI is a second-order operator matrix that is depending on the material coefficients Cijkl . Its computation is described, along with the definition of the scalar f7 , in Appendix A. 3.1. Second-order derivatives of NURBS functions Setup of the operator matrix DI requires the shape functions’ second derivatives with respect to the global coordinates. In order to obtain such expressions, we start from the common first-order derivatives with respect to global coordinates xi :

NI;i ¼ NI;nb nb;i ¼ NI;b nb;i ; where nb denotes the parametric variables describing a NURBS patch. Thereby, indices i; j, and k refer to global coordinates, whereas indices a; b, and c pertain to parametric coordinates. Another derivation with respect to global coordinates yields

NI;ij ¼ NI;bc nb;i nc;j þ NI;b nb;ij : The second-order derivatives N I;bc can be computed from the tensor product of univariate B-Spline functions’ derivatives, followed by the NURBS projection [56]. The first-order derivatives nb;i of the parametric coordinates with respect to global coordinates are evaluated as in default FEM shape function routines. Computation of the factors nb;ij , however, requires additional considerations. As stated by the isoparametric concept the physical coordinates xi result from the mapping

xi ¼ NI xIi : Using the relation

ð16Þ

P. Stein, B. Xu / Comput. Methods Appl. Mech. Engrg. 268 (2014) 225–244

231

@ xr ¼ drs ; @ xs the derivative of (16) with respect to the coordinate xj gives

dij ¼ NI;b nb;j xIi : Deriving this expression with respect to xk yields

0 ¼ NI;bc nb;j nc;k xIi þ NI;b xIi nb;jk : Due to (16),

NI;b xIi ¼ xi;b holds, leading to the equation system

xi;b nb;jk ¼ N I;bc nb;j nc;k xIi ;

ð17Þ

where the matrix xi;b is simply the Jacobian of the physical coordinates with respect to the parametric NURBS coordinates. These equations have to be solved at each Gauss point during the shape function evaluation. They have a unique solution as long as the isoparametric mapping (16) remains non-singular, which is the case even for highly distorted meshes [55]. However, above routine incurs high computational costs, as (17) has to be evaluated once for each second-order derivative nb;jk . Alternative procedures for their evaluation are described by Collier et al. [56]. 3.2. Consistent tangent stiffness Insertion of above ansatz into (11) and (12) produces the respective mechanical and chemical internal nodal forces for node I,

Z  u T PuI ¼  BI r dV; Z B Z  c T NI c_ dV þ BI rc dV PcI ¼ B Z B  c T  BI c rrh dV:

ð18aÞ

ð18bÞ

B

The external nodal forces arising from boundary tractions, flux, and body forces are considered in the global residual vector by means of the FE software. The nodal contributions to the consistent tangent matrix then result from

Z

Kuu IJ ¼ 

@PuI ¼ @uJ

Kuc IJ ¼ 

@PuI 1 ¼ 3 @cJ

Kcu IJ ¼ 

@PcI 1 ¼ @uJ 3

Z

Kcc IJ ¼ 

@PcI ¼ @cJ

Z



B

T BuI C BuJ dV;

B

B

Z



B

ð19aÞ

T BuI C / NJ dV;

ð19bÞ

 c T BI DJ c dV;  c T c BI BJ dV þ

ð19cÞ Z

 c T 1 BI NJ rrh dV  3 B

Z B

 c T c BI BJ c f 7 dV:

ð19dÞ

The residual expressions and the tangent stiffness matrices have been implemented in a beta version of the FEM software FEAP [57] that supports Isogeometric Analysis and NURBS shape functions. The nonlinear equations are solved using a Newton–Raphson iteration together with a backward Euler time integration scheme and the damping matrix

Dcc IJ ¼ 

@PcI 1 ¼ @ c_ J Dt

Z B

NcI NcJ dV:

ð20Þ

All other components of the damping matrix are zero. 4. Simulation of electrode particles during charging The particles that constitute a powder-structured, porous electrode are commonly described by spherical or ellipsoidal shapes. Although this provides—considering the particles’ irregular shape—only a coarse approximation, it allows to formulate analytical solutions for the governing equations described in Section 2. In this section, our approach is firstly verified using the analytical solution for a partially decoupled model. Secondly, the impact of the full coupling on the mechanical

232

P. Stein, B. Xu / Comput. Methods Appl. Mech. Engrg. 268 (2014) 225–244

behavior of a particle is demonstrated along with the influence of simulation parameters such as material stiffness, applied charge rate, and particle shape. 4.1. Spherical particle In order to describe the spherical particle, only an octant of the mesh has to be provided along with symmetry boundary conditions, namely the restriction of normal displacements and ion flux across the symmetry planes. The particle is described by a triquadratic NURBS mesh generated using the concept of operator- and template-based modeling [43]. Thereby, a NURBS line segment parallel to the x-axis is swept radially around the y-axis at the origin, resulting in the quadrant of a circular disc in the xz-plane. A sweep of said disc around the z-axis then yields the volumetric NURBS mesh. This coarse mesh of merely 27 control points is subsequently refined by knot refinement, whereby the tensor product structure of the patch, along with the underlying swept model structure, allows for independent refinement, e.g., in radial or circumferential direction. A rendering of the resulting mesh is illustrated in Fig. 2(a). Further modeling operations such as boundary extraction provide the NURBS surfaces for the definition of the mechanical and chemical boundary conditions, in particular for the integration of the scalar flux over the particle’s boundary. This can be easily cast into a template program taking as input the geometric parameters and the boundary values, and producing the desired mesh. Remark: The meshes created by this template are degenerated along the axis of revolution as well as in the origin. In these regions several control points coincide. However, we have not observed adverse effects on the simulations: as a consequence of their ‘‘variation-diminishing’’ property, NURBS meshes possess a higher robustness against distortions [55]. Furthermore, the integrals (18) and (19) are not evaluated directly at these singularities as the Gauss points are slightly offset from the element boundaries. The analytical solution given in [13] has been used as reference for verification. It describes the one-dimensional diffusion in a spherical particle, albeit not with the full coupling effect as given in (6). Instead, it is based on classical Fick diffusion in a sphere, where the concentration field affects the stress field according to (2), but experiences no feedback from the stresses. This model’s behavior can be readily reproduced by reducing the chemical residual (18b) to

PcI ¼

Z

NI c_ dV þ

B

Z B

 c T BI rc dV;

ð21Þ

and employing the terms

Kcu IJ ¼ 0;

ð22aÞ

Fig. 2. (a) NURBS mesh of a spherical electrode particle (b)–(d) Contour plots of the concentration field, the radial stresses, and the circumferential stresses predicted by the fully coupled model at t  0:25, i.e., when the stresses in the particle are at their respective peak levels.

P. Stein, B. Xu / Comput. Methods Appl. Mech. Engrg. 268 (2014) 225–244

Kcc IJ ¼ 

Z X

 c T c BI BJ dV

233

ð22bÞ

for the nodal tangent stiffness matrix. The mechanical residual (18a) and its corresponding tangent matrix terms (19a) and (19b) are not affected by this modification. Potentiostatic boundary conditions cðr ¼ R; tÞ ¼ cR have been prescribed on the particle’s exterior boundary. The particle is assumed to be embedded in a compliant matrix and to be unrestrained by the surrounding porous electrode [5]. Hence, stress-free boundary conditions are applied together with the material parameters given in Table 1. Given that the typical size of particles is in the range of microns or lower [5], the body forces are neglected. The simulation runs have been using a time step size of D t ¼ 5  104 , a particle size of k ¼ 106 m, and a zero reference concentration. Fig. 3 shows the concentration profiles as predicted by our partially decoupled model together with the dimensionless solution extracted from [13]. One can recognize a good correspondence of the numerical results with the analytical model. A slight deviation can be seen for t ¼ 103 which is the result of oscillations in the solution during the first time steps. The reduced model with its one-way coupling of the solution fields is attractive in that it allows problem treatment with C 0 -FEM. However, the stresses acting in a particle can have a notable effect on the diffusion process. Therefore, a constant charge current density in ¼ 2 A m2 has been applied on the boundary of the particle instead of the constant concentration cR , which, eventually leading to a homogeneous concentration in the particle, allows to observe only a weak coupling effect. The simulation has been performed for a slightly larger particle with a radius of k ¼ 5  106 m and with a time step size of D t ¼ 5  103 . All other parameters remain as before. In order to facilitate convergence during the very first time steps, the magnitude of the flux boundary condition has been scaled with a ramp function,

 ^f ¼ a Xk in ; D F

ð23Þ

where

a¼



t 0:1

for 0 6 t < 0:1;

1

for 0:1 6 t:

Fig. 4 compares the evolution of the concentration field predicted by the fully coupled and the simplified model. In the initial stages of charging, both models produce nearly identical concentration profiles. The diffusion wave reaches the center at t  0:1, around the time when the particle is charged with the full rate. Starting at t  0:25, the concentration curves detach slightly from each other. From this point on, the fully coupled model produces, compared with the decoupled model, lower concentrations at the particle’s exterior and higher concentrations at its center. Over time the gradients of the concentration field decay. The profile curves predicted by the simplified model at different time steps, however, are similar. Figs. 5 and 6 show the evolution of the radial stresses and circumferential stresses, respectively. They allow to observe the emergence of a tensile core and a compressive shell in the particle. Christensen and Newman [12] relate this to the lattice mismatch that is incurred by the concentration gradients between the lithium-rich shell and the relatively lithium-poor particle core. That is, the respective lattices try to accommodate their concentration-dependent lattice constants whereby the lithium-rich phase with its larger lattice constant pulls apart the lithium-poor phase. At the same time, the shell expands outward but, due to the spherical constraints, it is held back by the particle core, leading to large radial stresses in the core region [15].

1.0

Model Cheng (2008) Decoupled Time step t = 0.0010 t = 0.0100 t = 0.0574 t = 0.1000 t = 0.2000 t = 0.4000

c/cR [−]

0.8

0.6

0.4

0.2

0.0 0.0

0.2

0.4

0.6

0.8

1.0

r/R [−] Fig. 3. Time evolution of the concentration field over the radius of a spherical electrode particle. Solid lines show the solution extracted from [13], whereas dashed lines denote the corresponding results obtained from our modified model.

234

P. Stein, B. Xu / Comput. Methods Appl. Mech. Engrg. 268 (2014) 225–244

50.0

Model Coupled Decoupled Time step t = 0.05 t = 0.10 t = 0.25 t = 0.55 t = 1.00

3

Concentration [kmol/m ]

45.0 40.0 35.0 30.0 25.0 20.0 15.0 10.0 5.0 1.0 0.0

0.2

0.4

0.6

0.8

1.0

r/R [−] Fig. 4. Evolution of the concentration field in a spherical electrode particle during charging. Dashed lines denote the modified, partially decoupled model and solid lines show the solution considering the stress effect.

Radial stress σr [MPa]

50.0

Model Coupled Decoupled Time step t = 0.05 t = 0.10 t = 0.25 t = 0.55 t = 1.00

40.0

30.0

20.0

10.0

0.0 0.0

0.2

0.4

0.6

0.8

1.0

r/R [−] Fig. 5. Evolution of radial stresses in a spherical electrode particle under galvanostatic charging.

Circumferential stress σφ [MPa]

50.0

Model Coupled Decoupled Time step t = 0.05 t = 0.10 t = 0.25 t = 0.55 t = 1.00

40.0 30.0 20.0 10.0 0.0 −10.0 −20.0 −30.0 −40.0 −50.0 0.0

0.2

0.4

0.6

0.8

1.0

r/R [−] Fig. 6. Time evolution of circumferential stresses in a spherical electrode particle.

It is this stress distribution that causes the deviation in the concentration profiles in Fig. 4. According to (4), the ions are driven away not only from regions of high concentration but also away from regions of high compressive hydrostatic stress. The stress distribution emerging in the particle consequently increases the ion diffusion away from the particle’s compressed

P. Stein, B. Xu / Comput. Methods Appl. Mech. Engrg. 268 (2014) 225–244

235

von Mises stress σV [MPa]

shell towards its core, resulting in lower concentration gradients within the particle. The partially decoupled model, lacking the stresses’ effect on the diffusion, maintains the gradients in the concentration field. This introduces a characteristic relaxation behavior [12,33]. The radial and circumferential stress of the decoupled model reach their maxima at t  0:55. They maintain these values for the remainder of the simulation. The corresponding stresses in the fully coupled model, on the other hand, reach their maximum values at an earlier stage, around t  0:25. More important, though, is that they exhibit lower absolute values and that they slowly reduce over time. The stresses in the particle depend, as a consequence of the stress-free boundary conditions, solely on the gradients of the concentration field. The enhanced diffusion produces smooth, nearly uniform concentration fields, thus resulting in lower stress levels within the material. This effect becomes stronger with both increasing concentrations and stresses [12]. The simplified model, however, is unable to reduce the concentration gradients and therefore remains on the plateau of the highest occurring stresses. A uniform concentration field in the particle, on the other hand, effectively represents a uniform eigenstrain that–the particle’s expansion is not restricted–can be equilibrated by deformation of the particle, hence producing a stress-free particle. Main driving force for the relaxation are the stresses emerging in the particle. Fig. 7 illustrates this through the evolution of the maximum von Mises stresses during the charging process as a function of the electrode’s stiffness. With increasing Young’s modulus one can observe not only steeper growth of the stresses and higher peak stresses, but also a stronger post-peak relaxation behavior. The softer electrodes, in turn, show lower stress levels. Instead of characteristic stress peaks, they exhibit more uniform stress levels throughout charging, with relaxation being but weakly present. Based on this evidence it has to be concluded that the stresses arising in soft electrodes remain too low in order to develop an effect on the diffusion. Hence, the mechanical behavior effectively reduces to that of a partially decoupled model where diffusion is driven purely by the concentration gradient. A secondary factor contributing to the relaxation is the applied current density. This can be observed in Fig. 8 which depicts the time evolution of the maximum von Mises stresses in a particle for different charge rates. Here, similar to the influence of the material stiffness, both the peak in the stress curve and the relaxation behavior become more pronounced for higher charge rates, with the initial slope of the stress evolution depending on the current density. For low charge rates, however, one can recognize only a slow increase of the stresses, combined with an almost imperceptible relaxation. Caused by the finite diffusion speed in the material, ions accumulate in the particle’s boundary layers for higher charge rates [1]. According to previous reasoning, the strong concentration gradients effect high stress levels at the particle’s shell

160.0

Stiffness 2 GPa 5 GPa 10 GPa 20 GPa 50 GPa

140.0 120.0 100.0 80.0 60.0 40.0 20.0 0.0 0.0

0.2

0.4

0.6

0.8

1.0

t [−]

von Mises stress σV [MPa]

Fig. 7. Evolution of the maximum von Mises stress in a spherical particle for different material stiffnesses.

160.0

Current 2 −1 A/m −2 A/m2 2 −3 A/m −5 A/m2

140.0 120.0 100.0 80.0 60.0 40.0 20.0 0.0 0.0

0.2

0.4

0.6

0.8

1.0

t [−] Fig. 8. Evolution of the maximum von Mises stress in a spherical particle over time for different charge current densities.

236

P. Stein, B. Xu / Comput. Methods Appl. Mech. Engrg. 268 (2014) 225–244

that, reaching a certain threshold, gain a notable influence on the diffusion process. As a result, ions are driven towards the core of the electrode particle at a higher rate. This stronger diffusion reduces concentration gradients over the particle’s radius, resulting in lower stresses.

4.2. Ellipsoidal particle Using the template for the spherical mesh, the ellipsoidal particle can be readily generated by a non-uniform scaling of the control points’ coordinates. Thereby, the ellipsoid’s half-axes a; b; c are pairwise aligned with the global coordinate axes x; y; z and a fixed ratio b : c ¼ 1 is maintained. Straightforward application of ða; b; cÞ ¼ ða; 1; 1Þ for various values of a produces particles with variable volume. In order to maintain a volume of 4=3 pk3 for all ellipsoidal particles–corresponding to the volume of a spherical particle with radius pffiffiffi k–the half axes b; c are computed for a given a from b ¼ c ¼ 1= a. This transformation can be directly integrated into the template for the spherical particle. The code defining the boundary conditions remains unchanged since the model’s structure and topology is not affected by the affine transformation. The resulting NURBS mesh is displayed in Fig. 9 along with contour plots of the concentration and the stress fields in an ellipsoidal particle. A prolate spheroidal coordinate system ðw; q; #Þ is imposed in order to relate Cartesian stresses to the curved geometry, whereby its axis of revolution is aligned with the x-axis. Fig. 10 illustrates the coordinate surfaces corresponding to w; q, and #, respectively. The details of the coordinate system, that is, the bijection between Cartesian and spheroidal coordinates, is detailed in Appendix B along with the transformation relations for stress components. For the following simulations a particle with an aspect ratio of 1:5 : 1 : 1 has been used where the individual semi-axes have been scaled by the characteristic length scale k ¼ 5  106 m. The particles are subject to a current density in ¼ 2 A m2 that is controlled by the ramp function (23). All other parameters and boundary conditions follow those given in Section 4.1. Being an object of revolution, the spheroid exhibits a circular cross-section along its semi-minor axis, i.e., within the plane defined by x ¼ q ¼ 0. The two tangential stresses rq and r# along this axis are plotted in Fig. 11. Similar to a spherical particle, a compressive shell and a tensile core arises in the ellipsoid. Other than the two spherical circumferential stresses r/ and rH , the two spheroidal stress components show notable deviations as the stress rq along the semi-major axis exceeds r# approximately by a factor of 1:5. The latter, in turn, possesses a nearly identical distribution and magnitude as the

Fig. 9. (a) NURBS mesh of an ellipsoidal electrode particle with the half-axes 1:5 : 1 : 1 (b) – (d) Contour plots of the concentration field and the Cartesian stresses rx and rz , respectively, for t  0:3, i.e., when the stresses in the particle attain their respective maxima.

P. Stein, B. Xu / Comput. Methods Appl. Mech. Engrg. 268 (2014) 225–244

237

Fig. 10. The prolate spheroidal coordinate system used in the transformation from Cartesian to ellipsoidal stresses. Surfaces of constant w describe confocal ellipsoids, whereas surfaces of constant q form hyperboloids of revolution.

Stresses semi−minor axis [MPa]

60.0

Stress σρ σϑ Time step t=0.07 t=0.35 t=0.70 t=1.00

40.0 20.0 0.0 −20.0 −40.0 −60.0 0.0

0.2

0.4

0.6

0.8

1.0

r/R [−] Fig. 11. Evolution of the stresses

rq and r# along the semi-minor axis of an ellipsoidal particle.

circumferential stresses in a spherical particle (shown in Fig. 6). This suggests that the spheroidal shape retains the stress state of the spherical particle within the cross-section x ¼ q ¼ 0. In the vicinity of the equator the particle exhibits only a low curvature. Hence, the diffusion-induced in-plane stresses act nearly parallel to the x-axis, effectively pushing the particle tip outwards. This expansion is constrained by the bulk of the material around the particle’s core, causing large tensile stresses in longitudinal direction, as can be seen in Fig. 12. The curvature increases with growing distance to the equator, though. As a result, the in-plane stresses are more and more deflected from the longitudinal axis. Accordingly, the horizontal components of the diffusion-induced stresses reduce steadily, effecting the decrease in the stresses along the x-axis. The evolution of the ellipsoidal stresses rq and r# on the particle surface is depicted in Fig. 13. Thereby, the stresses are tracked along curves of constant ðw; #Þ. Starting at the equator, that is, at q ¼ 0, one can observe a gradual reduction of the stresses’ magnitude, accompanied by a convergence of their respective absolute values until, at the very tip, both stress components agree. They remain in the compressive regime. Material points in the equator region are constrained in a manner similar to the spherical particle. The closer a material point is to the tip, the less constrained it becomes. This asymmetry causes comparatively strong deformations along their longitudinal axis, with the magnitude of the deformation being highest at the ellipsoid’s tip. This expansion equilibrates, and hence reduces, the stress levels in this region, as can be witnessed in Fig. 13. The radial stress components rw are zero

238

P. Stein, B. Xu / Comput. Methods Appl. Mech. Engrg. 268 (2014) 225–244

σx along the x−axis [MPa]

60.0

Timestep t = 0.07 t = 0.35 t = 0.7 t = 1.0

50.0 40.0 30.0 20.0 10.0 0.0 0.0

0.2

0.4

0.6

0.8

1.0

r/R [−]

Stresses along outer surface [MPa]

Fig. 12. Evolution of

rx along the x-axis of the ellipsoidal particle.

−10.0

Stress σρ σϑ Time step t=0.07 t=0.35 t=0.70 t=1.00

−15.0 −20.0 −25.0 −30.0 −35.0 −40.0 −45.0 −50.0 −55.0 0.0

0.2

0.4

0.6

0.8

1.0

ρ [−] Fig. 13. Evolution of the stresses

rq and r# along the curved outer surface of an ellipsoidal particle.

50.0

Axis x−axis z−axis Timestep t = 0.07 t = 0.35 t = 0.70 t = 1.00

3 Concentration [kmol/m ]

45.0 40.0 35.0 30.0 25.0 20.0 15.0 10.0 5.0 0.0 0.0

0.2

0.4

0.6

0.8

1.0

r/R [−] Fig. 14. Comparison of the concentration field along two semi-axes in an ellipsoidal particle for different time steps.

on the surface by virtue of the stress-free boundary conditions. The stress state at the tip is hence determined by the two inplane stresses, which under-run their counterparts at the equator approximately by a factor of 2. Consequently, the stress term in the diffusion Eq. (9b) is less pronounced at the tip than along the semi-minor axis. This can be observed in

P. Stein, B. Xu / Comput. Methods Appl. Mech. Engrg. 268 (2014) 225–244

239

Fig. 14: in the first stages of charging, the concentration profiles along both half-axes are nearly aligned. With increasing time, the concentration gradients along the semi-major axis increase slowly, whereas the gradients along the short half-axis decay. This indicates a higher diffusion rate along the short semi-axis. Similar to the spherical model, the ellipsoidal particle develops a tensile core-compressive shell structure, with the highest stresses now acting in a ‘‘belt’’ around the particle’s equator. This is also the region in which the maximum von Mises stresses arise–for all ellipsoidal particles observed in this study, both with fixed and variable volume. Upon variation of the particle size, namely the size of the semi-major axis, differences emerge in the magnitude of the stresses. This can be seen in Figs. 15 and 16. Together, they compare the evolution of the maximum von Mises stress in dependence of the particle aspect ratio. Similar to the spherical case, one can observe a distinct stress relaxation over time. For particles with constant short half-axis the peak stress levels increase with larger aspect ratios. They stay significantly above those observed in a spherical particle and converge towards the stress state in a cylinder [25]. The particles with variable equatorial radius, in contrast, show a different behavior. Here, the peak stress levels increase until reaching a maximum for a  3=2 (corresponding to an aspect ratio of 1:837). For higher aspect ratios, the peak stress levels fall below those appearing in a spherical particle. In a similar manner do the stress peaks become more pronounced for larger aspect ratios, followed by an almost instantaneous relaxation. This behavior can be related to the concentration gradients that develop along the short semi-axis. Maintenance of a constant radius, as in Fig. 15, preserves the gradients as arising in a spherical particle, and hence the principal stress distribution. Due to the transition from a sphere towards an ellipsoid, the constraints’ character around the equator approximates those governing in an infinitely long cylindrical particle. Accordingly, the peak stresses converge to those corresponding to a plane strain state. The same holds for particles with fixed volume. Here, however, the short semi-axis is reduced with increasing particle length. Thereby, the diffusion paths from a particle’s shell to its core are shortened. Consequently, the concentration gradients in the particle are reduced, producing, in turn, lower stress levels.

von Mises stress σV [MPa]

60.0

Aspect 1.0 1.25 1.5 2.0 2.5 3.0 4.0 5.0 6.0

50.0 40.0 30.0 20.0 10.0 0.0 0.0

0.2

0.4

0.6

0.8

1.0

t [−] Fig. 15. Evolution of the maximum von Mises stress of an ellipsoidal particle during charging. The particle has the semi-axes a : 1 : 1 where the parameter a takes on the shown values, resulting in a variable particle volume.

von Mises stress σV [MPa]

60.0

Aspect 1.0 1.25 1.5 2.0 2.5 3.0 4.0 5.0 6.0

50.0 40.0 30.0 20.0 10.0 0.0 0.0

0.2

0.4

0.6

0.8

1.0

t [−] pffiffiffi pffiffiffi Fig. 16. Evolution of the maximum von Mises stress of an ellipsoidal particle during charging. The particle has the semi-axes a : 1= a : 1= a where the parameter a takes on the shown values, thereby maintaining a fixed particle volume.

240

P. Stein, B. Xu / Comput. Methods Appl. Mech. Engrg. 268 (2014) 225–244

Fig. 17. Von Mises stress distributions in oblate and prolate ellipsoidal particles for different parameters a. All particles have the volume of a spherical particle with radius 5  106 m. The models are shown at the respective time steps corresponding to the peak stress levels. All particles exhibit the semi-axes pffiffiffi pffiffiffi a : 1= a : 1= a, with a indicated in each sub-figure.

In addition to spherical and prolate, i.e., needle-like ellipsoidal particles, the simulation framework has been applied to the study of oblate, that is, lens-shaped electrode particles. A comparison of the stress distributions arising in particles of varying aspect ratio is given in Fig. 17. As described before, all prolate particles exhibit a region of large stresses around their equatorial plane. The particles’ tips, on the other hand, become nearly stress-free. As depicted in Fig. 16, the peak stress levels are reached at earlier time steps with increasing aspect ratio. The overall stress levels are consistently below those appearing in a spherical particle, suggesting their suitability for electrode applications. With the exception of the particle in Fig. 17(c), the oblate particles also show a ‘‘belt’’ of high stresses around their equator, which is now their semi-major axis. However, they do not only reach their maximum stress levels at later stages (again with the exception of Fig. 17(c)), but they also experience much higher stresses acting in the material. This is particularly striking as the flat particles possess only short diffusion paths along their semi-minor axis–which has lead to a reduction of overall stresses in the prolate ellipsoidal particles. As to be expected, the concentration gradients along the short half-axis are lower than along the semi-major axes, where they experience a steep increase towards the outermost layers of the particle. Due to the particle shape, the least constrained expansion of the particle can take place along the out-of-plane direction, i.e., along the short half-axis. All in-plane displacements of the material are, on the other hand, subject to comparatively high rotationally symmetric constraints. This allows only for a limited equilibration of diffusion-induced stresses, causing the extreme stress levels in the particle, and producing, accordingly, an increased diffusion rate towards the particle’s core. The particle in Fig. 17(c) represents an intermediate stage of stress distribution between the flat, oblate particles and the spherical particle. Here, the bulk of the material is concentrated in an ellipsoid with a semi-minor axis a ¼ 0:75k, whereas its semi-major axes b ¼ c  1:154k only slightly deviate from the radius of the sphere, r ¼ k. Consequently, the concentration gradients along the different axes show only minor differences. The low curvature around the particle’s pole causes large stresses acting along the horizontal plane in that region which, as a consequence of the particle’s compact shape, cannot be reduced through deformation. Compared with prolate ellipsoids, the higher stress levels in the oblate spheroids mean that the latter are more likely to experience material fatigue and mechanical degradation. Hence, it must be concluded that oblate spheroid particles are illsuited for electrode applications. 5. Summary and future work In this article we have described a model for coupled diffusion processes in Lithium-ion batteries. Its implementation is based on Isogeometric Analysis, motivated by the gradient of the stress field that occurs in the concentration’s evolution

P. Stein, B. Xu / Comput. Methods Appl. Mech. Engrg. 268 (2014) 225–244

241

equation. Using IGA, this gradient can be expressed purely in terms of the displacements and the concentration field, and the corresponding differential operator can be readily constructed from the constitutive law. Meshes for the spherical and ellipsoidal particles can be generated using the concept of operator- and template-based modeling [43], along with appropriate boundary conditions. The resulting template models allow to investigate and explain the influence of material stiffness, charge current density, and particle shape on the resulting stress distributions. This provides a–so far missing–explanation for the reduced stress levels in thin, needle-like ellipsoidal electrode particles, and hence their higher resilience against diffusion-induced stresses. A common simplification for the coupled diffusion processes is the negligence of the stress gradients in the ion flux. This permits a simple implementation even with C 0 -FEM, but it suppresses not only a significant factor for the magnitude of the stress levels, but also a relaxation effect, i.e., a steady decay of the stresses due to the enhanced diffusion. This stress effect is particularly relevant for material of high stiffness and for high (dis-) charge current densities. It can also be significant, though, for soft materials that show a stiffening effect with increasing ion concentration [15]. With our model we aim at gaining insight into the mechanical processes taking place in an electrode during cycled charge processes. More specifically, we are interested in the influence of particle morphology, anisotropy, surface mechanical effects, and damage on an electrode’s functionality and longevity. The results could be used to set up guidelines for the design and manufacturing of more robust and longer-lasting battery electrodes. Being in the first stages of this endeavor, we focused on the behavior of a single particle under varying shape. We made explicit the common simplifying assumptions used in the simulation of intercalation-induced stresses, the most severe of which is probably the restriction to small-strain kinematics that has been chosen in this article for ease of implementation. The extension of the hydrostatic stress gradient-operator to finite strain kinematics should not be overly problematic. The model is already able to regard anisotropic elastic effects. However, it neglected anisotropies in both the lattice expansion and the ion diffusion, which must be attributed to the scarcity of experimentally determined material parameters. The immediate next steps in our research are their investigation and the corresponding modification of our model. Experiments indicate a higher robustness of nano-structured battery electrodes against mechanical degradation [13,22,58]. Such size effects cannot be reproduced with our current model. In order to overcome this situation, two approaches can be pursued, the first being the incorporation of surface mechanical effects. This has been done by Cheng and Verbrugge [13] for an analytical model of a spherical electrode particle. An alternative measure for the description of size effects are strain-gradient models in which the elastic energy depends on both strains and their gradients. An isogeometric treatment of such models has been undertaken in [59] for the 2D case of a purely mechanical problem. Here, the higher-order continuous shape functions allowed for a straightforward implementation of the model equations. However, particular care has to be taken in the description of additional boundary conditions such as linear dependencies among control point displacements. We believe that our modeling framework will be of great support here, as it allows to reference model entities such as boundary curves and surfaces, and it keeps track of the topological structure of the model. Certain electrode materials, for instance graphite, show several staged phase transformations during the charge process [6]. Our model neglects these transformations. Instead, it assumes constant, concentration-independent elastic and chemical material parameters. The description of the phase transformations in Lithium-ion batteries can be achieved by means of phase-field models [17]. Here, the governing equations are derived from expressions of the free energy, given in terms of a so-called order parameter–for instance the current concentration in a point. For the treatment of the higher-order derivatives that show up in the equations, the concept of Isogeometric Analysis seems to be well suited [46]. This implementation provides a basis for more elaborate material formulations and particle morphologies. With minimal adaption it should also be applicable to the description of defect/dopant diffusion in semiconductors or the interdiffusion of hydrogen into storage tank walls, both of which are known to produce diffusion-induced stresses. For such applications, the high geometric flexibility of Isogeometric Analysis together with operator- and template-based modeling provides an excellent basis for parameter studies.

Acknowledgments We would like to express our gratitude towards Professor R. L. Taylor for providing us with a developer version of FEAP that supports Isogeometric Analysis.

Appendix A. Operator matrix for the hydrostatic stress gradient Section 3 introduced a matrix relating the nodal displacements to the hydrostatic stress gradient. In order to compute this matrix, (14) is evaluated row-wise. The three normal stresses are then summed up, giving the hydrostatic stress for smallstrain kinematics in terms of the strain components and the concentration. This expression can be expanded with respect to the displacement components’ partial derivatives. A subsequent derivation with respect to global coordinates then yields an expression for the hydrostatic stress gradient. This expression can be transformed into an operator matrix that takes the nodal displacements and concentrations into the sought quantity. To that end, the factors

242

P. Stein, B. Xu / Comput. Methods Appl. Mech. Engrg. 268 (2014) 225–244

f1 ¼ C1111 þ C2211 þ C3311 ; f2 ¼ C1122 þ C2222 þ C3322 ; f3 ¼ C1133 þ C2233 þ C3333 ; f4 ¼ C1123 þ C2223 þ C3323 ; f5 ¼ C1112 þ C2212 þ C3312 ;

ðA:1Þ

f6 ¼ C1113 þ C2213 þ C3313 ; 1 1 f7 ¼ ðf1 þ f2 þ f3 Þ ¼ Ciijj ; 3 3 are defined. They immediately follow from collecting the stiffness coefficients corresponding to the individual second-order partial derivatives of single displacement components as well as to the concentration. This allows to set up the asymmetric, second-order operator matrix

2

D11

6 DI ¼ 4 D21 D31

D12

D13

3

D22

7 D23 5

D32

D33

that represents the contribution of a single set of nodal displacements to the hydrostatic stress gradient. It has the components

D11 ¼ f1 NI;xx þ f4 NI;xy þ f6 NI;xz ; D12 ¼ f2 NI;xy þ f4 NI;xx þ f5 NI;xz ; D13 ¼ f3 NI;xz þ f5 NI;xy þ f6 NI;xx ; D21 ¼ f1 NI;xy þ f4 NI;yy þ f6 NI;yz ; D22 ¼ f2 NI;yy þ f4 NI;xy þ f5 NI;yz ; D23 ¼ f3 NI;yz þ f5 NI;yy þ f6 NI;xy ; D31 ¼ f1 NI;xz þ f4 NI;yz þ f6 NI;zz ; D32 ¼ f2 NI;yz þ f4 NI;xz þ f5 NI;zz ; D33 ¼ f3 NI;zz þ f5 NI;yz þ f6 NI;xz : Appendix B. Prolate spheroidal coordinates Section 4.2 required the transformation of stresses from the global Cartesian coordinate system into a curvilinear, ellipsoidal system. To that end, the following prolate spheroidal coordinates are employed:

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi y2 þ z2 þ ðx þ aÞ2 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 þ y2 þ z2 þ ðx  aÞ ; qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 q¼  y2 þ z2 þ ðx þ aÞ2 2a qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  y2 þ z2 þ ðx  aÞ2 ;

y # ¼ arctan ; z w P 1; q 2 ½1; þ1 ; # 2 ½0; 2pÞ; w¼

1  2a

ðB:1Þ

 is the distance of the ellipsoid’s foci to the origin along the focal axis, in this case the x-axis. It can be computed where a  and c of the semi-major and semi-minor axis, respectively, by from the sizes b

¼ a

pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2  c2 : b

This curvilinear coordinate system has the inverse relations

 w q; x¼a qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  2   w  1 ð1  q2 Þ cos #; y¼a qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  2   w  1 ð1  q2 Þ sin #; z¼a leading to the tangent vectors

ðB:2Þ

P. Stein, B. Xu / Comput. Methods Appl. Mech. Engrg. 268 (2014) 225–244

243

 q; x;w ¼ a

 w; x;# ¼ 0; x; q ¼ a qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  2   w cos # w  1 ð1  q2 Þ a ; y;w ¼  2  w 1 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi    q cos # w2  1 ð1  q2 Þ a y ;q ¼  ; ð1  q2 Þ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi    sin # w2  1 ð1  q2 Þ; y;# ¼ a qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi   w sin # w2  1 ð1  q2 Þ a ; z;w ¼  2  w 1 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi    q sin # w2  1 ð1  q2 Þ a z;q ¼  ; ð1  q2 Þ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi    cos # w2  1 ð1  q2 Þ: z;# ¼ a

ðB:3Þ

The associated metric coefficients are

2 g 11 ¼ a

w2  q2 w2  1

;

w2  q2 ; 1  q2     2 w2  1 1  q2 : ¼a

ðB:4Þ

2 g 22 ¼ a g 33

All other metric coefficients vanish due to the orthogonality of the basis vectors. Eventually, the transformation from Cartesian to spheroidal stress components can be achieved by

2

3

2

3

rw rwq rw# rx rxy rxz 6 7 7 T6 4 rwq rq rq# 5 ¼ T 4 rxy ry ryz 5T; rw# rq# r# rxz ryz rz where the rows of T consist of the normalized covariant base vectors of the ðw; q; #Þ-coordinate system:

TT ¼

h

x;w pffiffiffiffiffi g 11

x ;q pffiffiffiffi g2

x;# pffiffiffiffiffi g 33

i

:

ðB:5Þ

Remark: From (B.3) it can be seen that the stress transformation becomes singular both along the line w ¼ 1 and for the points q ¼ 1, corresponding to the line between the ellipsoid’s focal points. The tangent vectors’ normalization by the metric coefficients, as performed in (B.5), cancels out the line singularity but retains the singularity at the focal points, i.e., for ðw; q; #Þ ¼ ð1; 1; #Þ. For this reason, above stress transformation has only been used for determining the spheroidal stresses along the boundary of a given particle and along its semi-minor axis. Along the x-axis, only Cartesian stresses have been regarded for the analysis of the particle’s mechanical behavior. References [1] S. Golmon, K. Maute, S.-H. Lee, M.L. Dunn, Stress generation in silicon particles during lithium insertion, Appl. Phys. Lett. 97 (2010) 033111. [2] J. Vetter, P. Novák, M.R. Wagner, C. Veit, K.-C. Möller, J.O. Besenhard, M. Winter, M. Wohlfahrt-Mehrens, C. Vogler, A. Hammouche, Ageing mechanisms in Lithium-ion batteries, J. Power Sources 147 (2005) 269–281. [3] M. Park, X. Zhang, M. Chung, G.B. Less, A.M. Sastry, A review of conduction phenomena in Li-ion batteries, J. Power Sources 195 (2010) 7904–7929. [4] U. Kasavajjula, C. Wang, A.J. Appleby, Nano- and bulk-silicon-based insertion anodes for Lithium-ion secondary cells, J. Power Sources 163 (2007) 1003–1039. [5] K. Zhao, M. Pharr, J.J. Vlassak, Z. Suo, Fracture of electrodes in Lithium-ion batteries caused by fast charging, J. Appl. Phys. 108 (2010) 073517. [6] Y. Qi, H. Guo, L.G. Hector Jr., A. Timmons, Threefold increase in the Young’s modulus of graphite negative electrode during lithium intercalation, J. Electrochem. Soc. 157 (2010) A558–A566. [7] R.E. García, Y.-M. Chiang, W.C. Carter, P. Limthongkul, C.M. Bishop, Microstructural modeling and design of rechargeable Lithium-ion batteries, J. Electrochem. Soc. 152 (2005) A255–A263. [8] R.E. García, Y.-M. Chiang, Spatially resolved modeling of microstructurally complex battery architectures, J. Electrochem. Soc. 154 (2007) A856–A864. [9] R. Purkayastha, R.M. McMeeking, A linearized model for lithium ion batteries and maps for their performance and failure, J. Appl. Mech. T. ASME 79 (2012) 031021. [10] A. Bower, P.R. Guduru, V. Sethuraman, A finite strain model of stress, diffusion, plastic flow and electrochemical reactions in a Lithium-ion half-cell, J. Mech. Phys. Solids 59 (2011) 804–828. [11] P. Ramadass, B. Haran, R. White, B.N. Popov, Mathematical modeling of the capacity fade of Li-ion cells, J. Power Sources 123 (2003) 230–240. [12] J. Christensen, J. Newman, Stress generation and fracture in lithium insertion materials, J. Solid State Electrochem. 10 (2006) 293–319. [13] Y.-T. Cheng, M.W. Verbrugge, The influence of surface mechanics on diffusion induced stresses within spherical nanoparticles, J. Appl. Phys. 10 (2008) 083521. [14] Y.-T. Cheng, M.W. Verbrugge, Evolution of stress within a spherical insertion electrode particle under potentiostatic and galvanostatic operation, J. Power Sources 190 (2009) 453–460.

244

P. Stein, B. Xu / Comput. Methods Appl. Mech. Engrg. 268 (2014) 225–244

[15] R. Deshpande, Y. Qi, Y.-T. Cheng, Effects of concentration-dependent elastic modulus on diffusion-induced stresses for battery applications, J. Electrochem. Soc. 157 (2010) A967–A971. [16] R. Deshpande, Y.-T. Cheng, M.W. Verbrugge, A. Timmons, Diffusion induced stresses and strain energy in a phase-transforming spherical electrode particle, J. Electrochem. Soc. 158 (2011) A718–A724. [17] M. Huttin, M. Kamlah, Phase-field modeling of stress generation in electrode particles of Lithium ion batteries, Appl. Phys. Lett. 101 (2012) 133902. [18] K. Zhao, M. Pharr, J.J. Vlassak, Z. Suo, Inelastic hosts as electrodes for high-capacity Lithium-ion batteries, J. Appl. Phys. 109 (2011) 016110. [19] K. Zhao, M. Pharr, Q. Wan, W.L. Wang, E. Kaxiras, J.J. Vlassak, Z. Suo, Concurrent reaction and plasticity during initial lithiation of crystalline silicon in Lithium-ion batteries, J. Electrochem. Soc. 159 (2012) A238–A243. [20] K. Bhandakkar, H. Tanmay, Gao, Cohesive modeling of crack nucleation in a cylindrical electrode under axisymmetric diffusion induced stresses, Int. J. Solids Struct. 48 (2011) 2304–2309. [21] J.H. Seo, M. Chung, M. Park, S.W. Han, X. Zhang, A.M. Sastry, Generation of realistic particle structures and simulations of internal stress: a numerical/ AFM study of LiMN2O4 particles, J. Electrochem. Soc. 158 (2011) A434–A442. [22] C.K. Chan, H. Peng, G. Liu, K. McIllwrath, X.F. Zhang, R.A. Huggins, Y. Cui, High-performance lithium battery anodes using silicon nanowires, Nature Nanotechnol. 3 (2008) 31–35. [23] L.Y. Beaulieu, K.W. Eberman, R.L. Turner, L.J. Krause, J.R. Dahn, Colossal reversible volume changes in Lithium alloys, Electrochem. Solid-State Lett. 4 (2001) A137–A140. [24] Y.F. Gao, M. Zhou, Strong stress-enhanced diffusion in amorphous lithium alloy nanowire electrodes, J. Appl. Phys. 109 (2011) 014310. [25] X. Zhang, W. Shyy, A.M. Sastry, Numerical simulation of intercalation-induced stress in Li-ion battery electrode particles, J. Electrochem. Soc. 154 (2007) A910–A916. [26] C.M. DeLuca, K. Maute, M.L. Dunn, Effects of electrode particle morphology on stress generation in silicon during lithium insertion, J. Power Sources 196 (2011) 9672–9681. [27] X. Zhang, A.M. Sastry, W. Shyy, Intercalation-induced stress and heat generation within single Lithium-ion battery cathode particles, J. Electrochem. Soc. 155 (2008) A542–A552. [28] J. Park, W. Lu, A.M. Sastry, Numerical simulation of stress evolution in lithium manganese dioxide particles due to coupled phase transition and intercalation, J. Electrochem. Soc. 158 (2011) A201–A206. [29] S. Kalnaus, K. Rhodes, C. Daniel, A study of lithium ion intercalation induced fracture of silicon particles used as anode material in Li-ion battery, J. Power Sources 196 (2011) 8116–8124. [30] M. Zhu, J. Park, A.M. Sastry, Fracture analysis of the cathode in Li-ion batteries: a simulation study, J. Electrochem. Soc. 159 (2012) A492–A498. [31] W.H. Woodford, Y.-M. Chiang, W.C. Carter, Electromechanical shock of intercalation electrodes: a fracture mechanics analysis, J. Electrochem. Soc. 157 (2010) A1052–A1059. [32] K.-J. Bathe, Finite Element Procedures, Prentice Hall, 1996. [33] Y.F. Gao, M. Cho, M. Zhou, Stress relaxation through interdiffusion in amorphous lithium alloy electrodes, J. Mech. Phys. Solids 61 (2013) 579–596. [34] T.J.R. Hughes, J.A. Cottrell, Y. Bazilevs, Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement, Comput. Methods Appl. Mech. Engrg. 194 (2005) 4135–4195. [35] Y. Bazilevs, V. Calo, J. Cottrell, J. Evans, T. Hughes, S. Lipton, M. Scott, T. Sederberg, Isogeometric analysis using T-Splines, Comput. Methods Appl. Mech. Engrg. 199 (2010) 229–263. [36] J.A. Cottrell, A. Reali, Y. Bazilevs, T.J.R. Hughes, Isogeometric analysis of structural vibrations, Comput. Methods Appl. Mech. Engrg. 195 (2006) 5257– 5296. [37] J.A. Cottrell, T.J.R. Hughes, A. Reali, Studies of refinement and continuity in isogeometric structural analysis, Comput. Methods Appl. Mech. Engrg. 196 (2007) 4160–4183. [38] D. Benson, Y. Bazilevs, M.-C. Hsu, T. Hughes, Isogeometric shell analysis: the Reissner–Mindlin shell, Comput. Methods Appl. Mech. Engrg. 199 (2010) 276–289. [39] J. Kiendl, Y. Bazilevs, M.-C. Hsu, R. Wüchner, K.-U. Bletzinger, The bending strip method for isogeometric analysis of Kirchhoff–Love shell structures comprised of multiple patches, Comput. Methods Appl. Mech. Engrg. 199 (2010) 2403–2416. [40] Y. Bazilevs, C. Michler, V. Calo, T.J.R. Hughes, Isogeometric variational multiscale modeling of wall-bounded turbulent flows with weakly enforced boundary conditions on unstretched meshes, Comput. Methods Appl. Mech. Engrg. 199 (2010) 780–790. [41] Y. Bazilevs, M.-C. Hsu, I. Akkerman, S. Wright, K. Takizawa, B. Henicke, T. Spielman, T.E. Tezduyar, 3D simulation of wind turbine rotors at full scale. Part I: geometry modeling and aerodynamics, Int. J. Numer. Methods Fluids 65 (2011) 207–235. [42] Y. Bazilevs, M.-C. Hsu, J. Kiendl, R. Wüchner, K.-U. Bletzinger, 3D simulation of wind turbine rotors at full scale. Part II: fluid-structure interaction modeling with composite blades, Int J. Numer. Methods Fluids 65 (2011) 236–253. [43] P. Stein, M.-C. Hsu, Y. Bazilevs, K. Beucke, Operator- and template-based modeling of solid geometry for Isogeometric Analysis with application to Vertical Axis Wind Turbine simulation, Comput. Methods Appl. Mech. Engrg. 213 (2012) 71–83. [44] S. Cho, S.-H. Ha, Isogeometric shape design optimization: exact geometry and enhanced sensitivity, Struct. Multi. Optim. 38 (2008) 53–70. [45] W.A. Wall, M.A. Frenzel, C. Cyron, Isogeometric structural shape optimization, Comput. Methods Appl. Mech. Engrg. 197 (2008) 2976–2988. [46] H. Gómez, V.M. Calo, Y. Bazilevs, T.J.R. Hughes, Isogeometric analysis of the Cahn–Hilliard phase-field model, Comput. Methods Appl. Mech. Engrg. 49– 50 (2008) 4333–4352. [47] M.J. Borden, C.V. Verhoosel, M.A. Scott, T.J.R. Hughes, C.M. Landis, A phase-field description of dynamic brittle failure, Comput. Methods Appl. Mech. Engrg. 217–220 (2012) 77–95. [48] D. Anders, C. Hesch, K. Weinberg, Computational modeling of phase separation and coarsening in solder alloys, Int. J. Solids Struct. 49 (2012) 1557– 1572. [49] F. Yang, Interaction between diffusion and chemical stresses, Mater. Sci. Engrg. A 409 (2005) 153–159. [50] S. Prussin, Generation and distribution of dislocations by solute diffusion, J. Appl. Phys. 32 (1961) 1876–1881. [51] J.C.-M. Li, Physical chemistry of some microstructural phenomena, Metall. Mater. Trans. A 9 (1978) 1353–1380. [52] G.K. Singh, G. Ceder, M.Z. Bazant, Intercalation dynamics in rechargeable battery materials: general theory and phase-transformation waves in LiFePO4, Electrochim. Acta 53 (2008) 7599–7613. [53] M.E. Stournara, P.R. Guduru, V.B. Shenoy, Elastic behavior of crystalline LiSn phases with increasing Li concentration, J. Power Sources 208 (2012) 165– 169. [54] L. Piegl, W. Tiller, The NURBS Book, second ed., Springer, 1997. [55] S. Lipton, J. Evans, Y. Bazilevs, T. Elguedj, T. Hughes, Robustness of isogeometric structural discretizations under severe mesh distortion, Comput. Methods Appl. Mech. Engrg. 199 (2010) 357–373. [56] N. Collier, L. Dalcin, V.M. Calo, Petiga: high-performance Isogeometric Analysis, ArXiv e-prints arXiv:1305.4452 [cs.MS] (2013). [57] FEAP–a finite element analysis program, . [58] C.-Y. Sun, Y.-G. Zhu, T.-J. Zhu, J. Xie, G.-S. Cao, X.-B. Zhao, S.-C. Zhang, LiMn2O4 microspheres secondary structure of nanoparticles/plates as cathodes for Li-ion batteries, Mater. Res. 28 (2013) 1343–1348. [59] P. Fischer, M. Klassen, J. Mergheim, P. Steinmann, R. Müller, Isogeometric analysis of 2D gradient elasticity, Comput. Mech. 47 (2011) 325–334.