The effect of mechanical stress on lithium distribution and geometry optimisation for multi-material lithium-ion anodes

TL;DR

This paper presents a coupled mechanical-chemical model for lithium-ion anodes with multi-material structures, revealing how stress influences lithium distribution and voltage, and enabling optimization of anode design.

Contribution

It introduces a stress-dependent chemical potential model for multi-material anodes and demonstrates optimization of silicon core volume considering stress effects.

Findings

Stress significantly affects lithium distribution and OCV in multi-material anodes.

The model enables optimization of silicon core volume to improve anode performance.

Stress-coupling impacts lithium intercalation and maximum stress in the anode.

Abstract

A model is presented for predicting the open-circuit voltage (OCV) and lithium distribution within lithium-ion anodes containing multiple materials, coupling linear elasticity with a stress-dependent chemical potential. The model is applied to a spherical radially-symmetric nano-particle with a silicon core and a graphite shell, highlighting the large effect on lithium distribution and OCV caused by the stress-coupling. Various performance measures based on the expanded volume, the amount of lithium intercalated and the maximum stress induced, are calculated for a silicon core with a graphite shell to enable optimisation of the volume of the silicon core.

Equations

$$\textbf{E}^{e}=\frac{1}{2}\Big{[}\big{(}\boldsymbol{\nabla}^{*}\boldsymbol{u}^{*}\big{)}^{\text{T}}+\boldsymbol{\nabla}^{*}\boldsymbol{u}^{*}\Big{]}-\eta_{a}V^{m}_{a}c^{*}_{a}\boldsymbol{1}\quad\text{ in }\Omega_{a},$$

$${\mathbf{\sigma }}^{\ast }={\mathbb{C}}_a^{\ast }(c_a^{\ast }):{\text{E}}^e\text{ in }{\Omega }_a,$$

$${\mathbb{C}}_a^{\ast }(c_a^{\ast })=c_{ijkl}^{i,\ast }(c_a^{\ast })={\lambda }_a^{\ast }(c_a^{\ast }){\delta }_{ij}{\delta }_{kl}+G_a^{\ast }(c_a^{\ast }){\delta }_{ik}{\delta }_{jl}+G_a^{\ast }(c_a^{\ast }){\delta }_{il}{\delta }_{jk},$$

$$\boldsymbol{\sigma}^{*}=\mathbb{C}^{*}_{a}(c_{a}^{*}):\big{(}\boldsymbol{\nabla}^{*}\boldsymbol{u}^{*}-\eta_{a}V^{m}_{a}c^{*}_{a}\boldsymbol{1}\big{)}\quad\text{ in }\Omega_{a},$$

$$\sigma^{*}_{ij}=\lambda_{a}^{*}(c_{a}^{*})\delta_{ij}\frac{\partial u_{k}^{*}}{\partial x_{k}^{*}}+G_{a}^{*}(c_{a}^{*})\bigg{(}\frac{\partial u_{i}^{*}}{\partial x_{j}^{*}}+\frac{\partial u_{j}^{*}}{\partial x_{i}^{*}}\bigg{)}-\eta_{a}V^{m}_{a}c^{*}_{a}(3\lambda_{a}^{*}(c_{a}^{*})+2G_{a}^{*}(c_{a}^{*}))\delta_{ij}\quad\text{ in }\Omega_{a},$$

$${\nabla }^{\ast }\cdot {\mathbf{\sigma }}^{\ast }=0\text{ in }\Omega .$$

$${\mathbf{\sigma }}^{\ast }\cdot {\text{n}}_e=0\text{ on }{\Gamma }_e,$$

$$[{\mathbf{\sigma }}^{\ast }\cdot {\text{n}}_{ab}{]}_{-}^{+}=[{\mathbf{u}}^{\ast }{]}_{-}^{+}=0\text{ on }{\Gamma }_{ab},$$

$${\text{j}}^{\ast }=-\frac{D_a}{R_{g}T}{c}_a^{\ast }{\nabla }^{\ast }{\mu }_a^{\ast }\text{ in }{\Omega }_a,$$

$$c_a^{\ast }{\nabla }^{\ast }{\mu }_a^{\ast }=0\text{ in }{\Omega }_a,$$

$$[{\mu }_a^{\ast }{]}_{-}^{+}=0\text{ on }{\Gamma }_e.$$

$$[{\mu }_a^{\ast }{]}_{-}^{+}=0\text{ on }{\Gamma }_{ab},$$

$$\mu^{*}=\tilde{\mu}_{a}^{\text{SF},*}(c_{a}^{*})-\eta_{a}V_{a}^{m}\text{tr}\big{(}\boldsymbol{\sigma}^{*}\big{)}\quad\text{ in }\Omega_{a},$$

$$c_{0}\sum_{a=1}^{n}\bigg{(}\int_{\Omega_{a}}c_{a}^{\text{max}}\;\;\text{d}V\bigg{)}=\sum_{a=1}^{n}\bigg{(}\int_{\Omega_{a}}c_{a}^{*}\;\;\text{d}V\bigg{)},$$

$${\text{Li}}_{\text{s}}\leftrightharpoons {\text{Li}}_{\text{aq}}^{+}+e_{\text{s}}^{-},$$

$${\tilde{\mu }}_{\text{Li},a}^{\text{SF},\ast }(c_{\text{Li},a}^{\ast })$$

$${\tilde{\mu }}_{\text{Li},c}^{\text{SF},\ast }(c_{\text{Li},c}^{\ast })$$

$$E^{\text{OC}}=-\frac{{\mu }_{e,a}^{\ast }(c_{e,a}^{\ast })}{e}+z{\varphi }_a^{\ast }+\frac{{\mu }_{e,c}^{\ast }(c_{e,c}^{\ast })}{e}-z{\varphi }_c^{\ast }.$$

$${\tilde{\mu }}_{\text{Li},a}^{\text{SF},\ast }(c_{\text{Li},a}^{\ast })=-e{E}^{\text{OC}}+{\mu }_{\text{Li},c}^{\text{SF},\ast }(c_{\text{Li},c}^{\ast }).$$

$${\mathbf{x}}^{\ast }$$

$$c_a^{\ast }=c_a^{\text{max}}{c}_a,$$

$$\nabla \cdot \mathbf{\sigma }$$

$$\mathbf{\sigma }$$

$$\displaystyle=\lambda_{a}(c_{a})\delta_{ij}\frac{\partial u_{k}}{\partial x_{k}}+G_{a}(c_{a})\bigg{(}\frac{\partial u_{i}}{\partial x_{j}}+\frac{\partial u_{j}}{\partial x_{i}}\bigg{)}-\gamma_{a}c_{a}(3\lambda_{a}(c_{a})+2G_{a}(c_{a}))\delta_{ij}\quad\text{ in }\Omega_{a},$$

$$\mu=\mu_{a}^{\text{SF}}(c_{a})-S_{a}^{d}\text{tr}\big{(}\boldsymbol{\sigma}\big{)}\quad\text{ in }\Omega_{a},$$

$${\lambda }_a(c_a)=\frac{{\lambda }_a^{\ast }(c_a^{\ast })}{G_1^{\ast }(0)},G_a(c_a)=\frac{G_a^{\ast }(c_a^{\ast })}{G_1^{\ast }(0)},{\gamma }_a=\frac{{\eta }_a{V}_a^m{c}_a^{\text{max}}}{{\eta }_1{V}_1^m{c}_1^{\text{max}}},\text{ and }{S}_a^d=\frac{{\eta }_a{\eta }_1{V}_a^m{V}_1^m{c}_1^{\text{max}}{G}_1^{\ast }(0)}{R_{g}T}.$$

$$\mathbf{\sigma }\cdot \text{n}$$

$$[\mathbf{\sigma }\cdot \text{n}{]}_{-}^{+}=[\mathbf{u}{]}_{-}^{+}$$

$$c_{0}\sum_{a=1}^{n}\bigg{(}\int_{\Omega_{a}}c_{a}^{\text{max}}\;\;\text{d}V\bigg{)}=\sum_{a=1}^{n}\bigg{(}\frac{c_{a}^{\text{max}}}{c_{1}^{\text{max}}}\int_{\Omega_{a}}c_{a}\;\;\text{d}V\bigg{)}.$$

$$\mathbf{\sigma }$$

$$\mu$$

$$\displaystyle\frac{\text{d}}{\text{d}r}\bigg{(}\frac{1}{r^{2}}\frac{\text{d}}{\text{d}r}(r^{2}u)\bigg{)}-\frac{\gamma_{a}(3\lambda_{a}+2G_{a})}{\lambda_{a}+2G_{a}}\frac{\text{d}c_{a}}{\text{d}r}=0,\quad$$

$$\frac{\text{d}u}{\text{d}r}+\frac{2u}{r}-\frac{{\gamma }_a(3{\lambda }_a+2G_a)}{{\lambda }_a+2G_a}{c}_a={\beta }_a\text{ in }{\Omega }_a,$$

Full text

††journal: Journal of The Electrochemical Society††∗ Corresponding author

Email addresses: i.roper1993@gmail.com (Ian P. E. Roper)

chapman@maths.ox.ac.uk (S. Jon Chapman)

please@maths.ox.ac.uk (Colin P. Please)

1 Introduction

Lithium-ion batteries are ubiquitous in modern technology, powering mobile phones, laptop computers and electric vehicles [48]. Many different materials are being tested for use in all components of lithium-ion batteries to achieve the largest energy capacity while remaining portable and light-weight. A key component for maximising energy capacity is the anode; the more lithium that can be intercalated into the anode, the more charge can be stored, and thus the more energy can be obtained from a single charge. However, upon lithiation many high capacity anodes increase in volume [56], due to changes in the crystal structure. This can cause mechanical stresses because of geometrical constraints within the anode and the rest of the battery, but also because of concentration gradients of lithium within the anode causing non-uniform expansion. The expansion of the anode can displace other components of the battery and the high stresses within the anode material itself can cause it to crack [2]. These effects can result in loss of connection between components within the battery, leading to capacity fade after several charge/discharge cycles.

One anode material with high expansion and a large capacity for lithium atoms which has received a lot of recent research attention is silicon. When fully lithiated, silicon can accommodate 3.75 lithium atoms per silicon atom, forming the alloy Li3.75Si at room temperature [17], giving silicon anodes a volumetric capacity of 3500 mAh/g. This is much greater than the volumetric capacity of 372 mAh/g achieved by commonly-used graphite anodes [25], due to graphite only being able to accommodate one lithium atom per six carbon atoms (LiC6). This makes silicon a very inviting anode material; however, when silicon is lithiated, it expands to around four times its original size [6], much more than the $\sim 10$% volume increase observed in fully lithiated graphite.

Experimental efforts to make silicon a viable anode material often involve using nano-structures within the anode design to minimise the variation in the lithium concentration. These designs are often based on nanoparticles [22, 27, 29, 40], nanowires [7, 18, 54] and nanotubes [41, 45]. Additionally, there are many nano-structure designs which attempts to constrain the expansion of the silicon, for example in a core–shell structure [40, 51, 52], or yolk–shell structure [26, 46], or using self-healing polymers [44]. These structures are often used in conjunction with each other, for example Wang et al. showed that for a silicon oxide shell which is thick enough, inwards growth can be induced in silicon nanotubes [45].

In determining which structures may be the most effective at mitigating the adverse effects of highly expanding anode materials, it is insightful to construct mathematical models for the mechanical behaviour of the anode and the behaviour of the lithium within. Several mathematical models for the expansion of anodes upon lithiation have been proposed. In general, the expansion is assumed to be due to the intercalation of lithium, and opposed by a mechanical response from the stiffness of the anode material, causing stresses. Linear elasticity models [8, 36, 57] have been used to model the mechanical response when strains are small, but more commonly finite strain models using geometrically nonlinear elasticity are used for silicon to capture the large strains that occur [4, 6, 9, 10, 15]. Sethuraman et al. also showed that silicon plastically flows during cycling [38] due to the high stresses induced. This plastic flow has been incorporated into several finite strain models [4, 6, 9], as well as linear models [50]. Several more simple models have neglected the elastic stresses altogether and modelled the mechanical response of the silicon as solely plastic [21, 58]. An additional phenomenon that has been the basis of several modelling studies is crack formation and propagation due to the high stresses within the anode [8, 13, 15, 37, 50].

Two types of model for the lithiation are commonly used: a single phase model in which the silicon is gradually lithiated [4, 6, 10], and two-phase models in which there is a sharp reaction front between the lithiated phase and the unlithiated phase [20, 50, 55]. The former can be simply modelled using diffusion of the lithium atoms, whereas the lithium kinetics in the latter are often described using a Cahn-Hilliard phase-field model. Experimentally, the two-phase model has been shown to likely be more physically accurate [9, 29]. The time-dependence of these models for the evolution of the lithium concentration often cause the main focus of these studies to be on the stresses inside the nano-structures due to silicon non-uniformity, rather than due to material heterogeneity [8, 9, 10, 36].

The coupling of the lithium distribution to the stress is a two-way phenomenon, since in addition to lithiation causing expansion, lithium diffusion is also affected by stress. Stress-assisted diffusion is typically included into lithium transport models by considering the chemical potential of the lithium in the anode to depend on the hydrostatic stress. This coupling affects the diffusion of the lithium through the anode due to the stress induced by lithiation. Several works [4, 8, 11, 14, 36] have included stress-assisted diffusion using the model formulated by Larche and Cahn [23] using a small-strain assumption. Works by Wu [49] and Cui et al. [10] have resulted in a more general framework for incorporating a stress-dependence into the chemical potential of the lithium. However, Cui’s general form of stress-dependent chemical potential reduces to those proposed by Wu, Larche and Cahn under the appropriate assumptions. This general stress-assisted diffusion framework has more recently been applied to finite strain models, for example in [6].

To simplify the mathematical models of these anode nano-structures, simplified geometries are often used. Nano-particles are approximated by spherical geometries [10, 50, 58], nano-wires and nano-rods are approximated to cylindrical geometries [6], and thin films of silicon are assumed to grow unidirectionally [5]. However, even in these simplified geometries, the use of nonlinear elasticity or other complex models requires the use of software packages to solve the resulting differential equations [6, 9, 10]. Despite these systems requiring numerical solvers, the chemical potential of the lithium is often approximated as an analytical expression of the lithium concentration [6, 10, 39]. These approximations fail to capture the large ‘steps’ in chemical potential that occur at transitions in the structure of the intercalated anode [3]. Using an interpolated function of the chemical potential of the lithium in a material would capture these steps and their use in these complicated models would only marginally increase the computational cost.

In this work, we consider a multi-material anode, which we call a hybrid anode. We present a linearly elastic mechanical model coupled to a stress-assisted diffusion model using the stress-dependent chemical potential proposed by Larche and Cahn with a stress-independent term depending on lithium concentration. We focus on a time-scale much slower than the diffusion of lithium through the nano-particle allowing us to use a quasi-static approximation, focusing on the stresses induced by the different expansions of the different materials in chemical and mechanical equilibrium. We present this model as a general chemo-mechanical model for anode lithiation. We then apply our model to a spherical, radially-symmetric geometry consisting of a core and shell of different materials, representing a nano-particle with a constraining shell.

Previous works have often modelled a constrained core–shell anode nano-particle using a zero displacement condition at the shell interface [10, 58] or by modelling the shell as a material that cannot be lithiated [16, 43]. One of our main focuses is the distribution of lithium between the core and the shell and thus we include the lithiation of both materials and the resulting displacement of the interface. We show how the distribution of lithium between the two materials depends heavily on both stress-assisted diffusion and the geometry of the nano-particle.

The lithium distribution predicted by the model is used to calculate the chemical potential of the lithium in the hybrid nano-particle, which we can then convert to an OCV. OCVs can be used to infer the state of charge (SOC) of a battery at a given voltage and are used in battery management systems and other battery models [12, 47]. While experimentally measured OCVs for single material anodes are easily found within the literature [32, 24], experimentally obtaining the OCV curve of a new anode design with multiple materials can be expensive and is very time-consuming. Therefore, a theoretical technique for predicting OCVs of multi-material anodes is of interest to manufacturers and researchers. Sethuraman showed that stresses can significantly affect the OCVs [39], thus the OCV of a hybrid anode is not only dependent on the individual materials but also on the stresses induced in each material by the others at equilibrium. We show the importance of including stress into the lithium concentration model by comparing chemical potentials when including and excluding stress-assisted diffusion.

Lastly, while the modelling of the stresses within lithium-ion anodes has been widely studied, very few of these models discuss optimal designs based on the results of the models. In the second half of this work, we discuss several different performance measures that could be used in defining optimality for the anode design. These are i) the amount of lithium intercalated, ii) the expanded volume of the anode, and iii) the maximum induced stress. We then use these performance measures to suggest three different optimality conditions for the hybrid anode design. We derive these conditions for the silicon core–graphite shell nano-particle design and show how the objective function varies with the volume of the silicon core and find the optimal design in each of the three cases.

The rest of this paper is organised as follows. In Section 2.1 we formulate the chemo-mechanical model we use to describe the lithiation of a lithium-ion anode for a general multi-material anode. In Section 2.2, we describe how the stress-free chemical potential of the lithium in each anode material can be calculated from the material’s OCV. We then simplify the general model to one for a radially symmetric nano-particle with a core of material 1 and a shell of material 2 and solve this model, showing that for certain simple geometries, the concentrations are uniform in each material, even with the inclusion of stress-assisted diffusion. In Section 3, we present the lithium concentrations and chemical potentials in a nano-particle with a silicon core and a graphite shell to show the importance of including stress-assisted diffusion into the lithium concentration model. In Section 4.1 we then discuss the three performance measures we will use to then find an optimal design of the silicon–graphite nano-particles in Section 4.2. Finally in Section 5, we discuss the conclusions that can be drawn from this work and the validity of the results in certain regimes.

2 Mathematical Model

In this section, we develop a mathematical model for the lithiation of a multi-material anode. We begin with a general geometry and arbitrary number of different anode materials. We then apply the model to a spherical, radially symmetric geometry with a core of one material and a shell of another material. The anode is the only battery component which we model in this work; we do not model the electrolyte, the cathode, binder, current collectors or any other components.

2.1 General Formulation

We denote the domain of the entire anode as $\Omega$ and the region occupied by each anode material as $\Omega_{a}$ with $a=1,\dots,n$, where $n$ is the number of different materials in the hybrid anode. We present boundary conditions on the boundary between the anode and the electrolyte $\Gamma_{e}$ and on boundaries between different anode materials $\Gamma_{ab}$ where $a\neq b$ and $a,b\in\{1,\dots,n\}$.

2.1.1 Mechanical Model

We suppose the deformation in the anode is due to the volumetric changes due to lithiation and is counteracted by an elastic response from the material due to its stiffness. For simplicity, we assume the strains are small enough that we may use linear elasticity theory. Similar to previous linearly-elastic studies of anode expansion [8, 57], we write the elastic strain tensor as

[TABLE]

where $\boldsymbol{u}^{*}=(u_{1}^{*},u_{2}^{*},u_{3}^{*})$ is the displacement (m), $c_{a}^{*}$ is the concentration of lithium (mol m*-3*) and $\boldsymbol{\nabla}^{*}$ is the gradient operator. The molar volume of material $a$, measured when zero lithium has been intercalated, is denoted as $V^{m}_{a}$ (m3 mol*-1*) and $\eta_{a}$ is the coefficient of compositional expansion (CCE) of material $a$, which is a measure of the volumetric expansion due to lithiation [42].

We relate the stress to the strain tensor by Hooke’s law,

[TABLE]

where $\mathbb{C}_{a}^{*}$ is the stiffness tensor of material $a$ which can vary with the lithiation state of the material. We assume that each anode material is isotropic, allowing us to write the stiffness tensor in suffix notation as

[TABLE]

where $\lambda^{*}_{a}$ is the first Lamé parameter of material $a$ and $G_{a}^{*}$ is the shear modulus of material $a$ [19], both of which can vary with the lithiation state of the material. We now write (2) as

[TABLE]

which in suffix notation, becomes

[TABLE]

where we use the summation convention for repeated indices. Lastly, the entire anode is in mechanical equilibrium, which is given by

[TABLE]

We prescribe a traction free boundary condition on $\Gamma_{e}$, giving

[TABLE]

where $\textbf{n}_{e}$ is the unit vector normal to the interface between the anode and the electrolyte. We also prescribe continuity of normal stress and displacement between anode materials

[TABLE]

where $\textbf{n}_{ab}$ is the unit vector normal to the interface between anode materials $a$ and $b$.

2.1.2 Chemical Model

We model the movement of lithium within the anode as a diffusive process, with flux $\textbf{j}^{*}$ given by

[TABLE]

where $\mu_{a}^{*}$ is the chemical potential of the lithium intercalated into material $a$, $D_{a}$ is the diffusion coefficient of lithium through material $a$, $R_{g}$ is the ideal gas constant 8.314 J mol*-1* K*-1* and $T$ is the temperature, which we assume to be constant, uniform and equal to 298 K. The time-scale we focus on is much slower than any dynamics in the system, so that we may make the quasi-steady approximation $j^{*}=0$. Thus

[TABLE]

implying the chemical potential of the lithium is uniform in each anode material for all non-trivial concentrations $c_{a}$. At the interface with the electrolyte we have

[TABLE]

Similarly, at the interface between anode materials $\mu^{*}_{a}$ must be continuous:

[TABLE]

and therefore $\mu^{*}_{a}$ takes the same value for $a=1,\dots n$ and we denote this value by $\mu^{*}$.

We define the chemical potential of the lithium intercalated into material $a$ using the stress-dependent chemical potential commonly used with linear elasticity [23], given by

[TABLE]

where, $\text{tr}(\boldsymbol{\sigma}^{*})$ denotes the trace of $\boldsymbol{\sigma}^{*}$ and $\tilde{\mu}_{a}^{\text{SF,*}}(c_{a}^{*})$ is the stress-free chemical potential of lithium in material $a$. We explain how $\tilde{\mu}_{a}^{\text{SF,*}}$ is obtained in Section 2.2.

To close the system, we prescribe the total amount of intercalated lithium. We describe this as a proportion of the maximum amount possible using the state of charge parameter $c_{0}\in[0,1]$, where $c_{0}=0$ denotes no lithium and $c_{0}=1$ denotes a fully lithiated anode. Thus we impose

[TABLE]

where $c_{a}^{\text{max}}$ is the maximum lithium concentration possible in material $a$.

In summary, the full dimensional mechanical model is given by the governing equation (6), the Cauchy stress (4) and the boundary conditions (7)-(8). The full dimensional chemical model is given by the uniform chemical potential (13) and state of charge condition (14). Given a value of $c_{0}\in[0,1]$, the concentration and displacement profile can be calculated from the model.

2.2 Using OCVs to Calculate Chemical Potential

In this section, we outline how the stress-free chemical potential of the lithium atoms $\tilde{\mu}_{a}^{\text{SF},*}(c_{a}^{*})$ in a single anode material $a$ can be determined from the OCV of that material. We follow Bazant [1] and Newman et al. [30] by equating the electrochemical potentials of the reactants and the products of the surface reaction occurring at $\Gamma_{e}$. The reversible reaction at the surface of each electrode during charging or discharging is

[TABLE]

where the subscript s denotes that the lithium atoms and electrons are in the solid electrode and the subscript aq denotes that the lithium ions are dissolved in the electrolyte.

As these reactions are in equilibrium during OCV measurement, we may equate the total electrochemical potential of the reactants and products. Since we are only concerned with single-material electrodes here, and the electrodes are assumed to be unconfined, there is no stress-dependent contribution to the chemical potential. The balances of electrochemical potentials at the surface of the anode and the cathode are thus given by

[TABLE]

respectively, where, $\phi^{*}$ represents the electrical potential. The first subscript of the chemical potentials $\mu^{*}$ in (16)-(17) denotes the chemical species that the chemical potential corresponds to. The subscript of the electrical potentials and the second subscript of the chemical potentials denote the phase which the species is in, $a$ denoting anode, $c$ denoting cathode and el denoting electrolyte. We have used the first subscript to avoid confusion with specifying the chemical potentials of different species, however, it should be noted that $\tilde{\mu}^{\text{SF},*}_{a}(c_{a}^{*})$ in the notation of Section 2.1.2 is $\tilde{\mu}^{\text{SF},*}_{\text{Li},a}(c_{\text{Li},a}^{*})$ in (16). Finally, the charge of the Li+ ions is denoted by $z$ and $e$ is the fundamental charge 1.60217$\times 10^{-19}$C.

The OCV, denoted by $E^{\text{OC}}$, measures the difference in the electrochemical potential of the electrons in the anode and the cathode, divided by the charge of the electron $-e$ [30], and hence is given by

[TABLE]

We can thus subtract (17) from (16), substitute (18), and rearrange to obtain an expression for the chemical potential of the lithium in the anode in terms of the OCV, given by

[TABLE]

Hence, we convert the OCV, $E^{\text{OC}}$, conventionally measured in volts, to chemical potential $\tilde{\mu}^{\text{SF},*}_{\text{Li},a}$, measured in J mol*-1*, by multiplying by the fundamental charge, $e$ (J V*-1*), and Avagadro’s number, 6.02214086 $\times 10^{23}$ mol*-1*.

The chemical potential for the lithium in the cathode $\mu^{\text{SF},*}_{\text{Li},c}(c_{\text{Li},c}^{*})$ in (19) does not concern this work for two reasons. Firstly, the OCV of the anode material is often measured against a cathode of lithium metal, thus the backwards reaction of (15) is lithium plating and so $\mu^{\text{SF},*}_{\text{Li},c}$ is independent of $c_{\text{Li},c}^{*}$. This makes $\mu^{\text{SF},*}_{\text{Li},c}$ an additive constant, only shifting the reference potential for $\tilde{\mu}^{\text{SF},*}_{\text{Li},a}$. Secondly, we are only concerned with comparing the potentials of the lithium in different anode materials or comparing the potentials of lithium in an anode material at different states of charge. Therefore, as long as the cathode used when calculating the OCV of the different materials is the same, this additive constant will cancel when comparing them.

2.3 Nondimensionalisation

We now nondimensionalise the variables in our chemo-mechanical model. The spatial coordinate $\boldsymbol{x}^{*}$, the displacement vector $\boldsymbol{u}^{*}$, the Cauchy stress $\boldsymbol{\sigma}^{*}$, the concentration in each material $c_{a}^{*}$ and the stress-free and stress-dependent chemical potentials $\mu_{a}^{\text{SF,*}}$ and $\mu_{a}^{*}$ are nondimensionalised by setting

[TABLE]

respectively. Here, $L$ is the representative length-scale of the anode particle, $G_{1}(0)$ is the shear modulus of material 1 at zero lithiation and we have chosen to scale each variable using the parameters relevant to material 1. We note that the linear forms of the stress (4) and the chemical potential (13) are only justified in situations where $\eta_{a}V_{a}^{m}c_{a}^{\text{max}}\ll 1$ for $a=1,\dots,n$, since then the linearisation of the elastic strain tensor in (1) is permitted. A full derivation from the fully nonlinear formulation to the mechanical model we use here can be found in [35].

Using this nondimensionalisation, the governing equation (6) is

[TABLE]

with

[TABLE]

while the uniform chemical potential is given by

[TABLE]

where the dimensionless parameters are given by

[TABLE]

Henceforth, we drop the explicitly written lithiation-dependence of the Lamé parameters and write $\lambda_{a}\equiv\lambda(c_{a})$ and $G_{a}\equiv G_{a}(c_{a})$. The boundary conditions are

[TABLE]

while the SOC condition is

[TABLE]

2.4 Spherical Geometry

We now focus on the specific geometry of a single, spherical, radially symmetric nano-particle comprising a core of anode material 1 enclosed in a shell of anode material 2. We denote the radius of the core made of material 1 by $R_{1}$, and thus material 1 occupies the domain $\Omega_{1}=\{r|0\leq r<R_{1}\}$. We denote the radius of the entire nano-particle by $R_{2}$ so that the shell made of material 2 has thickness $R_{2}-R_{1}$ and occupies the domain $\Omega_{2}=\{r|R_{1}<r\leq R_{2}\}$. This geometry is shown schematically in Figure 2. Our representative length-scale of the anode particle $L$ is now given by the radius, $R_{2}$.

We assume the displacement $\boldsymbol{u}$ is solely radial due to symmetry and thus we write $\boldsymbol{u}=u(r)\boldsymbol{e}_{r}$. The Cauchy stress and the chemical potential of lithium in material $a$, (22) and (23), are then given by

[TABLE]

respectively. The mechanical equilibrium equation (21) becomes

[TABLE]

Integrating (30) gives

[TABLE]

where the integration constant $\beta_{a}$ may be different in the two materials. The trace of the Cauchy stress tensor (28) is given by

[TABLE]

Using (31) gives

[TABLE]

Thus the trace of the Cauchy stress and therefore the chemical potential, (29), can be written in terms of $c_{a}$ only. Since $\mu_{a}$ is uniform in each material we conclude that $c_{a}$ is also uniform in each material.

Substituting a uniform $c_{a}$ into (31) and solving the resulting ODE yields

[TABLE]

where $A_{a}=3\beta_{a}+(3\gamma_{a}(3\lambda_{a}+2G_{a})c_{a})/(\lambda_{a}+2G_{a})$ is now the constant arising from the first integration, $B_{a}$ is the constant arising from the second integration and each integration constant may be different in the two materials. Substituting (34) into (29) gives

[TABLE]

We can substitute (34) into (28) and write the radial stress and hoop stresses as

[TABLE]

respectively. It remains to use the boundary conditions to determine the integration constants $A_{a}$ and $B_{a}$.

The boundary conditions (25) and (26) for the spherical geometry are written as

[TABLE]

respectively, where $R=R_{1}/R_{2}$. Additionally, we ensure that the displacement is bounded at the origin by prescribing

[TABLE]

Substituting (34) and (36) into (38)-(40) gives

[TABLE]

where

[TABLE]

Thus the displacements and stress are given once the concentrations $c_{1}$ and $c_{2}$ have been determined.

The SOC condition (27) becomes

[TABLE]

which we can rearrange as

[TABLE]

Finally, we use (35) to relate the concentrations $c_{1}$ and $c_{2}$. Combining this with (47) gives a single algebraic equation in $c_{2}$ for a given $c_{0}$ to be solved numerically222For very high values of $R$ and very low values of $c_{0}$, this function has multiple roots. For consistency, we take the lowest root in these cases.. As the concentrations $c_{a}$ must be between zero and one due to the nondimensionalisation (20), we place bounds on $c_{2}$, ensuring that both $c_{1},c_{2}\in[0,1]$. We use (47) to find that

[TABLE]

where the first and second elements of the maximum and minimum functions ensure $c_{2}\in[0,1]$ and $c_{1}\in[0,1]$, respectively.

3 Chemical Results

We illustrate the behaviour of a nano-particle consisting of a silicon core and a graphite shell described by the model and geometry explained in Section 2.4. The relevant mechanical and chemical parameters for these materials are given in Table 1, including the nondimensional parameters $\gamma_{a}$ and $S^{d}_{a}$. We calculate the Lamé parameters of silicon and graphite from the Young’s modulus and the Poisson’s ratio. The MATLAB code to produce the figures in this section is available publicly on GitHub [34].

We assume the Young’s modulus varies linearly with the lithiation state using the values at zero and full lithiation as reference values. Therefore, the dimensional Lamé parameters are given by

[TABLE]

where $E^{0*}_{a}$ is the Young’s modulus of material $a$ at zero lithiation, $\nu_{a}$ is the Poisson’s ratio of material $a$ and $\eta^{E}_{a}$ is the linear variation of the Young’s modulus of material $a$. As the Poisson’s ratios of these materials vary negligibly with lithiation state [10], we assume $\nu_{a}$ to be independent of $c_{a}$.

We plot the non-dimensionalised stress-free chemical potentials of the lithium intercalated into silicon and graphite, $\mu_{1}^{\text{SF}}$ and $\mu_{2}^{\text{SF}}$, in Figure 2. It can be seen that there are several regions where a small change in lithium concentration causes a large change in the chemical potential, as noted in [3].

In Figure 3 we plot the lithium concentration in each material, $c_{1}$ and $c_{2}$, against $c_{0}$ for different values of the volume fraction $\psi=R^{3}$ of the silicon core. It can be seen that the graphite is saturated ($c_{2}=1$) at fairly low $c_{0}$ values, especially for large silicon cores, whereas the lithium concentration in the silicon, $c_{1}$, remains relatively low until the graphite is saturated. After this point, all the lithium being intercalated must be intercalated into the silicon and thus $c_{1}$ is linear in $c_{0}$.

The saturation of the graphite, and the low $c_{1}$ values, at low SOC are due to the hydrostatic stresses induced in the silicon core and graphite shell. The trace of the Cauchy stress in each material, given by (33), is plotted against $c_{0}$ in Figure 4 for different silicon volume fractions, $\psi$. Since $\gamma_{1}\gg\gamma_{2}$, the expansion of the silicon is much greater than that of the graphite. Therefore, when lithiated, the expansion of the silicon is being constrained by the graphite shell, inducing a compressive stress in the silicon, causing $\text{tr}(\boldsymbol{\sigma})$ to be negative in $\Omega_{1}$. Conversely, the graphite is being stretched by the large expansion of silicon, inducing a tensile stress, causing $\text{tr}(\boldsymbol{\sigma})$ to be positive in $\Omega_{2}$. The tensile stress in the graphite lowers the chemical potential of the lithium whereas the compressive stress in the silicon increases the chemical potential of the lithium. Therefore, for the chemical potentials of the lithium in each material to be equal, the stress-free potential of the lithium in the silicon must therefore decrease and the stress-free potential of the lithium in the graphite must increase. The stress-free potential of the lithium in each material is a monotonically increasing function of the concentration, therefore $c_{1}$ must decrease and $c_{2}$ must increase.

The low $c_{1}$ values are present for a larger range of SOC for small silicon cores than for large ones since a thicker graphite shell is stronger, inducing a larger compressive stress on the silicon. Likewise, a small silicon core does not stretch the graphite shell as much and thus, for example, the $c_{2}$ plot for $\psi=0.05$ is more similar to the single material case than that for $\psi=0.5$.

In Figure 5 we plot the chemical potential of the lithium in the hybrid nano-particle against SOC for different values of $\psi$. It can be seen that for SOC values for which the graphite is saturated, the chemical potential is much greater than that of the individual materials. This is because if the graphite is saturated but more lithium is intercalated, the lithium concentration in the silicon must increase. This increases $\mu_{a}^{\text{SF}}$ but also decreases $\text{tr}(\boldsymbol{\sigma})$ further as seen in Figure 4, causing the chemical potential to increase rapidly with increasing SOC.

In Figures 6 and 7 we replicate Figures 3 and 5, respectively, but we do not account for the effect of stress-assisted diffusion (we set $S_{a}^{d}=0$) to show the importance of including the stress of the anode materials in the chemical potential of the lithium. It can be seen in Figure 6 that $c_{1}>c_{2}$ for $c_{0}\lessapprox 0.8$ when $S^{d}_{a}=0$, which is very different to the result found in Figure 3 in which $c_{2}>c_{1}$ for $c_{0}\gtrapprox 0.05$. This is solely due to the chemical potential of the lithium in the silicon being less than that for graphite for $c_{a}\lessapprox 0.8$, seen in Figure 2. In Figure 7, the chemical potential of the lithium in the hybrid nano-particle is approximately an interpolation between the chemical potentials of the lithium within the individual materials, which is also very different to the case in which $S_{a}^{d}\neq 0$.

The large differences in $c_{1}$ and $c_{2}$ between Figures 3 and 6 caused by the exclusion of stress-coupling demonstrates how expansion due to lithium will be dramatically affected by this mechanism. Furthermore, the large differences in the chemical potential of the lithium in the hybrid particle between Figures 5 and 7 show that stress-modelling is essential in predicting OCVs of hybrid anodes as a function of SOC.

4 Optimal Size of the Silicon Core

4.1 Performance Measures

We now discuss the performance measures that might be considered to find an optimal anode design, applying them to the nano-particle geometry shown in Figure 1. The three properties we consider as these performance measures are i) the total amount of lithium absorbed, ii) the expanded volume of the anode and iii) the maximum stress induced. We first derive each of these measures using the results of our model and then investigate how these might be applied in practice to determine the optimal volume fraction of the silicon core in the spherical core–shell nano-particle geometry.

4.1.1 Amount of Lithium

The capacity of a lithium-ion battery is closely related to the amount of lithium that the anode can accommodate and thus this is a very important performance measure of a lithium-ion battery anode. For a general anode with domain $\Omega$ and different anode materials $i=1,\dots,n$, each with domain $\Omega_{i}$, we measure the total amount of intercalated lithium relative to a fully lithiated anode of material 1 and domain $\Omega$. We denote this relative amount of lithium as

[TABLE]

For the core–shell geometry in Figure 1, $Q$ is given by

[TABLE]

where we recall $R=R_{1}/R_{2}$.

4.1.2 Relative Expanded Volume

The expansion of the anode can have serious adverse effects on the battery performance and so we consider this expansion by calculating the volume of the expanded nano-particle compared to the original volume, denoting this as the relative expanded volume, $V$. For a general geometry of anode, $V$ would usually have to be calculated numerically, however, for the spherical nano-particle geometry in Figure 1, $V$ is given by

[TABLE]

Here, we use the dimensional displacement $u^{*}(R_{2})=\eta_{1}V_{1}^{m}c_{1}^{\text{max}}R_{2}u(1)$ as in (20).

In Figure 8, we plot $V$ against $\psi$ for different states of charge. It can be seen that $V$ increases with the volume of the silicon core with this effect becoming more prominent for larger $\psi$. Additionally, a greater SOC causes a greater $V$ for all volumes of silicon core. The calculated $V$ value for a fully lithiated nano-particle with a large silicon core, is much greater than the observed expanded volume for a solely silicon nano-particle, given by $V=J^{c}_{1}=3.8$. This over expansion for large silicon cores is due to our adoption of linear elasticity in Section 2.1.1. This assumes that $\eta_{1}V_{1}^{m}c_{1}^{\text{max}}\ll 1$, whereas for silicon, $\eta_{1}V_{1}^{m}c_{1}^{\text{max}}=0.933375\sim 1$. Thus the nonlinear elasticity formulation should be retained to produce quantitatively accurate results for silicon at high concentrations.

4.1.3 Maximum Induced Stress

Capacity fade exhibited by expanding anode materials after cycling is often attributed to cracks in the anode material caused by high stresses. Thus the final performance measure we use to optimise anode performance is the maximum induced stress. We use the von Mises stress [19] as a scalar effective stress measure. For a general anode geometry, this is given by

[TABLE]

In radial symmetry, the von Mises stress (53) can be written as

[TABLE]

so that from (20), (36) and (37), we have

[TABLE]

for $a=1,2$. The condition at $r=0$ (40) gives $B_{1}=0$. Therefore, the effective stress is only non-zero in $\Omega_{2}$. The $r^{-3}$ dependence shows that this effective stress is the greatest at the minimum value of $r$ in $\Omega_{2}$, which is $r=R=R_{1}/R_{2}$. Therefore, the maximum induced stress is given by

[TABLE]

In Figure 9, we plot $\sigma_{\text{eff}}^{*}(R_{1})$ against $\psi$ for different SOC. For all states of charge, $\sigma_{\text{eff}}^{*}(R_{1})$ is minimised as $\psi\rightarrow 0$; therefore, a very small silicon core is optimal. Of course, a single-material anode nano-particle will not induce any stress in equilibrium and thus the two single-material designs in this case are globally optimal for this measure with $\sigma_{\text{eff}}^{*}(R_{1})=0$.

4.2 Optimisation

We now suggest objective functions and constraints that could be used to optimise the design of an anode using the three performance measures defined above. We derive these objective functions and constraints for the core–shell spherical nano-particle shown in Figure 1, and plot how the objective function varies with the volume of the silicon core with the aim of finding the optimal size of silicon core for each optimisation problem.

4.2.1 Amount of Lithium per Relative Expanded Volume

The performance of anode materials are often measured by their volumetric or gravimetric capacity. The first measure we use to optimise the anode geometry is the amount of lithium per expanded volume which is equivalent to the volumetric capacity. We calculate the amount of lithium per expanded volume by dividing the relative amount of lithium, $Q$ in (51), by the relative expanded volume, $V$ in (52), giving

[TABLE]

In Figure 10, we plot $Q/V$ against $\psi$ for the same states of charge as in Figures 8-9. We can see from this plot that $Q/V$ is maximised by a fully lithiated nano-particle with a silicon core of volume fraction $\approx 0.45$. The reason that $Q/V$ is not optimised by a nano-particle of a single material is that the amount of lithium $Q$ is linear in $\psi$ ($=R^{3}$), whereas the expanded volume $V$ increases much more rapidly with $\psi$ at large values of $\psi$ than for smaller values. Therefore, as $\psi$ becomes larger, the extra lithium that can be intercalated is out-weighed by the increased expansion and $Q/V$ begins to decrease with $\psi$. However, as explained in Section 4.1.2, the linear elasticity model overestimates $V$ at high SOC and high $\psi$ so this prediction needs to be validated with a nonlinear model. For small SOC values, $Q/V$ has a much more linear relationship with $\psi$, and thus if the nano-particle is only lithiated a small amount, a fully silicon nano-particle is the optimal design according to the $Q/V$ measure. Lastly, we observe that for large $\psi$, $Q/V$ is not monotonic with $c_{0}$ and so having a partially lithiated nano-particle yields a higher volumetric capacity than a fully lithiated one.

4.2.2 Maximising Amount of Lithium Constrained by Maximum Expanded Volume

We now consider the problem of maximising the amount of lithium the anode can hold subject to constraints on the expanded volume and the maximum effective stress. We begin with the expanded volume constraint. Thus, we want to maximise $Q$, given by (50), subject to $V\leq V_{\text{max}}$, for some prescribed $V_{\text{max}}$.

For a given $\psi$, $Q$ is maximised by maximising the SOC, $c_{0}$. From Figure 8, we see that for a given $\psi$ value, the expansion $V$ is monotonically increasing with $c_{0}$. Thus, the maximum viable $c_{0}$ value, which we denote as $\check{c}_{0}$ occurs when $V=V_{\text{max}}$, unless the fully lithiated volume is less than the constraint, in which case $\check{c}_{0}=1$. Therefore, for each $\psi\in[0,1]$, we find $\check{c}_{0}$ such that

[TABLE]

We then maximise the corresponding $Q_{\text{max}}=Q(\check{c}_{0})$ over $\psi$.

The value of $\psi$ at which $V_{c_{0}=\check{c}_{0}}$ changes from $V_{c_{0}=1}$ to $V_{\text{max}}$ in (58) can be calculated by substituting (34) into (52) for $c_{1}=c_{2}=1$ and equating to the prescribed $V_{\text{max}}$, giving

[TABLE]

where $\hat{A}_{2}$ and $\hat{B}_{2}$ are given by (43) and (44) with $c_{1}=c_{2}=1$. This can be solved to find a critical value of $R$, which we denote $\hat{R}$, which corresponds to the $\psi$ value at which $V_{c_{0}=1}=V_{\text{max}}$. For all values of $\psi$ such that $R<\hat{R}$, $Q_{\text{max}}$ is given by (51) with $c_{1}=c_{2}=1$ and $R=\psi^{1/3}$. However, for $R>\hat{R}$, we must numerically solve

[TABLE]

to find $\check{c}_{0}$ and $Q_{\text{max}}$. We both derive the expression for $\hat{R}$ and the bounds on $V_{\text{max}}$ such that $0<\hat{R}<1$ in A.

In Figure 11 we plot $Q_{\text{max}}$ against $\psi$ for four different values of $V_{\text{max}}$. The $\psi$ values for which $R<\hat{R}$ can be seen by the straight line on the left side of the plot. This is because for $c_{1}=c_{2}=1$, $Q$ is a linear function of $R^{3}=\psi$ in (51). Immediately to the right of this region, $Q_{\text{max}}$ decreases for $R>\hat{R}$. The $\psi$ value which gives the largest $Q_{\text{max}}$ is $\psi=\hat{R}^{3}$ as can be seen by the peaks in $Q_{\text{max}}$ after the linear region for small $\psi$.

4.2.3 Maximising Amount of Lithium Constrained by Maximum Stress

We now wish to maximise the amount of intercalated lithium under a constraint on the maximum induced stress, defined in (53). We wish to find the anode geometry which maximises (50) subject to the maximum induced stress in the anode (56) being less than some prescribed maximum $\sigma_{\text{max}}$.

As with the expanded volume constraint in Section 4.2.2, for each volume fraction of silicon, $\psi$, we must find the maximum $c_{0}$ such that the maximum induced stress $\sigma_{\text{eff}}^{*}(R_{1})$ is less than $\sigma_{\text{max}}$. In Figure 12, we plot $\sigma_{\text{eff}}^{*}(R_{1})$ against SOC and it can be seen that $\sigma_{\text{eff}}^{*}(R_{1})$ is a monotonically increasing function of $c_{0}$. Therefore, the maximal $c_{0}$ value, denoted by $\check{c}_{0}$, under the constraint $\sigma_{\text{eff}}^{*}(R_{1})<\sigma_{\text{max}}$ will be such that $\sigma_{\text{eff}}^{*}(R_{1})=\sigma_{\text{max}}$. We then calculate the maximum amount of lithium, $Q_{\text{max}}=Q(\check{c}_{0})$, for that $\psi$ value by substituting $\check{c}_{0}$ into (51).

In Figure 13, we plot $Q_{\text{max}}$ against $\psi$ for several different values of $\sigma_{\text{max}}$. It can be seen that $\psi=1$ gives the greatest $Q_{\text{max}}$ value for each value of $\sigma_{\text{max}}$. However, for $\sigma_{\text{max}}=4.0$ and 5.0 GPa, $Q_{\text{max}}$ is not monotonic in $\psi$, thus there are local maxima of $Q_{\text{max}}$. This is useful for design purposes given other constraints. For example if the maximum stress that is allowed is $\sigma_{\text{max}}=4.0$ GPa, but we must also restrict $\psi$ to be less than 0.25, having $\psi=0.05$ would allow a larger amount of lithium to be intercalated than with $\psi=0.2$.

As in Section 4.2.2, for sufficiently large $\sigma_{\text{max}}$, there are $\psi$ values for which a fully lithiated nano-particle ($\check{c}_{0}=1$) does not produce stresses as large as $\sigma_{\text{max}}$. Therefore, the the linear relationship for small $\psi$ in Figure 11 is also seen for large values of $\sigma_{\text{max}}$. However, the $\sigma_{\text{max}}$ values illustrated in Figure 13 are too small to see this behaviour; all of the curves there have $\check{c}_{0}<1$. We give the value of $\hat{R}$ (the core radius for which the maximum induced stress at $c_{0}=1$ is $\sigma_{\text{max}}$) and the bounds on $\sigma_{\text{max}}$ such that $0<\hat{R}<1$ (thus the linear relationship for small $\psi$ is observed) in A.

The $\sigma_{\text{max}}$ values plotted for here are much greater than the tensile strength of graphite, $\sigma_{\text{C}}^{Y}$, typically between 8 and 12 MPa [28]. The results for the constraint $\sigma_{\text{max}}=\sigma_{\text{C}}^{Y}$ are qualitatively similar to that of $\sigma_{\text{max}}=0.5$ GPa. Therefore, to avoid the yielding or cracking of the graphite at the interface, the optimal design is as large a silicon core as possible. For $\psi=0.99$, the silicon can only be lithiated to around $c_{1}=2.2\times 10^{-4}$ for the interfacial stress (56) to remain under the yield stress, while the graphite remains unlithiated. This achieves a $Q_{\text{max}}$ value also of around $2.2\times 10^{-4}$ (as $c_{2}=0$ in this case).

5 Conclusions and Discussion

In this paper, we present a model for the inclusion of mechanical stress in calculating lithium distribution for multi-material lithium-ion batteries. We show how the stress can be included into the concentration model for simple geometries; for example spheres, cylinders and 1-dimensional plates, and determine that the lithium concentration is uniform in each material when diffusion is fast compared to the charging rate. By applying the model to a radially symmetric spherical nano-particle with a silicon core and a graphite shell, we show there are large changes in lithium distribution if the chemo-mechanical coupling effects are neglected. Finally, we present a framework for finding the optimal design for the geometry of multi-material anodes and present three performance measures and three optimality conditions. We present results from optimising the core volume for a silicon core, graphite shell geometry as an example of the insight that can be obtained from using these measures.

There are several limitations to the static linear elasticity model we derive in Section 2. Firstly, while the quasistatic assumption is valid for small currents, the non-uniformity in the lithium concentration caused by diffusion through the anode during more strenuous battery function causes substantial stresses that are not captured by this model. Secondly, there are mechanical phenomena commonly observed experimentally, such as plasticity and cracking, that are also not included in this model. Our solely elastic mechanical model assumes that the materials are able to act elastically regardless of the stress induced by the lithiation and there is no yield stress at which the material either begins to plasticise or crack. Lastly, as noted in Section 2.3, the linear elasticity assumptions used to derive the mechanical model in Section 2.1.1 relies on the parameter $\eta_{a}V_{a}^{m}c_{a}^{\text{max}}\ll 1$ for $a=1,\dots,n$, thus restricting the materials this linear model is applicable to.

We also make several assumptions about the materials in our model which are significant simplifications for many anode materials. We assume the anode materials are isotropic to allow us to write the stiffness tensor in terms of two Lamé parameters, $\lambda$ and $G$, and this ignores the anisotropy of the crystal structure of the anode. While we have taken the lithium-concentration-dependence of the Young’s modulus into account, we still assume that the Poisson’s ratio and molar density are lithium-concentration-independent. Finally, we use parameters and a model suitable for bulk materials, while we are imagining nano-structures in our example. It has been experimentally shown that nano-sized materials have very different mechanical properties from bulk materials [59] and accounting for surface effects in nano-structures becomes increasingly important as the size of the material decreases [33].

We use the example of a silicon core and a graphite shell to show the results in Figures 3-13; however, for silicon, $\eta_{1}V_{1}^{m}c_{1}^{\text{max}}=0.933375$. Therefore, our application of this model to a silicon core surrounded by a graphite shell needs to be interpreted with some care. The expanded volume of the hybrid particle with a silicon core and graphite shell, shown in Figure 8, is larger than the expected values because of the linear assumption being violated. Additionally, the disparity between $\eta_{1}$ and $\eta_{2}$ in Table 1 means that the graphite shell is being stretched by a large amount, causing very large tensile stresses. As we have neglected cracking and plasticity from the model, this means the tensile stresses of the graphite are most likely exaggerated in the results shown here, which will not only affect the optimisation results in Section 4 but also the lithium distributions within the anode. Although the model is somewhat outside its range of validity when applied to silicon anodes, we use silicon as material 1 in our example to clearly show the necessity of including stress into the static lithium concentration model due to its large expansion. This model can still be applicable to silicon if we were to restrict the lithiation of the silicon such that $V^{m}_{1}c_{1}^{\text{max}}<3.75$, thus the results shown here are still accurate for small $c_{1}$ levels.

Despite the limitations of the model, the simplifications made in this work allow us to analytically show the consequences of including a stress-dependence into the chemical potential. One key result of this is to show that stress is induced solely by the presence of different materials, independently of the non-uniformity of the lithium concentration due to diffusion. While the differences in stress-dependent lithium concentrations and stress-independent concentrations will be less pronounced with materials with lower expansions than silicon, these stresses are still important to include for accurate OCV prediction and stress modelling. While several anode materials have non-isotropic crystal structures, causing the isotropy assumption to be inappropriate, the crystalline structure of anode materials usually become amorphous after the first charge, especially if the stresses are large. Additionally, the introduction of $\lambda$ and $G$ does not hugely simplify the model for general geometries and so removing this assumption for these cases would only slightly increase the complexity of the model. Finally, the simplifications made to the model have allowed to provide analytical insight into the key performance indicators of anode nano-particles that we define in Section 4. Therefore, this model can be used for insight into optimal anode design for lithium-ion batteries, but it must be adapted for more quantitatively accurate results and for more challenging design geometries to be used in practice.

6 Acknowledgements

This publication is based on work supported by the EPSRC Centre For Doctoral Training in Industrially Focused Mathematical Modelling (EP/L015803/1) in collaboration with Nexeon. Special acknowledgement goes to the industrial collaborator, Bill Macklin, for very insightful discussions and motivations for this work. JC and CP acknowledge funding from the Faraday Institution (EP/S003053/1).

Appendix A Calculation of $\hat{R}$ for Maximum Expansion and Maximum Stress Constraints

We present expressions for the critical radii $\hat{R}$ for the maximum expansion constraint and the maximum stress constraint where $\hat{R}$ is the radius of the silicon core such that the constraint is met when $c_{0}=1$. We also present the inequalities that the constraints $V_{\text{max}}$ and $\sigma_{\text{max}}$ must satisfy such that $0<\hat{R}<1$.

A.1 Maximum Expansion Constraint

Substituting (43)-(44) into (59) yields

[TABLE]

where $\Lambda_{a}=3\lambda_{a}+2G_{a}$ for $a=1,2$ and $\bar{\eta}_{1}=\eta_{1}V_{1}^{m}c_{1}^{\text{max}}$. This can be rearranged to give

[TABLE]

We rearrange the inequality $0<\hat{R}<1$ with $\hat{R}$ given by (62) to find

[TABLE]

A.2 Maximum Stress Constraint

Substituting (44) into (56) with $c_{1}=c_{2}=c_{0}=1$ yields

[TABLE]

which can be rearranged to give

[TABLE]

We rearrange $0<\hat{R}<1$ with $\hat{R}$ now given by (65) to find

[TABLE]

must be satisfied to satisfy $0<\hat{R}<1$.