A mathematical model for mechanically-induced deterioration of the binder in lithium-ion electrodes

# A mathematical model for mechanically-induced deterioration of the binder in lithium-ion electrodes

J. M. Foster Department of Mathematics & Statistics, McMaster University, Hamilton, Canada, L8S 4K1 (jamie.foster@port.ac.uk, bprotas@mcmaster.ca). Current address: Department of Mathematics, University of Portsmouth, Portsmouth, UK, PO1 2UP.    S. J. Chapman Mathematical Institute, University of Oxford, Oxford, UK, OX2 6GG.    G. Richardson School of Mathematics, University of Southampton, Southampton, UK, SO17 1BJ.    B. Protas.

## Abstract

This study is concerned with modeling detrimental deformations of the binder phase within lithium-ion batteries that occur during cell assembly and usage. A two-dimensional poroviscoelastic model for the mechanical behavior of porous electrodes is formulated and posed on a geometry corresponding to a thin rectangular electrode, with a regular square array of microscopic circular electrode particles, stuck to a rigid base formed by the current collector. Deformation is forced both by (i) electrolyte absorption driven binder swelling, and; (ii) cyclic growth and shrinkage of electrode particles as the battery is charged and discharged. In order to deal with the complexity of the geometry the governing equations are upscaled to obtain macroscopic effective-medium equations. A solution to these equations is obtained, in the asymptotic limit that the height of the rectangular electrode is much smaller than its width, that shows the macroscopic deformation is one-dimensional, with growth confined to the vertical direction. The confinement of macroscopic deformations to one dimension is used to obtain boundary conditions on the microscopic problem for the deformations in a ’unit cell’ centered on a single electrode particle. The resulting microscale problem is solved using numerical (finite element) techniques. The two different forcing mechanisms are found to cause distinctly different patterns of deformation within the microstructure. Swelling of the binder induces stresses that tend to lead to binder delamination from the electrode particle surfaces in a direction parallel to the current collector, whilst cycling causes stresses that tend to lead to delamination orthogonal to that caused by swelling. The differences between the cycling-induced damage in both: (i) anodes and cathodes, and; (ii) fast and slow cycling are discussed. Finally, the model predictions are compared to microscopy images of nickel manganese cobalt oxide cathodes and a qualitative agreement is found.

## 1 Introduction

Much current research is focused on the development of lithium-ion batteries for use in a variety of areas, ranging from consumer electronics to the automotive industry, with the current market for such devices being in excess of $20bn (USD) per annum and set to increase further [29]. Many previous studies have considered how to model the electrochemical processes occurring within such cells, with the aim of informing optimal design. The majority of these studies have been, at least partially, based on the seminal models of Newman et al. [16, 17, 22, 23], which account for the complex battery microstructure using ‘averaged’ quantities—an approach that has been subsequently formalized in [48, 12]. Providing an exhaustive list, or even a summary, of the extant work on lithium-ion batteries is almost impossible and instead we point to the following reviews [38, 61, 3]. Of particular interest is modeling that can inform design improvements in terms of increased energy density, facilitation of higher charge/discharge rates, and improved safety and cycling lifetime. Here, we will focus on the last of these, specifically the (detrimental) morphological changes that can occur during cell fabrication and usage as a result of the development of internal mechanical stresses in response to (i) polymer binder swelling due to electrolyte absorption, and; (ii) dilation (and contraction) of the electrode particles as the cell is cycled. Lithium-ion batteries comprise four key constituents: (i) the anode (the negative electrode during discharge); (ii) the cathode (the positive electrode during discharge); (iii) the electrolyte, and; (iv) a porous separator that ensures unhindered passage of ionic current (via the electrolyte) while electronically isolating the two electrodes, see Figure LABEL:coin_cell.1. This assembly is sandwiched between two metallic current collectors (CCs). In commercial cells both the anode and cathode have a complex morphology; both are typically composed of small particles (1-10m) of active material (AM) embedded in a porous polymer binder matrix which (in the case of the cathode) is doped with highly conducting acetylene black nano-particles. The pores within the polymer matrix (of typical size 10–100nm) form channels through which the electrolyte can penetrate — see Figure LABEL:example_slices. The binder performs two key tasks; firstly it ensures the structural integrity of the electrode and secondly it serves to connect electrode particles electronically to the current collector. After manufacture, when such electrodes are assembled as part of a cell, they are soaked in electrolytic fluid which is absorbed by the binder and results in significant swelling of this polymer phase. Commonly used binders include polyvinylidene fluoride (PVDF) and carboxymethyl cellulose (CMC), and each of these materials is known to exhibit expansions of around by volume ( by weight) after begin immersed in EC:DMC for a few hours/days [11, 36, 41]. Even though this binder/acetylene black phase only constitutes a relatively small fraction of the material making up an electrode — typical figures for the AM:filler ratio by weight (here filler denotes both polymer binder and acetylene black) range from 70:30 to 90:10 in cathodes [32, 34] and 90:10 to 95:5 in anodes [5, 60] — a 50% volumetric expansion has the potential to lead to significant deformations and stresses (and ultimately morphological damage to the electrode) [50]. After assembly, when the cell is under operating conditions, the chemical potential difference between the positive and the negative electrodes leads to charge transfer events (redox/(de-)intercalation reactions) in the electrode particles. Importantly, these reactions are accompanied by volumetric changes in the electrode particles. By far the most commonly used active material in commercial anodes is graphitic carbon (LiC). On lithiation (associated with battery charging) these anode particles exhibit a significant volumetric expansion, of around 10% [14, 46]. In contrast typical cathode active materials (e.g. cobalt oxide (LiCoO), iron phosphate (LiFePO) and nickel manganese cobalt oxide (LiNiMnCoO)) exhibit smaller, but nonetheless appreciable, volumetric expansions in the range of 2–4% [59, 39]. This swelling/contraction of the active material phase can lead to large mechanical stresses within the electrodes which, in turn, can lead to delamination of the binder from the electrode particles and the CC. Notably delamination, as seen in Figure LABEL:example_slices occurs preferentially (at least early on) in the plane parallel to the current collector and it is one of the goals of this work to explain this pattern of delamination. In unfavorable cases, structural change within the electrode has detrimental effects on cell performance. For example, delamination of the binder from the surfaces of the electrode particles and CC (as can be seen in Figure LABEL:example_slices) has been proposed as a cause of degradation, and has been shown to cause both a measurable decrease in the ‘connectivity’ and the bulk (or ‘effective’) electronic conductivity of cathodes [19, 26, 27, 30, 33, 49, 52]. In addition to decreasing the efficacy of electron conduction, the delamination process can increase current densities (by reducing the area of electronic contact) and in turn enhance Joule heating, an effect which has the potential to cause to thermal runaway [55]. Delamination also exposes additional areas of the surfaces of the electrode particles directly to the electrolyte leading to the formation of extra solid electrolyte interface (SEI) layers [43] and to the consumption of active lithium (from the electrolyte), an effect that results in decreased battery capacity. Whilst a good deal of work has been focused on: (i) characterizing the mechanical processes occurring within various different types of active material, e.g. phase changes and particle fracturing [44, 9, 10, 63, 57]; (ii) discerning the impact of cycling-induced compression of, and the consequent reduction of pore space in, separator membranes on long-term cell performance [42, 24, 25], and; (iii) understanding the relatinship between the stack pressure, in e.g. commerical pouch cells, and the health of the cells [35, 6, 7]; much less has been done to elucidate the mechanical behavior of individual electrodes (including the binder) as a whole. Previous modeling in this area includes [40, 58]. The former considers the effects of electrode particle expansion and contraction on electrical connectivity when the surrounding matrix is composed entirely of carbon black particles. However, the absence of a binder material means these results are irrelevant to commercial cells. In contrast, whilst [58] does model the binder material, it focuses almost entirely on predicting stresses within the electrode particles. Here we focus primarily on modeling the deformations of the binder itself in response to: (i) its swelling in response to electrolyte absorption, and; (ii) swelling/contraction of the electrode particles as the cell is cycled. We work within an idealized 2-dimensional geometry as illustrated in Figure LABEL:sketch. The electrode is assumed to be initially rectangular, with its bottom surface bound to a rigid current collector, while the electrode particles are taken to be circular cylinders with centers located on a regular square lattice. The mechanical evolution of the electrode is described by a model based on Biot’s theory of poroelasticity [54], but with the additional feature that the porous skeleton is treated as a linear viscoelastic. Closely related models have been proposed previously for studying the behavior and interaction of soil and groundwater in mining applications, see e.g. [1, 4]. Our primary goal here is to use this model to predict the normal stresses, in the binder, on the surfaces of the electrode particles and infer from these whether (and how) delamination occurs. The remainder of this work is as follows. In §LABEL:prior_snap, we formulate (and non-dimensionalize) a model for the mechanical behavior of an electrode (estimating the sizes of the dimensionless parameters from existing physical data). In §LABEL:the_planet, we carry out an ad-hoc upscaling of the governing equations to obtain homogenized equations for the macroscopic deformations of the electrode and exploit the small aspect ratio of the electrode (via an asymptotic analysis) to find an approximate one-dimensional solution to this homogenized model which we verify against a full numerical solution. In §LABEL:time_passes, we investigate the microscale problem around a single electrode particle imposing boundary data compatible with the solution to the macroscopic model. The forcing for this microscale problem arises from: (i) increases in binder volume as it absorbs electrolyte, and; (ii) from electrode particle volume changes in response to cell cycling. These two forcing protocols lead to markedly different forms of deformation corresponding to distinct microstructural damage which can be observed in real imaging data such as that shown in Figure LABEL:example_slices. Finally, in section §LABEL:discandconc we draw our conclusions. ## 2 Problem formulation Any model capable of predicting the mechanical evolution of a composite electrode should, in principle, be multiphase — capturing the deformations of the polymer binder and the electrode particles, as well as the flow of the electrolyte through the pores within the binder. Here we assume that electrode particles can be treated as rigid bodies which change volume in response to lithium intercalation and de-intercalation during operation of the cell. The assumption of particle rigidity can be justified by the much larger mechanical moduli of the active material (typically GPa [45, 13, 37, 15]) in comparison to the mechanical moduli of the binder (typically MPa, see table LABEL:numbers). Moreover, we restrict our interest to electrode particles that can be considered isotropic, at least on lengthscales of m, so that when the particles expand/contract they do so without changing shape. In a typical porous electrode microstructure, the flow of the fluid (electrolyte) through the pore space exerts negligibly small stresses on the porous viscoelastic medium through which it flows (the binder). This means that, to a good approximation, the deformation of the binder matrix decouples from the flow of the electrolyte which allows us to determine the mechanical state of the binder without tracking fluid flow. Throughout we assume infinitesimal strain. ### 2.1 Geometry A sketch of the model geometry considered is shown in Figure LABEL:sketch. The electrode is comprised of an array of square ‘unit’ tiles (whose sides have length and where ), each of which contains a rigid circle of active electrode material of radius at its center. We define the half-width and thickness of the whole electrode to be and respectively. Denoting the spatial coordinates by , where we adopt the convention of denoting dimensional quantities by a star, we define the following external boundary segments, illustrated in Figure LABEL:sketch:  Γ∗1 ={(x∗1,x∗2)|−L≤x∗1≤L, x∗2=0},Γ∗2={(x∗1,x∗2)|x∗1=L, 0 The location of the circular internal boundary segments, where the active electrode particles meet the porous binder and electrolyte, can be parametrized as follows:  Γ∗m,n5={(x∗1,x∗2)|(x∗1−x∗1,m)2+(x∗2−x∗2,n)2=r∗20}wherex∗1,m=mΔ,x∗2,n=(n−1/2)Δ, \hb@xt@.01(2.3) for and . For convenience, we also introduce the following shorthand for the complete collection of all of these internal boundaries  Γ∗5=∪m=Mm=−M∪n=Nn=1Γ∗m,n5. \hb@xt@.01(2.4) ### 2.2 Governing equations for the porous binder A suitable governing equation for the deformation of the solid matrix, in (see Figure LABEL:sketch), is found by considering a balance of forces on an representative elementary volume (REV) of the porous medium (small by comparison to , but large by comparison to a typical pore diameter). On neglecting inertial effects (justified by the large time scale for electrochemical cycling, hrs) and body forces we arrive at  ∂σ∗ij∂x∗i=∂p∗∂x∗j, \hb@xt@.01(2.5) where is the fluid pressure and are the components of the (symmetric) stress tensor with the indices and running over the values and , and we use the Einstein summation convention for repeated indices. Since we work with infinitesimal strain theory, it is not important to distinguish between Eulerian and Lagrangian coordinates — in fact, with small strain theory, these two systems coincide. The components of the symmetric Cauchy (infinitesimal) strain tensor (again, for a REV of the porous composite) are defined in terms of the components of the deformation vector, , by  ϵ∗ij=12(∂u∗i∂x∗j+∂u∗j∂x∗i). \hb@xt@.01(2.6) In order to define constitutive equations it is helpful to decompose both stress and strain tensors into their volumetric and deviatoric parts as follows:  s∗ij=σ∗ij−δijS∗,whereS∗=12σ∗kk, \hb@xt@.01(2.7) e∗ij=ϵ∗ij−δijE∗,whereE∗=12ϵ∗kk. \hb@xt@.01(2.8) Here, and are the components of the deviatoric (traceless) stress and strain tensors (respectively), and are the volumetric parts of the stress and strain (respectively). Polymer binder materials have viscoelastic properties [53] and although there is evidence that some binders behave nonlinearly we choose here to describe their properties using the standard linear model of viscoelasticity (SLM) [18, 62]. The resulting constitutive equations are  G∗τ∂s∗ij∂t∗+s∗ij=G∗2G∗τ∂e∗ij∂t∗+G∗1e∗ij, \hb@xt@.01(2.9) K∗τ∂S∗∂t∗+S∗=K∗2K∗τ∂∂t∗(E∗−β∗abs)+K∗1(E∗−β∗abs), \hb@xt@.01(2.10) where, and (for ) are material constants associated with the shear and bulk deformations of the material respectively. In particular, those with the subscripts 1 and 2 are the viscoelastic moduli (with dimensions of pressure), while those with the subscript are the relaxation time scales. Notably, these material constants are drained coefficients — the term “drained” being used to highlight the fact that these quantities are measured by subjecting a REV of the porous medium to an external load and measuring the resulting strains (as functions of time) whilst allowing the fluid to freely drain out of or enter in to the medium at constant pressure. Finally, the function is the volumetric expansion of the binder due to absorption of the electrolyte which is taken to be a function of time only owing to the (roughly) uniform distributions of the polymer and electrolyte throughout the electrode. It should be also noted that we do not use convective time derivatives (upper convective, Helmholtz or other) in our formulation of the SLM (LABEL:end)-(LABEL:end2) since we are working with infinitesimal strains. ### 2.3 Boundary and initial conditions for the porous binder Zero displacement (no slip) conditions are imposed on the interface between the electrode and CC, i.e. on , see Figure LABEL:sketch. Although it is not obvious that the electrode should remain adhered to the CC, the images shown in Figure LABEL:example_slices indicate that there is little or no motion there — indentations in the CC, resulting from the pressure applied to the electrode during manufacture, remain close to their associated electrode particles even after cycling. On the lateral extremities of the electrode, and , zero stress conditions are imposed. The rationale for these conditions is simply that the adjacent material (the electrolyte) is passive, and therefore there is no means to provide any load/traction to these surfaces, see Figure LABEL:coin_cell.1. The boundary conditions that should be applied on the upper surface of the electrode, , are less obvious. Here, the electrode is in contact with the separator, which is in turn in contact with the counter electrode, its CC and finally (in the case of a coin cell) a spring which is used to ensure contact is maintained between the electrochemical components. The stiffness of this spring is very small compared to the moduli of the polymer binder, so that it is unable to cause any appreciable deformation, and furthermore the separator, being composed of a fibrous material or a flexible plastic film, is incapable of supporting significant shear stresses. We therefore apply stress-free conditions on the electrode’s upper surface . Finally, we assume that both: (i) the state of charge of the electrode, and by extension the volume (or the volumetric expansion) of the electrode particles, and; (ii) the degree of electrolyte absorption, are known a priori and we therefore apply Dirichlet conditions on the deformation on the electrode particle surfaces, , and take (the volumetric expansion of the binder due to absorption of electrolyte) to be a known function of time. We remark that even though non-zero deformations are permissible on some of the boundary segments, since we are working within infinitesimal deformation theory, the model does not constitute a free boundary problem and the boundary conditions are applied at fixed locations in the spatial coordinates . In summary, the boundary conditions on the porous skeleton are \hb@xt@.01(2.11) \hb@xt@.01(2.12)  onΓ∗2:  σ∗11=0 \hb@xt@.01(2.13)  σ∗12=0, \hb@xt@.01(2.14)  onΓ∗3:  σ∗22=0 \hb@xt@.01(2.15)  σ∗12=0, \hb@xt@.01(2.16)  onΓ∗4:  σ∗11=0 \hb@xt@.01(2.17)  σ∗12=0, \hb@xt@.01(2.18)  onΓ∗m,n5:  to0.0pt$u∗1=v∗1,m,n(t)+(x∗1−x∗1,m)g∗(t∗),$\hb@xt@.01(2.19)  onΓ∗m,n5:  to0.0pt$u∗2=v∗2,m,n(t)+(x∗2−x∗2,n)g∗(t∗).\$ \hb@xt@.01(2.20)

Here the final two conditions represent the uniform growth of the circular electrode particles from an initial radius to a radius , and allow for a displacement of the centre of the -th cylinder. This displacement is determined by imposing that there is no net force on the cylinder, so that111If we did not allow for the displacement and fixed the position of the cylinders, then they would apply a non-zero net force to the binder. This is equivalent to imagining the ends of the cylinders to be fastened to some support, rather than allowing them to move freely.

 ∫Γ∗m,n5σ∗ijnjds=0, \hb@xt@.01(2.21)

where is the normal to . In principle, we should also allow a (linearised) rotation of the cylinder, determined by the condition that there is no net torque. However, we will see later that the torque induced on the binder by ignoring such rotations is small.

To close the solid-state component of the model, it remains to specify initial data for the relevant quantities. We assume initially, prior to the addition of the electrolyte, that the electrode is in a zero stress state so that
\hb@xt@.01(2.22) \hb@xt@.01(2.23)

### 2.4 Governing equations for the electrolyte

It is usual to relate the flux of a fluid through a porous medium to the pressure gradient within the fluid via Darcy’s Law which, for a deformable porous medium takes the form [21]

 ϕ∗f(w∗i−∂u∗i∂t∗)=−kμ∂p∗∂x∗i, \hb@xt@.01(2.24)

where , , , , are the volume fraction of fluid, a component of the fluid velocity, the permeability of the porous medium and the fluid viscosity. The problem for the fluid component is closed by making further statements on the conservation of mass of each phase (the fluid and the solid skeleton). On noting that, by definition, , where and are the volume fractions of the solid and fluid phases respectively, these conservation equations are

 ∂∂t∗(ϕ∗f)+∂∂x∗i(w∗iϕ∗f)=0, \hb@xt@.01(2.25) ∂∂t∗(1−ϕ∗f)+∂∂x∗i(∂u∗i∂t∗(1−ϕ∗f))=0. \hb@xt@.01(2.26)

Summing the equations above and substituting for from (LABEL:Peq), in the standard fashion, gives

 ∂2u∗i∂t∗∂x∗i=kμ∂2p∗∂x∗2i.

The typical pressure scale can be estimated from this equation in terms of typical time, displacement and length scales , and , respectively. It is given by

 P0=μLU0kτ

and since the flow occurs along the width of the electrode (it cannot flow through the impermeable current collectors) we take to be the electrode half-width, i.e. cm. On referring to table LABEL:numbers, we can compute a generous estimate of the pressure Pa. In order to work out whether the fluid pressure plays a significant part in the viscoelastic deformation of the binder, we use our estimate for in the force balance equation (LABEL:start) in order to compare the sizes of the mechanical stress and fluid pressure gradient terms on the left- and right-hand side of this equation. The typical size for the stress tensor is obtained from the constitutive equations (LABEL:end)-(LABEL:end2) using an estimate for the size of the strain tensor; on using a lengthscale m (i.e. the electrode thickness), this gives

 Σ0=G∗1U0h.

On referring to table LABEL:numbers we obtain an estimate of the size of the stress tensor Pa which is a factor of larger than and we conclude that it is reasonable to neglect the fluid pressure term in (LABEL:start), by replacing it by the following purely mechanical force balance equation:

 ∂σ∗ij∂x∗i=0. \hb@xt@.01(2.27)

### 2.5 Non-dimensionalization

We non-dimensionalize the model by setting:

 x∗i =hxi, t∗ =τt, x∗i,m =hxi,m \hb@xt@.01(2.28) β∗abs =β0βabs, u∗i =β0hui, g∗ =β0g, \hb@xt@.01(2.29) ϵ∗ij =β0ϵij, e∗ij =β0eij, E∗ =β0E, \hb@xt@.01(2.30) σ∗ij =β0G∗1σij, s∗ij =β0G∗1sij, S∗ =β0G∗1S \hb@xt@.01(2.31) v∗i,m,n =β0hvi,m,n, \hb@xt@.01(2.32)

where is the typical size of volumetric expansion of the binder owing to electrolyte absorption. The non-dimensionalization leads to a system characterized by the dimensionless parameters:

 Gτ=G∗ττ,G2=G∗2G∗1,K∗τ=K∗ττ,K1=K∗1G∗1,K2=K∗2G∗1,γ=hL,λ=Δh,r0=r∗0h. \hb@xt@.01(2.33)

Estimates of these parameters depend upon the different lengthscales in problem (, , and ), as well as the timescale for cell (dis-)charge. We base these, in turn, on the devices studied in [33, 19] from which we obtain the estimates shown in table LABEL:numbers. We base our estimate of the electrolyte viscosity, , on a 1 molar LiPF solution in EC/DMC — one of the most common electrolytes used in both commercial and research cells — its value shown in table LABEL:numbers is taken from a Sigma-Aldrich data sheet (the viscosities of most other common battery electrolytes are similar). The electrode permeability can be estimated from, for example, the Carman-Kozeny formula [8, 31]. Using suitable values for the porosity and approximations of the electrode geometry a value of m is derived — which is comparable to that of sandstone. We base our estimates of on the data provided in [11, 36, 41], where it is stated that the densities of EC:DMC, PVDF and CMC are g/cc, g/cc and g/cc respectively. Furthermore, according to [11, 36, 41], a typical volumetric expansion for a binder, either PVDF or CMC, is around 75% by weight. Thus, an estimate for is in the range 0.4–0.6. Typical values for mechanical coefficients are more difficult to estimate. It appears that they depend strongly on not only the polymer under consideration, but also the carbon black content and processing conditions. The situation is further complicated by the fact that these polymers tend to soften as they absorb electrolyte — contrast, e.g. the values given in [58] against those in [53]. However, the relaxation modulus arising from experimental measurements of softened (via electrolyte absorption) PVDF has been given in [58] and a six term Prony series was proposed. Here we have opted to model the binder behavior using a SLM which corresponds to a two term Prony series [18]. We therefore approximate the time dependent modulus in [58] using with MPa, MPa and s which gives rise to the estimates of and (for ) given in table LABEL:numbers.

### 2.6 The dimensionless problem

Applying the scalings (LABEL:puppet1)-(LABEL:puppets_bum) to equations (LABEL:start)-(LABEL:hello_giles2) leads to the following dimensionless system for dimensionless (unstarred) variables

 ∂σij∂xj=0, \hb@xt@.01(2.34)
 ϵij=12(∂ui∂xj+∂uj∂x1), \hb@xt@.01(2.35)
 eij=ϵij−δijE, \hb@xt@.01(2.36)
 sij=σij−δijS, \hb@xt@.01(2.37)
 E=12ϵkk, \hb@xt@.01(2.38)
 S=12σkk, \hb@xt@.01(2.39)
 Gτ∂sij∂t+sij=G2Gτ∂eij∂t+eij, \hb@xt@.01(2.40)
 Kτ∂S∂t+S=K2Kτ∂∂t(E−βabs)+K1(E−βabs) \hb@xt@.01(2.41)

with boundary and initial conditions

on :
 u1=0 \hb@xt@.01(2.42)
 u2=0, \hb@xt@.01(2.43)
on :
 σ11=0 \hb@xt@.01(2.44)
 σ12=0, \hb@xt@.01(2.45)
on :
 σ22=0 \hb@xt@.01(2.46)
 σ12=0, \hb@xt@.01(2.47)
on :
 σ11=0 \hb@xt@.01(2.48)
 σ12=0, \hb@xt@.01(2.49)
on :
 u1=v1,m,n(t)+(x1−x1,m)g(t) \hb@xt@.01(2.50)
 u2=v2,m,n(t)+(x2−x2,n)g(t), \hb@xt@.01(2.51)
at :
 ui|t=0=0 \hb@xt@.01(2.52)
 σij|t=0=0 \hb@xt@.01(2.53)

 ∫Γm,n5σijnjds=0. \hb@xt@.01(2.54)

In dimensionless variables, the locations of the various boundary segments are

 Γ1 ={(x1,x2)|−1/γ

where , for and .

## 3 Upscaling and the macroscale problem

There are two predominant scales in this problem: there is a microscale defined by the typical distance between electrode particles and a macroscale defined by the dimensions of the electrode222In fact it could be argued that the macroscale is subdivided into a mesoscale of , defined by the electrode height, and a genuine macroscale of , defined by the electrode width. However, for the current purposes it proves convenient to incorporate both these scales into a single macroscale. . On the microscale we consider the local deformation of the binder matrix around individual particles as a result of either: (i) the binder swelling in response to soaking-up electrolyte, or, (ii) the electrode particles expanding and contracting in response to lithiation and delithiation; whereas on the macroscale we consider the bulk response of the entire material, including binder and electrode particles. A proper treatment of the macroscale problem requires homogenization of the microscale equations to obtain effective medium equations approximating the macroscopic behavior of the electrode (see for example [47, 48]). Here however we are primarily interested in the results of the microscale problem because we wish to find the microscopic stress distribution around an electrode particle in order to see whether, and if so how, delamination of the binder from the electrode particle occurs. In order to accomplish this we consider a locally periodic array of electrode particles and solve the problem in a single periodic ‘unit’ cell around one of these particles. Nevertheless, the imposition of appropriate boundary conditions on the edge of this cell still requires knowledge of the macroscopic solution. In particular, we seek to demonstrate that, because of its large aspect ratio, the macroscopic deformation of the composite material is one-dimensional, occurring in the direction normal to the current collector , except close to the edges of the electrode. Since the information we require from the macroscale is limited, we look to accomplish this without resorting to a lengthy homogenization procedure. Instead we use physical intuition to write down a plausible set of macroscale equations, and leave the task of systematically determining the macroscopic parameters in terms of the microscale geometry and parameters to a further work.

The simplified homogeneous geometry for the macroscale problem is shown in Figure LABEL:macro_scheme. We note that the locations of the external boundaries of the electrode are unchanged by the upscaling procedure. Thus, we retain the notation for the external boundary segments from the full problem. Explicit definitions of for are given in (LABEL:diddly)-(LABEL:dodo).

It is clear that the effective stress and strain at the macroscale, denoted by the superscript “eff”, will still satisfy (LABEL:nd_start)-(LABEL:strain), and we can define effective deviatoric and volumetric stresses as usual, so that

 ∂σeffij∂xj=0, \hb@xt@.01(3.1)
 ϵeffij=12⎛⎝∂ueffi∂xj+∂ueffj∂xi⎞⎠, \hb@xt@.01(3.2)
 seffij=σeffij−δijSeff, \hb@xt@.01(3.3)
 Seff=12seffkk, \hb@xt@.01(3.4)
 eeffij=ϵeffij−δijEeff, \hb@xt@.01(3.5)
 Eeff=12ϵeffkk. \hb@xt@.01(3.6)

What is less clear is the constitutive equation which should hold for the homogenised material. There is no particular reason it should correspond to a SLM (or a two term Prony series) — more likely there will be a continuum of timescales involved. Nevertheless, for the small amount of information we require, we expect that approximating the macroscopic constitutive behaviour by a SLM with effective coefficients is sufficient.

The final part of the homogenised model concerns the macroscopic effect of the changing volumes of the electrode particles on the microscopic length scale. This is captured by an “effective” growth function — given by the ratio of the volume expansion of an electrode particle to the volume the periodic cell in which it is embedded — which appears in the equations in the same way that the binder swelling does. Notably, because we assume that particles are uniformly distributed in space, is a function of time only. Thus the remaining macroscale equations (to be solved on , see Figure LABEL:macro_scheme) are given by:
\hb@xt@.01(3.7)

 \hb@xt@.01(3.8)

The boundary conditions on the outer edges of the electrode remain unchanged, and the conditions on the interface between the electrode and the embedded electrode particles are no longer needed. The boundary conditions that close the macroscopic problem are therefore

 onΓ1:
 ueff1=0 \hb@xt@.01(3.9)
 ueff2=0, \hb@xt@.01(3.10)
 onΓ2:
 σeff11=0 \hb@xt@.01(3.11)
 σeff12=0, \hb@xt@.01(3.12)
 onΓ3:
 σeff22=0 \hb@xt@.01(3.13)
 σeff12=0, \hb@xt@.01(3.14)
 onΓ4:
 σeff11=0 \hb@xt@.01(3.15)
 σeff12=0, \hb@xt@.01(3.16)
 att=0:
 ueffi|t=0=0 \hb@xt@.01(3.17)
 σeffij|t=0=0. \hb@xt@.01(3.18)

We seek an asymptotic solution to (LABEL:macro_start)-(LABEL:end_macro_bcs), valid throughout the bulk of the electrode (away from the lateral boundaries and ) for aspect ratio . We rescale with and take the limit . On taking the leading-order terms, integrating with respect to and imposing the boundary conditions (LABEL:barf1) and (LABEL:barf2) we find that . Thus

 σeff=(σeff11000),Seff=σeff112,s=(σeff11/200−σeff11/2). \hb@xt@.01(3.19)

 ϵeff=⎛⎜ ⎜ ⎜ ⎜⎝012∂ueff1∂x212∂ueff1∂x2∂ueff2∂x2⎞⎟ ⎟ ⎟ ⎟⎠,Eeff=12∂ueff2∂x2,eeff=⎛⎜ ⎜ ⎜ ⎜⎝−12∂ueff2∂x212∂ueff1∂x212∂ueff1∂x212∂ueff2∂x2⎞⎟ ⎟ ⎟ ⎟⎠. \hb@xt@.01(3.20)

Substituting (LABEL:ozzyc) and (LABEL:gusc) into the component of (LABEL:eq:Geff) gives

 Geff2Geffτ∂∂t(∂ueff1∂x2)+∂ueff1∂x2=0. \hb@xt@.01(3.21)

The only solution to this problem satisfying both initial data and the boundary condition (LABEL:bob) is

 ueff1=0. \hb@xt@.01(3.22)

A pair of coupled equations for the remaining unknowns, and , are derived by substituting (LABEL:ozzy) and (LABEL:gus) into the (or ) component of (LABEL:eq:Geff) and (LABEL:eq:Keff); they are

 Geffτ∂σeff11∂t+σeff11+Geff2Geffτ∂∂t(∂ueff2∂x2)+∂ueff2∂x2=0, \hb@xt@.01(3.23) Keffτ∂σeff11∂t+σeff11=Keff2Keffτ∂∂t(∂ueff2∂x2−(βeffabs+βeffam))+Keff1(∂ueff2∂x2−(βeffabs+βeffam)). \hb@xt@.01(3.24)

On using Laplace transforms to solve the above equations, integrating the resulting expression for with respect to , and imposing the boundary condition (LABEL:kelly) we find

 ueff2 = \hb@xt@.01(3.25) σeff11 = \hb@xt@.01(3.26)

where denotes the inverse Laplace transform.

Most importantly, we have found that the thin geometry of the electrode and the zero displacement condition on the current collector forces the deformation to be essentially one-dimensional (in the -direction) — see (LABEL:deaf). Crucially, as we will see in the subsequent section, the relatively simple structure of the solution on the macroscopic scale will allow us to construct appropriate boundary conditions to apply at the microscopic scale.

We note that the one-dimensional solution we have obtained only applies throughout the bu lk of the electrode, and does not satisfy the stress-free conditions imposed on . This apparent discrepancy can be resolved by rescaling spatial variables appropriately and studying the solution in boundary layers local to the lateral edges of the electrode. In the interests of brevity we do not present the details of this analysis here. Instead, to provide evidence that the one-dimensional solution is indeed correct in the bulk, we solve the problem (LABEL:macro_start)-(LABEL:end_macro_bcs) numerically and verify agreement with the analytical solution (details of the numerical approach are presented in §LABEL:sec:numer). A selection of numerical solutions with increasing aspect ratio (decreasing ) is shown in Figure LABEL:macro_def_examples — as one can see, the agreement throughout the bulk improves as and the two solutions are almost indistinguishable for .

## 4 The microscale problem

We now revisit the microscale problem, rescaling the governing equations about a generic individual electrode particle at somewhere in the bulk of the electrode by writing

 x2−x2,n=λX2,x1−x1,m=λX1,u1=v1,m,n+λU1,u2=v2,m,n+λU2. \hb@xt@.01(4.1)

We aim to solve (LABEL:nd_start)-(LABEL:smell_my_cheese) for the microscale variables , on a single periodic unit cell , . Boundary conditions on the microscale cell problem need to be chosen in order that the micro-solution is compatible with the one-dimensional solution to the macroscale problem (LABEL:deaf). Since there is no macroscopic shear, we require on and . In addition, since there are no lateral macroscopic displacements, we require on . On the other hand, averaging the 22 component of the microscopic strain over the unit cell, which yields the net relative displacement of the horizontal boundaries of the cell, should be also equal to the macroscopic strain, so that

 ∫1/2−1/2∂U2∂X2dX2=U2∣∣X2=1/2−U2∣∣X2=−1/2=ϵeff22=∂ueff2∂x2=2ℓ(t).

Since we have not derived the effective coefficients of the macroscale problem, and therefore cannot accurately determine the macroscale strain , we do not impose directly but calculate it self-consistently using the condition that there is no normal macroscale stress in the vertical direction (), which implies

 ∫1/2−1/2σ22dX1=0 on X2=±1/2. \hb@xt@.01(4.2)

Finally we have the swelling condition that on the particle .

The microscale problem we have derived is symmetric in both and , so that we may solve it more efficiently by considering the quarter-cell domain, depicted in Figure LABEL:micro_scheme, with boundary segments

 Γmic1={X1,X2|r00,X2>0,X21+X22=r20}. \hb@xt@.01(4.5)

The boundary and symmetry conditions are then

 onΓmic1:
 U2=0 \hb@xt@.01(4.6)
 σ12=0, \hb@xt@.01(4.7)
 onΓmic2:
 U1=0 \hb@xt@.01(4.8)
 σ12=0, \hb@xt@.01(4.9)
 onΓmic3:
 U2=ℓ(t) \hb@xt@.01(4.10)
 σ12=0, \hb@xt@.01(4.11)
 onΓmic4:
 U1=0 \hb@xt@.01(4.12)
 σ12=0, \hb@xt@.01(4.13)
 onΓmic5:
 U1=X1g(t) \hb@xt@.01(4.14)
 U2=X2g(t), \hb@xt@.01(4.15)

where is determined by the requirement that

 ∫Γmic3σ22dX1=0. \hb@xt@.01(4.16)

Finally, we note that the condition (LABEL:smell_my_cheese) is automatically satisfied due to symmetry and also that there is no net torque on the particle, so that our earlier decision not to allow for rotations is justified.

Owing to the geometry, solution of the microscale problem must be obtained using numerical techniques. We employ the finite element method, and implement this approach using the open source software FreeFEM++ [28]. The problem has two non-standard features, namely: (I) the constitutive equations involve time derivatives of the stress- and strain-fields, and; (II) the problem is subject to an integral constraint, namely (LABEL:constraint). The first of these difficulties is tackled by discretizing the governing equations in time using an explicit approximation; this leads to an elliptic boundary-value problem for the state variables at each time step that are of similar form to those for an elastic medium — information from the previous time step is retained in the form of source terms. The second difficulty is circumvented by embedding a code that solves the problem for a specified value of (a tractable problem) within a root-finding loop that determines the value of at the new time step that leads to a stress field satisfying (LABEL:constraint). Since the variation of the average normal load on the top surface is a linear function of , only two ‘test’ problems need to be solved per time step in order to find the value of satisfying (LABEL:constraint). We now detail the method of numerical solution to this microscale problem.

### 4.1 Numerical approach

Before implementation in FreeFem++, the system of equations must be written in a suitable variational form. Given the non-standard nature of our problem, we briefly outline the procedure for deriving this variational formulation. First, we take the partial derivative of the force balance equation (LABEL:nd_start) with respect to time and then eliminate in favor of , and using equations (LABEL:ee3)-(LABEL:nd_eq_end). On doing so, and denoting a partial derivative in time with a dot, we arrive at

 \hb@xt@.01(4.17)

This equation is then discretized explicitly in time using a forward Euler approximation, multiplied by a test function and integrated by parts (in the usual way) to obtain the variational form

 ∫Γa1ϵij(u(n+1))vjni+a2ϵij(u(n))vjni+