# Influence of constraints on axial growth reduction of cylindrical Li-ion battery electrode particles

###### Abstract

Volumetric expansion of silicon anode particles in a lithium-ion battery during charging may lead to the generation of undesirable internal stresses. For a cylindrical particle such growth may also lead to failure by buckling if the expansion is constrained in the axial direction due to other particles or supporting structures. To mitigate this problem, the possibility of reducing axial growth is investigated theoretically by studying simple modifications of the solid cylinder geometry. First, an annular cylinder is considered with lithiation either from the inside or from the outside. In both cases, the reduction of axial growth is not found to be significant. Next, explicit physical constraints are studied by addition of a non-growing elasto-plastic material: first, an outer annular constraint on a solid silicon cylinder, and second a rod-like inner constraint for an annular silicon cylinder. In both cases, it is found that axial growth can be reduced if the yield stress of the constraining material is significantly higher than that of silicon and/or the thickness of the constraint is relatively high. Phase diagrams are presented for both the outer and the inner constraint cases to identify desirable operating zones. Finally, to interpret the phase diagrams and isolate the key physical principles two different simplified models are presented and are shown to recover important qualitative trends of the numerical simulation results.

## 1 Introduction

The lithium-ion battery (LIB) has established itself as the power source of choice for mobile phones, laptop computers, and a variety of hand-held electronic devices [1, 2, 3]. The reason for this wide-spread use of LIBs is that they are lightweight (primarily because lithium is the lightest metal) and they have large storage capacity (because they have relatively high energy density). Motivated by this success, LIBs have been put forward, for over a decade, as the best candidates for use in electric transportation systems. Such a vision has materialized in the recent years with the launch of a number of electric-hybrid and electric vehicles. While the portable electronic devices need just one or, at most, a few unit cells in a LIB, the battery pack in a single electric car may use as many as 6000 cells [4]. The development of these cars has, therefore, resulted in a tremendous increase in the demand for LIBs. Yet, currently, wide-spread use of electric vehicles is limited by the energy capacity limits.

The encouraging, albeit limited, success with electric cars in recent years has fostered hope that the automotive industry can indeed be driven towards a paradigm shift from a singular dependence on fossil fuels to a cleaner power source like LIBs. However, if this initial success is to be sustained and even the conservative predictions of growth trajectory for wide-spread commercial deployment are to be met, then LIBs need to have significantly higher energy capacities than current commecial ones. Only then prices may be brought down by using fewer cells in battery packs, and mileages may be improved – thus requiring fewer charging stations. For these developments to actually happen, the current materials in commercial LIBs need to be replaced with those that, through a fundamentally different chemistry, can provide higher limits on lithium storage. For the anode, the best candidate to emerge with such a property is silicon.

In the fully lithiated state, a single atom of Si can accommodate up to 4.4 atoms of Li resulting in the equilibrium (amorphous) phase, LiSi [5]. This is a significant improvement on the traditional anodic material, graphite, which gives LiC in the fully lithiated state. This translates into a theoretical specific capacity value of 4200 mAhg compared to only 372 mAhg for graphite. A higher energy capacity follows directly from this higher specific capacity.

Despite the exciting promise of such new chemistry, the use of silicon in a LIB comes with its own big challenge: lithiation of silicon results in significant volume change, which can be as high as 310% in the fully lithiated state [6]. Slow diffusion of Li in Si results in a spatially non-uniform Li concentration, leading to stress generation due to differential growth. A similar situation arises during delithation. Indeed, it has been found that cyclic charging and discharging causes pulverization of the Si anode particles which ultimately results in a decay of the specific capacity [7, 8, 9]. This problem may be overcome to some extent by using nanostructured Si anode particles [8, 10, 11, 12, 13, 14, 15].

Nevertheless, even with nanostructured particles, practical challenges persist. Although the Si anode particles do not pulverize, they can still grow. This growth may be problematic within the finite confines of the electrode (and, by extension, the whole battery) packaging. While it might seem that the immediate solution to this problem would be to allow for some “free” space within the packaging to allow for such volumetric expansion, in practice, this approach is untenable because it leads to loss of electrical contact within the electrode, and, thus to a severely impaired battery performance. Furthermore, these anode particles will be in physical contact with other particles and supporting substrates. Expansions in the presence of such constraints will lead to stresses and, possibly, mechanical failure even without fracture. Such a situation was investigated in [16] where we studied the possibility of a cylindrical Si anode particle failing through buckling (under axial confinement). Left unconstrained, however, such a cylindrical particle grows both axially and radially during lithiation. The percentage increase in length of such a cylinder with increasing state of charge, at different charging rates, is shown in Fig. 1. It is this increase in length which provides the motivation for our present study.

We observe in Fig. 1 that the percentage increases in length for three different influx rates (represented by ) are large and nearly equal to unconstrained, isotropic growth. Our objective here is to investigate the simplest modifications in geometry on a cylinder which might reduce axial growth. We first study the case of an annular cylinder to check if introducing a “free space” within the particle might be useful for limiting axial growth. Next, guided by the observation that a slight spatial non-uniformity in Li concentration as seen for higher influx rates leads to a decrease in the axial growth, we impose two simple constraints in the radial direction. The question we ask ourselves is: If spatial variations in lithiation can reduce the axial growth through the influence of “saturated” regions which do not grow any more, should it not be possible to achieve even greater reductions if we use explicit constraints? Such constraints would be made of other materials having negligible volumetric expansion compared to that of silicon.

The remainder of the paper is organized as follows. In Sec. 2, we recapitulate the formulation (in non-dimensional form) presented in our previous work [16] for a solid Si cylinder, adapting it for the general composite cylindrical geometry considered here; the specific geometries we study in the following sections may be obtained from this as special cases. Notably, we incorporate the two-way coupling between stress and concentration of Li in Si – meaning, that we account for the effects of both diffusion-induced stress and stress-enhanced diffusion. Our framework, while similar in spirit to a number of recent works [17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27] involving similar two-way coupling, is, however, based on the stress-dependent chemical potential presented in [28] for spherical particles. Additionally, we also for the possibility of plastic deformations in our formulation as an important departure from prior works on cylindrical Si particles. With this framework, we are able to study high charging rates which induce stresses strong enough to reach the yield stress of Si. We relegate the specification of the boundary conditions to each specific case described in the remaining sections. In Sec. 3, we discuss an annular cylinder of silicon with no explicit constraints examining two situations: lithiation occurs (i) only from the inner surface and (ii) only from the outer periphery. In Sec. 4, we discuss the case when a solid cylinder is constrained by an outer shell, examining the influence of the shell thickness and yield stress on axial growth. In Sec. 5, we discuss the case of an annular cylinder constrained from the inside by a concentric cylindrical rod, again examining the influence of the radius and yield stress (of the rod) on the axial growth. In Sec. 6, we present phase-diagrams of the percentage increases in length for both the outer and the inner constraint cases corresponding to different combinations of the radius of the constraint and the ratio of the yield-stress of the constraining material to that of silicon. In Sec. 7, we describe our first simplified model based on a simple force balance, and highlight two scaling relationships that characterize the phase-diagrams. In Sec. 8, we describe our second simplified model which can recover some of the qualitative trends of the phase-diagrams. Finally, in Sec. 9, we summarize our findings.

## 2 Mathematical Formulation

We consider a general circular cylinder made up of a number of concentric shells as shown in Fig. 2. In terms of a non-dimensional radial coordinate variable, , we have

The only geometric restriction we impose is that the volume of silicon in Region C is equal to that in a solid circular cylinder with radius , so that there is a consistent basis for comparison of the stresses and deformation induced by the lithation-driven expansion of silicon. Per unit length, this requirement translates as:

(1) |

We extend the formulation from our previous work on such a solid circular cylinder [16] to this situation. Any general point in this composite cylinder is, in cylindrical coordinates, , , , and the non-dimensionalized set of equations governing the lithiation-driven mechanical and chemical processes are:

(2) | ||||

(3) | ||||

(4) | ||||

(5) | ||||

(6) |

Eq. (2) governs the transport of Li in the material, with representing the non-dimensional measure of Li concentration, and the flux. Eq. (3) represents the mechanical equilibrium in the reference (undeformed) configuration in terms of the first Piola-Kirchhoff stresses, and , along the radial and circumferential (or, hoop) directions, respectively. Eqs. (4) and (5) govern the evolution of the plastic stretches, along the radial and hoop directions, respectively. The parameter is the characteristic strain rate, is the concentration-independent diffusivity of silicon, is a reference radius, and is the yield stress of silicon. Also, , , and are the deviatoric parts of the Cauchy stresses while is the effective stress; these stresses are defined later. Eq. (6) indicates that the plastic deformation is assumed to be volume-preserving (see [16] for details). The first Piola-Kirchhoff stresses may be expressed in terms of the strains, , , and , as

(7a) | ||||

(7b) | ||||

(7c) |

Here, and are the radial and the axial displacements respectively. Further, represents the volumetric expansion of Si or the constraining material due to lithiation, is the concentration-dependent non-dimensional modulus of elasticity, and is Poisson’s ratio. Note that the saturation level of Li in Si is set by the factor, . In the constraining material (of which Regions B and D are made), this saturation level is set at an arbitrarily small value of . It is this difference in the saturation levels which results in a lower volumetric expansion of Regions B and D compared to Region C, and imparts their functionality as external physical constraints against the expansion of silicon (since radial displacements at the interfaces of Regions B, C, and D must match). The parameter represents the coefficient of compositional expansion while represents the rate of change of modulus of elasticity with concentration; these values are taken to be the same for both silicon and the constraining material.

The stresses may be linked to the displacement by expressing the strains in terms of the displacement as

(8a) | ||||

(8b) | ||||

(8c) |

In Eq. (2) the non-dimensional flux is given by

(9) |

where

(10) |

with

(11a) | ||||

(11b) | ||||

(11c) |

Eqs. (4) and (5) are expressed in terms of the non-dimensionalized deviatoric parts of the Cauchy stress tensor given by

(12) |

where

(13a) | ||||

(13b) | ||||

(13c) |

Furthermore, the non-dimensional effective stress in Eqs. (4) and (5) is given by

(14) |

The boundary and the initial conditions associated with these equations depend on the particular configuration under consideration, and will be described on a case by case basis in the following sections. Common to all situations, however, is the physical condition that the cylinder ends are free. This implies that there can be no axial force acting on the ends, so that

(15) |

where and are the appropriate inner and outer limits of the radius.

Before proceeding, we note that in almost all the cases we investigate, the percentage increase in length is lower for higher values of when lithiation is from the outside (the exceptions are Figs. 4 (d) and 5 (d), which will be studied separately). The same trend is observed in the unconstrained solid cylinder case; see Fig. 1. This trend can be explained on the basis of the spatial heterogeneity of lithiation. For a relatively high influx rate, the amount of lithiation is higher in the periphery than in the inner region leading to a greater spatial heterogeneity of lithiated silicon over the cylinder radius compared to situations with a relatively low influx rate. The unlithiated region then constrains the lithiated periphery resulting in more plastic flow and lower axial expansion. The same explanation holds for annular silicon regions (irrespective of the nature of constraints) undergoing lithiation at different influx rates.

Material property or parameter | Value |
---|---|

, parameter used in activity constant | -29549 Jmol |

, parameter used in activity constant | -38618 Jmol |

, diffusivity of Si | 1 10 ms |

, characteristic strain rate for plastic flow in Si | 1 s |

, modulus of elasticity of pure Si | 90.13 GPa |

, stress exponent for plastic flow in Si | 4 |

, universal gas constant | 8.314 JKmol |

, initial radius of unlithiated Si electrode | 200 nm |

, temperature | 300 K |

, molar volume of Si | 1.2052 10 mmol |

, maximum concentration of Li in Si | 4.4 |

, coefficient of diffusivity | 0.18 |

, coefficient of compositional expansion | 0.2356 |

, rate of change of modulus of elasticity with concentration | -0.1464 |

, Poisson’s ratio of Si | 0.28 |

, initial yield stress of Si | 0.12 GPa |

[28] | |

[29] | |

[18] | |

[21] |

In what follows we look at special cases, focussing on the percentage increase in length that results from lithiation up to a 50% state of charge (SOC) for various values of the lithation rates. The motivation for choosing this particular value of the SOC is three-fold. First, in an actual battery, the anode particles which are farthest from the separator (nearest the current collector) might not reach a 100% lithiated state. So, a 50% SOC is representative for the entire collection of particles. Second, as the SOC builds up, there is increased resistance to further charging, and so going up to a 100% SOC is more time-consuming for our numerical simulations. Third, if the axial growth corresponding to 50% SOC is practically untenable, it certainly is so for higher values of the SOC. Thus, although our findings are based on the 50% SOC, our prescriptions to overcome undesirable scenarios are sufficiently general to be usefully applied for any value of the SOC. Unless otherwise stated, we use the values of the material properties and the parameters (taken from [16]) in Table 1. The constraint (either the inner or the outer) and the silicon regions are distinguished primarily by a difference in the volumetric expansion, and, in particular situations, also by a difference in the yield stress value.

## 3 Annular cylinder

To realize the annular cylinder geometry, we let the Regions B and D vanish, and retain Region C. Thus, . We consider two sub-cases: (a) lithiation only from the inside, and (b) lithiation only from the outside. The initial and boundary conditions corresponding to Eq. (2) are

(16) |

(17) | |||

(18) |

For both (a) and (b), the boundary conditions for the stresses are

(19) |

Again for both (a) and (b), the initial conditions corresponding to Eqs. (4) and (5) are:

(20) |

The predicted percentage increase in length corresponding to the two cases are shown in Fig. 3 (a) and (b), where . We choose this particular value of so that the radius of the annulus hole is significantly smaller than the silicon thickness. We observe that the percentage increase in length for different influx rates in (a) are independent of the charge rate unlike in (b) where such values are lower for higher influx rates. The reason for this trend is that strong hoop stresses are generated at the inner boundary of the annulus during lithiation which results in a higher value of the diffusion coefficient. However, since in (a) the influx is from the inner boundary near which diffusion is strong, there is greater spatial homogeneity of lithiated silicon. In (b) even though diffusion is faster near the inner boundary, spatial heterogeneity results at the outer boundary resulting in a decrease in axial extension.

We make two important remarks. First, the percentage increase in length when is almost the same as that when . Since increasing the value of also implies increasing the value of to preserve the silicon content, there is no practical benefit of having an annulus with larger inner radius. Second, we have also examined the case when the radial displacement at the outer periphery is constrained to be zero so that all radial growth is accommodated by filling in the hole of the annulus. However, such a set-up leads to extremely high () increases in length, and is, thus, antithetical to the overall objective of this investigation.

## 4 Outer constraint

We let Regions A and B vanish so that , and we have again a solid circular cylinder of silicon; Region C is thus defined by . The outer constraint (Region D) is retained; . Lithiation is considered to occur from the outside through Region D. The initial and boundary conditions are

(21) | |||

(22) | |||

(23) |

We now study the influence of the thickness and the yield stress of the constraining material on the percentage increase in length. We choose two different thicknesses such that the outer radius becomes and . The choice of is motivated by the desire to study the case where the constraining material thickness is one order of magnitude smaller than the radius of the silicon region; on the other hand, is chosen arbitrarily to look at the influence of a considerably larger thickness. Fig. 4 shows the predicted percentage increase in length for in (a) and (c), and for in (b) and (d). In (a) and (b) the yield stress value in Region D is the same as in Region C while in (c) and (d) it is 10 times higher in Region D than that in Region C.

We observe that an increase in reduces the percentage increase in length. The reason for this trend appears to be that Region D undergoes a lower volumetric expansion compared to Region C. However, displacements must be continuous across the interface between Regions D and C. This condition, therefore, constrains the axial deformation of Region C. Such a constraining influence, of course, gets stronger as Region D gets thicker. Additionally, increasing the yield stress value of Region D constrains the overall axial deformation even further. Thus, as seen in (d), the percentage increase in length is significantly reduced when the yield stress of Region D is 10 times higher than that of Region C and the thickness of Region D is high. An interesting observation is that, unlike other cases, for (d), the percentage increase in length corresponding to is slightly smaller than that corresponding to for sufficiently high values of the SOC.

## 5 Inner constraint

For this case, we let Regions A and D vanish. Thus the inner constraint is now a solid circular cylinder defined in the region, . Region C retains its annular shape, and is defined in the region . Lithiation is considered to occur only from the outside and not through Region B. The initial and boundary conditions are

(24) | |||

(25) |

(26) |

(27) |

As in the case of the outer constraint, we study the influence of the thickness (radius) and the yield stress of the constraining material on the percentage increase in length. We choose two different radii of the constraining material: and . Fig. 5 shows the percentage increase in length for in (a) and (c), and for in (b) and (d). In (a) and (b), yield stress value in Region B is the same as that of Region C, while in (c) and (d) it is 10 times higher in Region B than in Region C.

We observe that an increase in decreases the percentage increase in length because Region B undergoes a lower volumetric expansion compared to Region C. The continuity of displacements across the interface of Regions B and C imposes a constraint on the extent to which the silicon annulus (Region C) can expand in the axial direction. The greater the radius of the Region B, the higher is the constraining influence. Most interestingly, the percentage increase in length shows a significant reduction when the yield stress in Region B is 10 times higher than that in Region C and the radius of Region B is also comparatively high (see panel (d)). To isolate and understand the influence of this increased yield stress, we compare (b) and (d) since the radius of Region B is the same in both. Considering a particular lithiation rate, , for instance, we observe that in panel (d), there is a distinct transition in the percentage increase of length at SOC mAhgand another one much further on, just before SOC mAhg. In contrast to this trend, in panel (b), there is only one transition in the percentage increase of length at SOC mAhg. We probe the stress fields and the plastic stretches corresponding to these values of the SOC. We find that for the case of panel (d), in the region before the first transition (at SOC mAhg), Region C starts to yield from the outside, and the yield region grows in size; Region B, however, does not yield at all. Beyond the transition point, the whole of Region C is in a state of yield, while Region B does not yield yet. Thus, after the first transition, axial deformations are delimited by the small elastic deformations of the inner constraining cylinder. At the same time, further volumetric expansion in the silicon of Region C is easily accommodated through plastic flow in the radial direction without the need to increase the length. This plastic flow significantly reduces the rate of length increases (shown by a drop in the gradient of the plot). For panel (b), however, the only transition occurring at SOC mAhg corresponds to a situation where both Regions B and C have yielded. Following this transition, further volumetric expansion due to increasing SOC values is easily accommodated by plastic flow in both Regions B and C which results in comparatively higher percentage increase in length. This influence of the higher yield stress is however hardly noticeable when the inner radius is very small as in panel (c) with .

## 6 Phase diagrams for the explicitly constrained cases

Taking a cue from Figs. 4 (d) and 5 (d) which show that higher values of the yield stress and greater geometrical extent of the constraint result in smaller length increase, we study the combined influence of such factors through phase diagrams showing the percentage increase in length for the outer constraint as well as for the inner constraint case corresponding to different combinations of the thickness of the constraining material and the ratio of the yield stress of the constraining material to that of silicon. We do not present similar results for the simple annular cylinder because in that case we cannot define a ratio of the yield stresses. We first present the outer constraint case, and then the inner constraint case.

We map out in Figs. 6 and 7 a phase diagram of such percentage increases in length corresponding to various combinations of or and the ratio of the yield stress of the constraining material to that of silicon. We observe that the percentage increase in length is large for low values of the yield stress ratio and small geometrical extent of the constraint. Additionally, a decrease in the area of constraining region must be accompanied by an increase in the yield stress ratio in order to maintain the same percentage increase in length. In the case of the inner constraint, this finding is an important one for practical battery design because decreasing the inner radius allows decreasing the outer radius of the silicon annulus which is desirable from the perspective of reducing the overall volume.

When plotted on a log-log scale (with the intercepts along removed), the various contours of Fig. 7 are found to collapse, at least up to a yield stress ratio of around , on to a straight line with slope approximately equal to . This is clearly seen in Fig. 8.

To interpret the phase diagram presented in Fig. 7, we consider some simplified models which we describe now.

## 7 First simplified model for the explicitly constrained cases

We posit that the condition of zero net force in the axial direction is achieved through a simple force balance over the entire cross-section of the “composite” cylinder (in the current or deformed configuration) in such a way that a net tensile force in one region is balanced exactly by a compressive force in the other region. In this simplified picture, we assume no spatial variation of the axial stress in each of the regions (silicon and constraint), and that (after yield) the axial stresses in each region are simply equal to the corresponding yield stresses. Further, the evolution of growth (and, hence, the stresses) is not tracked; rather a relation between the radius and the length change is obtained in the final deformed configuration with the assumption of negligible contributions from elastic deformation. The percentage increase in volume is taken to be 155% corresponding to 50% SOC.

### 7.1 Outer constraint

The fractional increase in volume of the silicon cylinder (Region C) is given by

(28) | ||||

(29) |

where , and where we have used . Region D is assumed to undergo negligible volume change; thus

(30) | ||||

(31) |

The axial force balance then gives

(32) |

where with and being the constant (magnitudes of the) stresses in Regions C and D, and where we have used both (29) and (31).

To specify and we assume that both Regions C and D have yielded plastically by the time 50% SOC is reached. For such a situation, both and may be expected to be at their respective yield stress values. Guided by Eq. (32), a log-log plot of the raw simulation data shows a collapse (barring two outliers) of the vs points on to a straight line (at least up to ) with a slope of as seen in Fig. 9. This is in good agreement with the scaling relationship from Eq. (32) of .

### 7.2 Inner constraint

The fractional increase in volume of the silicon annulus (Region C) is given by

(33) | ||||

(34) |

where we have used . With the assumption that the radius of the constraint remains practically unchanged (based on the fact that its volumetric expansion is taken to be negligibly small compared to the silicon), the force balance gives

(35) | ||||

(36) |

where and are the constant (magnitudes of the) axial stresses in Regions B and C respectively, and . To specify and we need to consider different regimes of mechanical behaviour of Regions B and C.

We expect that when the yield stress of Region B is not too high, then both Regions B and C have yielded plastically by the time 50% SOC is reached. Both and may then be expected to be at their respective yield stress values. Therefore, we take which approximates the left-most section of each contour in Fig. 7. For high enough yield stresses in Region B, it might not yield even though Region C does. For such a situation, ( being the modulus of elasticity of the constraining material which is considered to be the same as that of Si, being the fractional increase in length), and this is a good approximation for relatively low values of . For Region C in a state of yield, . Using this in Eq. 36 gives horizontal lines which approximate the right-most part of each contour in Fig. 7. However the model is not quantitatively accurate in predicting axial length increases because we have assumed that the axial stresses are equal to the yield stress, whereas the charging rates are such that in the simulations, the axial stresses in both Regions C and D increase much beyond the respective yield stress values. Through Eq. (36), this simple model has predicted the correct scaling ratio, , observed in the collapsed plots of Fig. 8.

## 8 Second simplified model

We develop a second simplified model to capture the evolution of the stresses and the plastic stretches for the case of the annular silicon cylinder with an outer or an inner constraint. Here, we model the annular cylinder-constraint system by an equivalent two cylinder model, as shown in Fig. 10. Thus, Regions B, C, and D are each modelled by solid circular cylinders. Guided by the phase-diagrams and the first simplified model, we consider two regimes. In the first one, we consider that both the cylinders are flowing plastically, whereas in the second one, we consider that while the cylinder made of silicon has yielded, the cylinder made of the constraining material has not (i.e. it stays elastic). Note that the first (plastic-plastic regime) is applicable to the inner constraint case as well as the outer constraint one. The second regime, however, is applicable only to the inner constraint case; in other words, the system with the outer constraint is assumed never to operate within the elastic-plastic regime. Further, we also assume that the lithium concentration in the silicon cylinder is not affected by the stresses, and is, therefore, known independent of any mechanical considerations. In what follows, we consider the two regimes (plastic-plastic and elastic-plastic) separately. Within the first regime, we discuss the outer and the inner constraint cases, and within the second regime, we discuss only the inner constraint case.

### 8.1 Plastic constraint-Plastic Si regime

This regime is applicable for both the outer and the inner constraint cases. We assume separable solutions in the form of

(37a) | ||||

(37b) |

Further assuming small strains, we have

(38a) | ||||

(38b) | ||||

(38c) |

Substituting Eq. (37) in Eq. (38a) and (b), we have

(39) | |||

(40) |

and from Eq. (38c), we have

(41) |

where we have used .

Again using the smallness of the strains, we have for the Piola-Kirchhoff stresses:

(42a) | ||||

(42b) | ||||

(42c) |

Similarly for the Cauchy stresses, we have:

(43a) | ||||

(43b) | ||||

(43c) |

Comparing Eq. (42) and Eq. (43), we have:

(44) |

The plastic stretch evolution equations in the radial and the hoop directions are, respectively,

(45a) | |||

(45b) |

where . Since and they have the same initial condition of , Eq. (45) implies that . That is, . Therefore, from Eq. (44), we get . Using this in the mechanical equilibrium equation, we have:

(46) |

which implies that is spatially invariant. If we use the traction-free boundary condition then we have . And, this, in turn, implies from Eq. (44) that . In this situation, we have the following:

(47) |

Further,

(48) | ||||

(49) |

Then the (only independent) platic stretch evolution equation becomes

(50) |

Since we are assuming that both regions have yielded, the Heaviside function may be dropped.

In the following, we use the subscript “con” to denote Region B or D. It is also important to note that in the constraint cylinder there is no lithium influx; therefore, in these regions we have . Thus for the two cylinder system we have the following equations:

Region B or Region D

(51) | ||||

(52) |

Region C

(53) | ||||

(54) |

In this system, we have six unknowns: , , , , , and , and four equations. We, therefore, need two more equations to close the system. The first of these equations is given by the condition that the percentage increase in length of both cylinders should be the same; thus:

(55) |

The final equation is given by the condition of force-balance in the axial direction of the two-cylinder system which models the fact that the net force in the axial direction in the original annular-constraint system is zero. Thus, in the reference configuration we have

(56) |

Here, is either or , and are the areas of the two cylinders in the reference configuration, and is the ratio of the two reference areas.

Comparing Eq. (52) and Eq. (54), we have

(57) |

Taking logarithms and differentiating with respect to time, we obtain

(58) |

If we let , then using Eq. (56) in Eq. (58) leads to

(59) |

where . The sign functions may be resolved from a physical understanding of the situation. The length increase of the Region C cylinder is due to lithiation-induced growth, and under the condition of zero net axial force this growth is constrained; therefore, the axial stresses are expected to be compressive, implying . The length increase in the Region D or Region B cylinder is to due to stretching so that Eq. (55) is satisfied; therefore, . Eq. (59) is a non-linear algebraic equation in , and is solved iteratively for each time step given the values of and . This value of is then used to advance the plastic stretch using Eq. (53), following which the increase in length, , of the Region C cylinder (and, hence, of the two-cylinder system) is found from Eq. (54).

Note from Eq. (59) that depends on and only through . Thus this model also predicts that contours of constant axial extension should satisfy the scaling law .

### 8.2 Elastic constraint-Plastic Si regime

We consider the Region B cylinder to be in the elastic regime, and the Region C cylinder to have yielded. Again, we start by assuming separable solutions of the form Eq. (37). But since in Region B, we cannot take to be vanishingly small to obtain relations akin to Eq. (38) because that would result in (note that in Region B). However, we observe that

(60) |

Then, the expressions of the Piola-Kirchhoff stresses in the radial and the hoop directions, respectively, become:

(61a) | |||

(61b) |

from which we infer . From the mechanical equilibrium equation, this equality, in turn, implies

(62) |

Again using the traction-free boundary condition, this gives us . Therefore, , too. Then, from Eq. (61), we have

(63) |

Using this in the expression for and simplifying, we obtain