Sharp-Interface Formation during Lithium Intercalation into Silicon [original]

Sharp-Interface Formation during Lithium Intercalation into

Silicon

Esteban Meca∗Andreas M¨unch†Barbara Wagner∗‡

February 25, 2016

Abstract

In this study we present a phase-field model that describes the process of in-

tercalation of Li ions into a layer of an amorphous solid such as a-Si. The gov-

erning equations couple a viscous Cahn-Hilliard-Reaction model with elasticity in

the framework of the Cahn-Larch´e system. We discuss the parameter settings and

flux conditions at the free boundary that lead to the formation of phase boundaries

having a sharp gradient in ion concentration between the initial state of the solid

layer and the intercalated region. We carry out a matched asymptotic analysis to

derive the corresponding sharp-interface model that also takes into account the dy-

namics of triple points where the sharp interface in the bulk of the layer intersects

the free boundary. We numerically compare the interface motion predicted by the

sharp-interface model with the long-time dynamics of the phase-field model.

Key words. Asymptotic Analysis, Phase-Field Model, Interface Dynamics, Numer-

ical Methods

AMS subject classification 74N20, 35Q74, 35B25, 74S10

1 Introduction

Silicon electrodes for lithium-ion batteries are currently the subject of very intense re-

search to make silicon a practical alternative to graphite. Silicon can store a large

amount of lithium, but the large stresses that these electrodes undergo tend to cause

their fracture and pulverization after a few cycles. Various investigations to circumvent

these problems have led to designs of nanostructured electrodes such as arrays of pillars

or nanowires.

However, in order to carry out systematic and knowledge-based optimizations a fun-

damental understanding of the mechanisms involved in the intercalation process of silicon

∗Weierstrass Institute, Mohrenstraße 39, 10117 Berlin, Germany

†Mathematical Institute, University of Oxford, Andrew Wiles Building, Woodstock Road, Oxford,

OX2 6GG

‡Technische Universit¨at Berlin, Institute of Mathematics, Straße des 17. Juni 136, 10623 Berlin,

Germany

Sharp-Interface Formation during Lithium Intercalation into Silicon 2

itself is needed. This has lead to a number of experimental and theoretical investiga-

tions, partly with, seemingly, contradictory results. In particular, the experiments of

Sethuraman et al. [28] on the stress evolution of a silicon electrode during intercalation

have had a great impact on the modeling of silicon electrodes, as they have been taken as

a proof of the surprising result that lithiated amorphous silicon behaves plastically, see

for example the analysis in Bower et al. [2]. Their analysis is based on the assumption

that the film of amorphous silicon is thin enough with respect to the substrate to which

it is attached, so that the stress can be assumed uniform in the whole amorphous silicon

film. However, while the thickness of the film can be controlled, the uniformity of the

stress is an unknown In fact, it has been shown recently in Huang et al. [14] and Ngo et

al.[24]that this assumption uniformity of stress can indeed obscure the interpretation of

the results.

In addition, as has been shown more recently in Levitas et al. [18], the theoretical

yield stress is never reached which has lead to some phenomenological modeling to

explain yielding [33], sometimes also discarding plasticity [18]. On the other hand, the

models used to take into account plasticity have problems in the determination of the

parameters, sometimes suggesting a power law with exponent as high as 50 for the

constitutive law [4], which can be taken as a hint that the model may not be complete.

Further analysis for different electrode geometries involving cylindrical and spherical

silicon particles as well as annular structures have been carried out to investigate stress

evolution and deformation of lithiated silicon, see [9, 6, 7].

On the other hand, regarding the stuctural properties of silicon, it is known that after

the first intercalation cycle, the original crystalline silicon becomes amorphous, see the

review by McDowell et al. [21] and references therein. Moreover, as it has been shown

in McDowell et al.[20] and Wang et al.[30], the first intercalation of crystalline as well as

amorphous silicon occurs through a two-phase mechanism. The question whether this

two-phase process also occurs in amorphous silicon in subsequent intercalation cycles

has been raised in [20, 19] and a recent study by Cubuk et al. [8] based on molecular

dynamics simultations relates the two-phase lithiation of amorphous silicon with a sharp

structural transition in amorphous LixSi for x≈2 between a phase in which the Li atoms

are embedded in a covalent silicon matrix to a phase in which Si atoms are packed in

small clusters with few covalent bonds between them, surrounded by a dense, amorphous

structure of Li atoms, suggesting a transition from a mechanically Si-like phase to a softer

Li-like phase. Moreover, experiments show that this two-phase lithiation is self-limiting

in nanowires [19], presumedbly due to the high stresses generated.

In Meca et al. [22] a phase-field model that couples a viscous Cahn-Hilliard-Reaction

model with elasticity has been derived to describe the process of intercalation of Li-ions

into a layer of an amorphous silicon and investigate conditions leading in the long-time

limit to the formation of sharp phase boundaries between the original silicon phase and

the lithiated phase. The analysis of the long time dynamics of the emerging sharp phase

boundaries is the subject of the present study.

After presenting the phase-field formulation in Section 2, we show in the remainder

of this paper how an intercalation process is nucleated, for example due to some defect

Sharp-Interface Formation during Lithium Intercalation into Silicon 3

on the free absorbing boundary of a layer (or nanowire), and upon reaching a critical

lithium concentration, a phase transition sets in, which in turn leads to a sharp phase

boundary that then moves into the bulk of the layer and along its free boundary.

In Section 3 we derive for this regime a sharp-interface model using matched asymp-

totic expansions, whereby the analysis is divided into three regions, the sharp interface

analysis in the bulk of the layer, the analysis close to the free boundary and the analysis

of the triple points, where the sharp interface in the bulk of the layer intersects the free

boundary. All three regions have to be matched and finally yield the sharp-interface

model. Finally, we investigate in Section 4 numerically the dynamics of the long-time

limit of the phase-field model and compare evolution of the emerging phase boundaries

with expression for the velocities found from matched asymptotic analysis.

2 Formulation of the phase-field model

The system we consider consists of a thin layer of amorphous silicon resting on an

undeformable substrate. Lithium enters in the layer as a consequence of the difference of

electrochemical potential with the electrolyte. As the lithium concentration is increased,

the layer experiences a phase transformation from a poorly lithiated phase (a-Si) to a

heavily lithiated one (LixSi with x>2). Lithium insertion causes a stress-free strain.

For a discussion of the role of phase transformation on stress for this process we have

formulated a mathematical model in Meca et al. [22] and we will briefly introduce the

complete phase-field model here.

To facilitate the discussion of the results we confine our description of elasticity to

linear elasticity and follow the standard approach to the coupling of phase transitions

with linear elasticity, see e.g. [10] for details. From Hooke’s law the elastic energy is

given by

W=1

2Cijkl ij −0

ijkl −0

kl(2.1)

where the fourth order tensor Cijkl is defined as the elasticity tensor or as the stiffness

tensor,ij denotes the strain tensor

ij =1

2(∂jui+∂iuj) (2.2)

defined by the displacement field u, as the difference between the actual position of a

material point and the position in the undeformed material x(reference configuration).

The stress-free strain tensor, that is, the strain due e.g. to composition changes in the

absence of stress, is assumed to grow linearly with the concentration:

0

ij =α(c−¯c)δij (2.3)

where αand ¯care constants, possibly depending on the phase. Note that this choice

implies an isotropic strain change with concentration, but this needs not to be the case

and more general relations are possible.

Sharp-Interface Formation during Lithium Intercalation into Silicon 4

Using the symmetries of Cijkl the choice of elastic energy implies for the stress tensor

σij =Cijkl kl −0

kl,(2.4)

and assuming that the timescale of the elastic relaxation is much faster than that of

diffusion or the phase transformation, elastic equilibrium implies that the divergence of

the stress tensor is zero:

∂jσij = 0.(2.5)

Equation (2.5) will have to be fulfilled separately at each phase in the layer and it can

be written explicitly in terms of the displacement field using

σij =E

1 + νM

ij +ν

1−2νM

kkδij,(2.6)

where Eis Young’s modulus,νis Poisson’s ratio and M

ij =ij −0

ij is the mechanical

strain, the difference between the strain tensor and the stress-free strain tensor.

For the present problem the elastic properties of both phases are assumed different.

We limit this difference to Young’s modulus, since according to Shenoy et al. [29] there

is no clear tendency in the variation of νand obtain the equation:

E∇c∇u+∇uT−E0

1 + ν

1−2να(c−¯c)∇c+E0

1−2ν(∇ · u)∇c

2∇2u+1

2(1 −2ν)∇(∇ · u)−α1 + ν

1−2ν∇c= 0

(2.7)

where E0is the derivative of Young’s modulus with respect to concentration c.

For the transport of concentration cwe use the viscous Cahn-Hilliard model

∂tc=∇ · (M(c)∇(µ+χ ε ∂tc)) ,(2.8)

where M(c) is the mobility function, which, in the present study, is taken to be a constant.

The last term is the viscous term, see [25] and χcorresponds to a parameter with

dimensions of viscosity. The chemical potential

µ=1

NΩ

δF

δc =−γε ∇2c+γ

2c(1 −c)(1 −2c) + ∂cW(ij, c),(2.9)

is the variational derivative of the free energy

F=NΩZΩ

2γε |∇c|2+γ

εf(c) + W(ij, c) (2.10)

where f(c) is the usual double-well potential (free energy per particle) and W(ij, c) is

the elastic energy as defined in Eq. (2.1) for an isotropic elasticity tensor. The constant γ

carries the dimensions of energy times length, and NΩis the (global) number of particles

per surface. The parameter εis related with the interface width and the interfacial

energy between the lithiated and the unlithiated phases.

We nondimensionalize eqs. (2.7) and (2.8) via

Sharp-Interface Formation during Lithium Intercalation into Silicon 5

µ=µ∗γH−1

0, x =x∗H0, y =y∗H0, t =t∗H3

0M−1γ−1, ε =ε∗H0,

denoting with ∗the nondimensionalized variables. The characteristic length scale H0

corresponds to the height of the layer and Mis the constant mobility. For the elastic

variables, we apply the scalings

Cijkl =C∗

ijkl

ESi

2(1 + ν), ui=u∗

iαH0, ij =∗

ijα, σij =σ∗

αESi

2(1 + ν).

After dropping the ∗, the nondimensionalized problem can be written as

∂tc=∇2(µ+εβ ∂tc),(2.11a)

µ=−ε∇2c+1

εf0(c) + δ ∂cW(ij, c),(2.11b)

∂jσij = 0,(2.11c)

σij = 2Gij −0

ij+2ν

1−2νGkk −0

kkδij,(2.11d)

where the constitutive laws for the nondimensional shear modulus G=E(c)/ESi and

stress-free strain 0

ij are specified as

G= 1 + g(c)ELixSi

ESi

−1, 0

ij =h(c),

and the derivative of the nondimensional elastic energy takes the form

∂cW(ij, c) = (1 −ν)G0

1−2ν∂1u2

1+∂2u2

2+1

2G0(∂1u2+∂2u1)2

+2νG0

1−2ν∂1u1∂2u2−2(1 + ν)

1−2ν(h(c)G)0∇ · u

+3(1 + ν)

1−2νh(c)2G0.(2.11e)

Here, h(c) and g(c) are interpolating functions such that g(0) = h(0) = 0 and g(1) =

h(1) = 1. Note that we have slightly generalized our previous choice (2.3) for 0

ij, while

letting ¯c= 0. Also, note that we have not defined a scaling for c, but this should define

0 and 1 to be the two equilibrium concentrations. For the boundaries in contact with

the substrate, we will take a no-flux/no-deformation boundary condition:

u= 0,n· ∇c= 0,n· ∇µ= 0,(2.11f)

where nis the normal vector to the surface. In the case of the boundaries in contact

with the electrolyte, we take a no-traction boundary condition and, following [5], assume

a consistent no-flux condition for c(also known as variational boundary condition), to-

gether with a Butler-Volmer type absorption condition for the chemical potential[32]

σ·n= 0,n· ∇c= 0,n· ∇µ=K(µ) = k1−eQ(µ−µext),(2.11g)

Loading more pages...