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\lx@sectionsign,Department of Mathematics & Statistics, McMaster University, Hamilton, Canada, L8S 4K1 (jamie.foster@port.ac.uk, bprotas@mcmaster.ca). \lx@sectionsignCurrent 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-10μm) 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 50% by volume (75% 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 (LiC6). 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 (LiCoO2), iron
phosphate (LiFePO4) and nickel manganese cobalt oxide
(LiNiMnCoO2)) 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 O(1−100)GPa [45, 13, 37, 15]) in comparison to the
mechanical moduli of the binder (typically O(1)MPa, see table LABEL:numbers). Moreover, we restrict our interest to electrode particles that can be considered isotropic, at least on lengthscales of O(1)μ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 2M+1×N square ‘unit’ tiles (whose sides have length Δ and
where M,N∈N), each of which contains a rigid circle
of active electrode material of radius r∗0 at its center. We
define the half-width and thickness of the whole electrode to be L=(M+1/2)Δ and h=NΔ respectively. Denoting the spatial
coordinates by (x∗1,x∗2), where we adopt the convention of
denoting dimensional quantities by a star, we define the following
external boundary segments, illustrated in Figure LABEL:sketch:
The location of the 2M+1×N circular internal boundary
segments, where the active electrode particles meet the porous binder
and electrolyte, can be parametrized as follows:
for m=−M,...,M and n=1,...,N. 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, τ∼10hrs) and body forces we arrive at
∂σ∗ij∂x∗i=∂p∗∂x∗j,
\hb@xt@.01(2.5)
where p∗ is the fluid pressure and σ∗ij are the
components of the (symmetric) stress tensor with the indices i
and j running over the values 1 and 2, 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
ϵ∗ij (again, for a REV of the porous composite)
are defined in terms of the components of the deformation
vector, u∗i, 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, s∗ij and e∗ij are the components of the
deviatoric (traceless) stress and strain tensors (respectively), S∗
and E∗ 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, G∗k and K∗k (for k=τ,1,2) 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
β∗abs 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 Γ∗1, 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, Γ∗2 and Γ∗4, 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,
Γ∗3, 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 Γ∗3. 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, Γ∗5, and take
βabs (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 xi. In summary, the boundary
conditions on the porous skeleton are
onΓ∗1:u∗1=0\hb@xt@.01(2.11)u∗2=0,\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 r∗0 to a radius
r∗0(1+g∗(t∗)), and allow for a displacement v∗i,m,n(t∗) of the centre of
the (m,n)-th cylinder. This displacement is determined by imposing that
there is no net force on the cylinder, so that^{1}^{1}1If we did not
allow for the displacement vi,m,n(t∗) 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 (n1,n2) is the normal to Γ∗m,n5.
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
att∗=0u∗i|t∗=0=0\hb@xt@.01(2.22)σ∗ij|t∗=0=0.\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 ϕ∗f, w∗i, k, μ, 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, ϕ∗f+ϕ∗s=1, where ϕ∗s and
ϕ∗f 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 ϕ∗f(w∗i−∂u∗i/∂t∗) from (LABEL:Peq),
in the standard fashion, gives
∂2u∗i∂t∗∂x∗i=kμ∂2p∗∂x∗2i.
The typical pressure scale P0 can be estimated from this equation in
terms of typical time, displacement and length scales τ, U0
and L, 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 L to be the
electrode half-width, i.e.L∼1cm. On referring to table
LABEL:numbers, we can compute a generous estimate of the pressure
P0∼101Pa. 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 P0 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 Σ0 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 h∼100μ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
Σ0∼105Pa which is a factor of 104 larger than P0
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:
Ratio of timescales for shear relaxation and cell (dis-)charge
O(10−2)
G2
Ratio of the second and first viscoelastic shear moduli
O(1)
Kτ
Ratio of timescale for bulk relaxation and cell (dis-)charge
O(10−2)
K1
Ratio of the first bulk modulus to the first shear modulus
O(1)
K2
Ratio of the second bulk modulus to the first shear modulus
O(1)
γ
Ratio of the electrode thickness and width
O(10−2)
λ
Ratio of the length of a ‘unit tile’ to electrode thickness
O(10−1−10−2)
r0
Ratio of the radius of an electrode particle to length of a ‘unit tile’
O(10−1−10−2)
Table 2.1: Descriptions and estimates of both the dimensional (upper portion) and dimensionless parameters (lower portion).
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 β0 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:
Estimates of these parameters depend upon the different lengthscales
in problem (Δ, r∗0, h and L), 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 LiPF6 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 k=O(1015)m2 is derived — which is
comparable to that of sandstone. We base our estimates of β0 on
the data provided in [11, 36, 41], where it is stated that the
densities of EC:DMC, PVDF and CMC are 1.2g/cc, 1.8g/cc and
1.6g/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 β0 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 E(t)=E∞+E0exp(−t/t0) with E∞=0.52MPa, E0=1.03MPa and t0=700s which gives rise to the estimates of G∗i and K∗i (for
i=1,2,τ) 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 Γ1:
u1=0
\hb@xt@.01(2.42)
u2=0,
\hb@xt@.01(2.43)
on Γ2:
σ11=0
\hb@xt@.01(2.44)
σ12=0,
\hb@xt@.01(2.45)
on Γ3:
σ22=0
\hb@xt@.01(2.46)
σ12=0,
\hb@xt@.01(2.47)
on Γ4:
σ11=0
\hb@xt@.01(2.48)
σ12=0,
\hb@xt@.01(2.49)
on Γ5:
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 t=0:
ui|t=0=0
\hb@xt@.01(2.52)
σij|t=0=0
\hb@xt@.01(2.53)
and the additional constraint
∫Γm,n5σijnjds=0.
\hb@xt@.01(2.54)
In dimensionless variables, the locations of the various boundary segments are
where x1,m=mλ, x2,n=(n−1/2)λ for m=−M,…,M and
n=1,…,N.
3 Upscaling and the macroscale problem
There are two predominant scales in this problem: there is a
microscale O(1−10μm) defined by the typical distance
between electrode particles and a macroscale defined by the dimensions
of the electrode^{2}^{2}2In fact it could be argued that the
macroscale is subdivided into a mesoscale of O(100μm),
defined by the electrode height, and a genuine macroscale of O(1cm), defined by the electrode width. However, for the
current purposes it proves convenient to incorporate both these
scales into a single macroscale.O(1mm). 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
Γ1, 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
Γk for k=1,…,4 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 βeffam — 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, βeffam(t)
is a function of time only. Thus the remaining macroscale equations
(to be solved on
Ωmac, see Figure LABEL:macro_scheme) are given by:
Geffτ∂seffij∂t+seffij=Geff2Geffτ∂eeffij∂t+eeffij,\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 Γ2 and
Γ4) for aspect ratio γ≪1. We rescale x1 with
1/γ and take the limit γ→0. On taking the
leading-order terms, integrating with respect to x2 and imposing
the boundary conditions (LABEL:barf1) and (LABEL:barf2) we find that
σeff12=σeff22=0. Thus
Substituting (LABEL:ozzyc) and (LABEL:gusc) into the 12 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,
σeff11 and ∂ueff2/∂x2, are derived by substituting
(LABEL:ozzy) and (LABEL:gus) into the 11 (or 22) component of
(LABEL:eq:Geff) and (LABEL:eq:Keff); they are
On using Laplace transforms to solve the above equations, integrating
the resulting expression for ∂u2/∂x2 with
respect to x2, and imposing the boundary condition (LABEL:kelly) we
find
ueff2
=
\hb@xt@.01(3.25)
σeff11
=
\hb@xt@.01(3.26)
where L−1 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
x2-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 Γ2,4. 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 γ→0+ and the two solutions are almost indistinguishable for 1/γ=10.
Figure 3.2: Panels (a)-(c) show numerical solutions for the deformation
field (u1,u2) to the macroscale equations
(LABEL:macro_start)-(LABEL:end_macro_bcs) for different values of
1/γ=1,5,10 respectively. Convergence to a one-dimensional
deformation is clearly observed as γ→0+. Panel (d) shows a comparison between
σ11|x1=0,x2=1/2(t) from the numerical simulations
and the asymptotic analysis, i.e. (LABEL:macro_sol_end).
4 The microscale problem
We now revisit the microscale problem, rescaling the governing
equations about a generic individual electrode particle at
(x1,x2)=(x1,m,x2,n) somewhere in the bulk of the electrode
by writing
We aim to solve (LABEL:nd_start)-(LABEL:smell_my_cheese) for the
microscale variables Xi, Ui on a single periodic unit cell
−1/2<X1<1/2, −1/2<X2<1/2. 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
σ12=0 on X1=±1/2 and X2=±1/2. In
addition, since there are no lateral macroscopic displacements, we
require U1=0 on X2=±1/2. 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
Since we have not derived the effective coefficients of the macroscale
problem, and therefore cannot accurately determine the macroscale
strain ϵeff22,
we do not impose ℓ(t) directly but calculate it self-consistently using the
condition that there is no normal macroscale stress in the vertical
direction (σeff22=0), which implies
∫1/2−1/2σ22dX1=0 on X2=±1/2.
\hb@xt@.01(4.2)
Finally we have the swelling condition that Ui=Xig(t) on the
particle Γm,n5.
The microscale problem we have derived is symmetric in both X1 and
X2, so that we may solve it more efficiently by considering
the
quarter-cell domain, depicted in Figure LABEL:micro_scheme, with
boundary segments
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 ℓ(t) (a
tractable problem) within a root-finding loop that determines the
value of ℓ(t) 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 ℓ(t), only
two ‘test’ problems need to be solved per time step in order to find
the value of ℓ(t) 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 σij in favor of ϵij, E and S
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 vi∈C∞0(Ω) and integrated by
parts (in the usual way) to obtain the variational form