Nomenclature

c,   lithium molar concentration in the reference state;

c_max,   saturation concentration at the stoichiometric limit;

Ω,   partial molar volume;

h₁,   initial thickness of the active layer;

H₁,   thickness of the deformed active layer;

h_c,   initial thickness of the current collector;

H_c,   thickness of the deformed current collector;

λ_θ,   in-plane stretch ratio;

λ_r,   transverse stretch ratio;

ε^s,   chemical volumetric eigen-strain;

E,   Young's modulus;

v,   Poisson's ratio;

σ_θ,   in-plane Cauchy stress;

σ_r,   transverse Cauchy stress;

σ_eq,   von Mises stress;

σ_Y,   yield stress;

F,   Faraday's constant;

R,   gas constant;

T,   temperature;

i_n,   electrical current density in the reference state;

D,   diffusion coefficient;

L_θ,   in-plane logarithmic stress;

L_r,   transverse logarithmic stresses;

Q,   state of charge (SOC).

1 Introduction

The bilayer electrodes made by bonding an active layer to a current collector layer have been widely used in lithium ion batteries, e.g., coin cells. Bilayer electrodes with very thick substrate (current collector) have been experimentally applied to the in-situ measurement of lithiation induced stress with multi-beam optical sensor (MOS)^[1-2]. Recently, a bilayer electrode with comparable thicknesses was proposed to measure the stress (strain) evolution with lithiation^[3-4] and the evolution of both the material property and the stress^[5]. Large pure bending of bilayer silicon composite electrodes has been observed during lithiation/delithiation cycles^[5]. Lang et al.^[6] showed that the LiFePO₄/silicon system could drive a large load with a rapid response to the electrochemical process. The lithium diffusion and diffusion induced stress in the bilayer electrode upon electrochemical cycling has been investigated by modeling in some early research^[7], e.g., the research on the effects of composition-dependent elastic modulus on diffusion^[8] and the research on the effects of hydrostatic stress on diffusion^[9]. However, these research was only conducted under the frame of elastic theory with small deformation, and thus was incapable of accounting for the full coupling between lithium diffusion and finite elastoplastic deformation.

Coupled diffusion finite deformation has been widely modeled in a broader range of electrode structures. Haftbaradaran et al.^[10] and Soni et al.^[11] proposed the stress-assisted diffusion model including the finite deformation hydrostatics stress for spherical particles and a thin film bonded to rigid substrate. Zhao et al.^[12] proposed a large plastic deformation model for the spherical particles in a lithiated electrode, considering the finite elastoplastic stretch caused by the diffusion induced stress and the influence of hydrostatic stress on diffusion. Bower et al.^[13] and Bower and Guduru^[14] developed a fully coupled constitutive theory to describe the finite deformation, plasticity, lithium ion diffusion, and electrochemical reactions in the electrodes with finite element implementation, and considered the simple applications of thin film bonded to rigid substrate. Cui et al.^[15-16] proposed a finite deformation stress-dependent diffusion theory, considering the interaction between the bulk diffusion and the chemical reaction for the applications of spherical particles in electrodes. Jia and Li^[17] conducted the chemo-mechanical modeling of spherical silicon to explain the experimentally observed stress mitigation due to the elastic softening during two-step lithiation. Anand^[18] developed a thermodynamically-consistent constitutive theory coupling the Cahn-Hilliard-type diffusion with large elastic-plastic deformation and accounting for the swelling and possible phase segregation caused by diffusing species. Di Leo et al.^[19] investigated the role of plastic deformation on electrochemical performance, in which the substrate was much thicker than the active layer. With the large elastic-plastic deformation method, Peigney^[20] studied the steady-state behavior of the medium coupling plasticity with diffusion under cyclic chemo-mechanical loadings.

Focusing on the two-way electrochemical-mechanical coupling behaviors of bilayer electrodes, this work aims to: (ⅰ) develop a fully coupling model for the diffusion induced finite elastoplastic bending of bilayer electrodes; (ⅱ) identify the key parameters governing the diffusion and finite elastoplastic bending; (ⅲ) analyze the finite plasticity behaviors of the bilayer electrode induced by lithium diffusion; (ⅳ) examine the effects of finite plasticity relaxation on lithium diffusion.

2 Finite elastoplastic bending induced by lithiation 2.1 Model and basic equations

We consider a bilayer plate electrode in lithium ion batteries (see Fig. 1(a)), where an active layer with the thickness h₁ and a current collector with the thickness h_c are perfectly bonded. Lithium-ions are inserted into the bilayer electrode from the free side surface of the active layer, causing the active layer to expand. This expansion is restricted by the current collector immune to lithiation, resulting in the finite pure bending of the bilayer electrode into a partial spherical shell (see Fig. 1(b)). The deformation follows the kinematic relationships as follows:

(1a)

(1b)

where λ_θ and λ_φ are in-plane stretch ratios in the two orthogonal directions, and are equal to each other, and λ_r is the stretch ratio in the thickness direction. θ₀ denotes the rotation of the cross-section per unit initial length. It does not vary spatially, but changes with the charging time t. Z is the thickness coordinate of the initial state, originating from the interface. r is the radial coordinate of the bent bilayer, and r_c, r_i, and r_a are the radial coordinates of the current collector surface, interface, and active layer surface, respectively. Equations (1a) and (1b) imply the relationship between the in-plane and the transverse stretch ratios, i.e.,

(2)

The volume expansion ratio of a lithiated fresh active material element without stress, which is denoted by J^s, is generally a function of the lithium molar concentration c defined as the total number of the lithium ions per unit volume of the initial state. If the partial molar volume Ω does not vary with the lithium concentration, the volume ratio, by definition, is given by

(3a)

(3b)

where ε^s is the volumetric strain due to the lithiation under the stress-free condition, and can be considered as the chemical eigen-strain. The chemical stretch ratio due to lithiation for isotropic materials can be expressed as

(4)

Following the multiplicative decomposition method^{[12-13, 18, 21-22]}, each of the total stretch ratios can be decomposed into three parts, i.e., the elastic stretch ratio, the plastic stretch ratio, and the chemical stretch ratio λ^s, as follows:

(5a)

(5b)

It is also assumed that the plastic deformation does not lead to any change in the volume, i.e.,

(6)

With Eqs. (5a) and (6), the elastic stretch ratios and elastic volume ratio can be expressed in terms of the total stretch ratios, the plastic stretch ratio, and the chemical stretch ratio as follows:

(7a)

(7b)

(7c)

The elastic behavior is assumed to follow a linear relationship between the logarithmic elastic strains and the conjugated logarithmic stresses, leading to

(8a)

(8b)

where E and ν are Young's modulus and Poisson's ratio, respectively. σ_r and σ_θ denote the transverse and in-plane Cauchy stresses, respectively. Substituting Eqs. (7a), (7b), and (7c) into Eqs. (8a) and (8b) gives

(9a)

(9b)

This is another form of the elastic constitutive equation, in which the influence of the lithium concentration on the stress is explicitly seen. The last term represents the lithiation induced stress, which should be removed for the current collector.

The plastic deformation is modeled with the perfectly plastic theory without strain hardening. The von Mises yield condition σ_eq=σ_Y for the finite pure bending deformation can be written as follows:

(10)

where σ_Y is the yield stress. Substituting Eqs. (9a) and (9b) into Eq. (10) yields

(11)

where ζ=sgn(σ_r-σ_θ) is the sign function. Hence, the plastic stretch can be expressed in terms of the total stretch ratio and the chemical eigen-strain as follows:

(12)

This expression is used only when the yield condition, i.e., Eq. (10), is satisfied. Otherwise, the material in the plastic zones undergoes elastic unloading deformation, and the plastic stretch ratio keeps unchanged. The equation of equilibrium gives

(13)

In this work, we focus on the plasticity of the active material in electrodes. We assume that the current collector is elastic. Then, all equations derived above are applicable for the current collector by prescribing the plastic stretch λ_θ^p=λ_r^p=1 and the chemical eigen-strain ε^s=0.

2.2 Active layer

The bilayer electrode bends, and may undergo elastoplastic deformation during lithiation. The plastic yield of the active layer may initiate from either the free surface or the interface, and grows into the interior. Now, we use the basic equations given in the preceding section to derive the governing equations for the active layer.

In the plastic zone, the solution for the Cauchy stresses can be obtained by combining the equilibrium equation (13) along with the yield criterion given by Eq. (10). If the plastic yielding starts from the free surface, we have

(14a)

(14b)

where Z=Z/h₁, the superscript (1) with parenthesis indicates the active layer, and λ_θ⁽¹⁾ (1) denotes the in-plane stretch of the free surface (Z=1). Substituting Eqs. (12), (14a), and (14b) into Eqs. (9a) and (9b) leads to a single independent nonlinear algebra equation where the total stretch ratios in the plastic zone starting from the free surface must be satisfied, i.e.,

(15)

where σ_Y=σ_Y/E₁.

If the plastic yielding starts from the interface, the Cauchy stresses can be expressed as follows:

(16a)

(16b)

where σ_r⁽¹⁾ (0) is the interfacial normal stress (Z=0). The total stretch ratios in the plastic zone starting from the interface must satisfy the following nonlinear algebra equation:

(17a)

(17b)

In the elastic zone, the equilibrium equation in terms of the total stretch ratios at a given charging time can be obtained by substituting Eqs. (9a) and (9b) with λ_θ^p(1)=1 into Eq. (13), i.e.,

(18)

where θ₀=h₁θ₀. The coefficients k₁, k₂, and k₃ are functions of the total stretch ratios λ_θ⁽¹⁾ and λ_r⁽¹⁾ and the chemical eigen-strain ε^s (see Appendix A).

Equation (2) for the kinematic relationship between the in-plane stretch ratio and the transverse stretch ratio is valid for both the elastic zone and the plastic zone in the active layer, and is rewritten here as follows:

(19)

2.3 Current collector

We focus on the plasticity of the active material in electrodes. We assume that the current collector is always elastic. With Eqs. (9a) and (9b), for the current collector, we prescribe the plastic stretch λ_r^p=1 and the chemical eigen-strain ε^s=0. Then, we can derive the governing equation for the current collector at a given charging time as follows:

(20)

(21)

where h_c=h_c/h₁. The superscript (c) with parenthesis indicates the current collector. The coefficients k₄ and k₅ are functions of the total stretch ratios λ_θ^(c) and λ_r^(c) (see Appendix A).

2.4 Boundary conditions and interface continuity conditions

For the elastic current collector, the stress-free surface condition, i.e., σ_r^(c)(-h_c)=0 with the substitution of Eq. (9a), leads to

(22)

For the active material, the stress-free surface condition, i.e., σ_r⁽¹⁾ (1)=0, has different mathematical expressions for elastic loading and plastic loading. For the elastic loading case, the boundary condition is

(23a)

In the plastic loading case, σ_r⁽¹⁾ (1)=0 with the substitution of Eqs. (9b) and (12) results in

(23b)

The displacement and traction must be continuous across the interface between the current collector and the active layer. With Eq. (2), we can see that the interface displacement continuity suggests that the in-plane total stretch ratio must be continuous across the interface, i.e.,

(24)

The traction continuity of the interface, i.e., σ_r⁽¹⁾(0)=σ_r^(c) (0), with the substitution of Eq. (9a) leads to

(25a)

for the elastic interface, where R_E=E₁/E_c, and

(25b)

for the interface where the active material is undergoing plastic yield.

Remember that the kinematic equation (1a) for in-plane stretch must hold true for the whole bilayer including the free-surfaces. It has been enforced by its derivatives via Eqs. (19) and (21) within both layers excluding the free surfaces. It needs to be specified at the free-surfaces of the current collector and the active layer that

(26a)

(26b)

For the pure bending of bilayer without external mechanical loading, the global balance could be enforced with zero resultant moment acting on a cross-section, i.e.,

(27)

Substituting Eqs.(9a) and (9b) into Eq.(27)yields

(28)

It is worth mentioning that the global balance of the resultant force acting on a cross-section is satisfied automatically via an integration of Eq. (12). Therefore, it is not needed to be specified again.

3 Diffusion equation coupled with finite elastoplastic bending

The transportation of lithium-ions in electrodes is modeled as a diffusion process, driven by the chemical potential gradient. The Gibbs free energy per unit volume of the reference state, which is the sum of the chemical potential energy G_chem and the mechanical potential energy, can be written as follows:

(29)

where G₀ is the free energy at the reference state, G_elastic is the elastic energy, and the last two terms are the work done by the system to its surroundings. L_θ and L_r are the logarithmic stresses conjugated to the total logarithmic strains ln λ_θ and ln λ_r. G_chem represents the change in the chemical potential energy during the lithiation without mechanical stress, and reads as^[18]

(30)

where T is the temperature, R is the gas constant, c is the molar concentration of lithium, and c_max is the saturation concentration at the stoichiometric limit. With the elastic constitutive equation, the Gibbs free energy can be expressed as follows:

(31)

The chemical potential or partial molar Gibbs free energy reads as

(32a)

(32b)

(32c)

where μ₀ is the chemical potential at the reference state, and μ_chem and μ_mech are parts of the chemical potential attributed to the lithium concentration and mechanical stress, respectively. Since the diffusion is one-dimensional in the thickness direction, the lithium flux is

(33)

where M is the mobility of lithium, and D=MRT is the diffusion coefficient. p₁, p₂, and p₃ are functions of the total stretch ratios λ_θ⁽¹⁾ and λ_r⁽¹⁾, the plastic stretch ratio λ_θ^p(1), and the chemical eigen-strain ε^s (see Appendix A). Combining Eq. (33) with the conservation condition we can derive the diffusion equation and write it in terms of the chemical eigen-strain ε^s as follows:

(34)

where t=Dt/h₁² is the dimensionless charging time. p₄, p₅, ..., p₉ are functions of the total stretch ratios λ_θ⁽¹⁾ and λ_r⁽¹⁾, the plastic stretch ratio λ_θ^p(1), and the chemical eigen-strain ε^s (see Appendix A). It is worth mentioning that our constitutive assumption, Eq. (33), is given for the initial state, and it is interchangeable with models based on the current deformation state^[14].

The initial molar concentration of lithium is assumed to be zero in the electrode, i.e.,

(35)

We consider the galvanostatic charging operation, and assume that the electrical current density is uniform over the electrode surface since the electrode is very wide. Let i_n denote the electrical current per unit area of the initial surface of the active layer, which is the electrical current density in the reference state. The lithium flux must satisfy at the surface and J=0 at the interface, where F=96 485.3 C/mol is Faraday's constant. With Eq. (33), we obtain the boundary conditions for lithium diffusion as follows:

(36a)

(36b)

where is the electrochemical load factor^[23], and characterizes the magnitude of the diffusion induced load applied to the electrode to cause mechanical deformation.

In summary, we have derived five governing equations for five unknown functions, i.e., ε^s, λ_θ⁽¹⁾, λ_r⁽¹⁾, λ_θ ^(c), and λ_r^(c), including one differential equation (34) for lithium diffusion, two kinematic equations (19) and (21) for the pure bending of the bilayer, one differential equation of equilibrium for the elastic current collector (20), and one equation of equilibrium for the active layer taking the form of one of Eqs. (15) and (17a). The unknown constant θ₀ is determined by the balance of the resultant moment equation (15) or (17a). The boundary conditions consist of eight equations, including two equations (36a) and (36b) for diffusion, two equations for the stress-free surface given by Eq. (22) and one of Eqs. (23a) and (23b), two interfacial mechanical continuity conditions given by Eq. (24) and one of Eqs. (25a) and (25b), two equations (26a) and (26b) for the kinematic relation at the free-surfaces. The initial condition for diffusion is given by Eq. (35). They describe the behavior of the fully coupling problems of diffusion and the finite elastoplastic bending of the bilayer electrodes. It is found that eight dimensionless parameters have emerged in the governing equations and the initial and boundary conditions, including six dimensionless material parameters, i.e., the maximum volumetric strain Ωc_max, the modulus ratio E₁/E_c, Poisson's ratios ν₁ and ν_c, the coupling coefficient the normalized yield stress σ_Y/E₁, a single dimensionless geometrical parameter, i.e., the thickness ratio h_c/h₁, and a single load parameter, i.e., the electrochemical load factor The governing equations are highly nonlinear, and are solved by the finite difference method along with the Newton-Raphson iterative scheme.

4 Results and discussion

In this section, we will use the fully coupled diffusion-elastoplastic model to examine the evolutional chemo-mechanical behaviors of the bilayer electrodes under galvanostatic charging operation. The state of charge (SOC) is proportional to the changing time by for the galvanostatic charging condition. The values of the dimensionless parameters in the model are given in Table 1 for numerical calculation. It needs to be clarified that some parameters of the lithiated silicon are concentration dependent, e.g., the Young's modulus would decrease from approximate 160 GPa of the initial state to 40 GPa of the fully lithiated state^[24-25] and the yield strength would correspondingly decrease from approximate 2 GPa to 1 GPa^{[1, 24-25]}. For simplification and demonstration, the material properties are assumed to be constant. The dimensionless parameters are chosen according to silicon with Young's modulus 80 GPa and the yield strength 1.5 GPa, which correspond to the median values and have also been used in other relevant calculations^{[12, 26]}.

As the lithiation starts and goes on, the bilayer electrode undergoes elastic pure bending till a critical value of Q when the active layer starts to yield at either the interface or the free surface. For the bilayer electrode with the parameters given in Table 1, Fig. 2(a) shows that the plastic yield starts from the interface since the von Mises stress at the interface attains the yield stress earlier than the von Mises stress at the free-surface reaches the yield stress. The plasticity zone starting from the interface grows towards the free-surface, and extends over a partial thickness of the active layer as lithiation goes on (see Fig. 2(b)). When Q=6.1%, the plasticity zone covers the whole thickness of the active layer. The active layer enters the plastic flow state. Therefore, the whole charge duration can be divided, according to deformation, into three periods, i.e., the elastic period, the partially plastic period, and the fully plastic period (see Fig. 2(a)). The small tensile in-plane stress at the interface during the very beginning charge time is observed in Fig. 2(a). This is because that the neutral surface of the bent bilayer is located near the free surface at the initial charge stage. The detailed discussion will be given in Subsection 4.2 for the temporal-spatial evolution of the stresses.

4.1 Temporal-spatial evolution of lithium concentration and chemical potential

First, let us examine how the elastic bending affects the temporal-spatial evolution of lithium in the active layer. The normalized molar concentration of lithium, which can be evaluated from the obtained chemical eigen-strain by c/c_max=ε^s/(Ωc_max ), is used for illustration. Figure 3(a) shows the lithium concentration distribution across the electrode thicknesat three typical values of Q during the elastic period. The results predicted by the diffusion model without the consideration of the stress effect are also included in Fig. 3(a) (see the dashed lines) in order to demonstrate the effect of the mechanical stress on the diffusion. The comparison between the coupled and uncoupled models indicates that the mechanical stress makes the concentration distribution less non-uniform at Q=0.5% and Q=1.8%, but the converse is true at Q=3.2%. This suggests that the mechanical stress facilitates the lithium diffusion, and drives more lithium ions into the interior region of the electrodes at the early stage of elastic bending, whereas the stress makes the lithium diffusion more difficult and impedes the lithium movement into the material interior at the later stage of elastic bending. This phenomenon can be explained by the chemical potential distribution shown in Fig. 3(b).

Remember that the spatial gradient of the chemical potential is the driving force of diffusion of lithium ions which always move towards the location of the low chemical potential from the location of high chemical potential. The total chemical potential given by Eq. (32a) includes the chemical and mechanical parts associated with the concentration and the mechanical stress. The normalized chemical part associated with concentration, μ_chem=μ_chem/(E₁Ω), having the highest value at the free-surface (lithium inlet) and the lowest value at the interface, decreases monotonically with the increase in the distance from the lithium inlet for any given Q. The normalized mechanical part of the chemical potential μ_mech=μ_mech/(E₁Ω) varies along the thickness in the same trend of the chemical part when Q=0.5% but in the opposite trend of the chemical part when Q=3.2%. Hence, in the early stage of the elastic period, both the chemical and mechanical parts of the chemical potential drive lithium ions from the lithium inlet towards the interface, while in the later stage of the elastic period, in contrast, the mechanical part will counteract the effect of the chemical part and make the lithium ions move backwards the lithium inlet. Therefore, the mechanical stress facilitates lithium diffusion in the early stage of the elastic period, while impedes the lithium diffusion at the later stage of the elastic period. The mechanical part of the chemical potential gradient is proportional to the stress gradient which will be discussed in Subsection 4.2.

Now, let us examine how the plastic yielding affects the temporal-spatial evolution of lithium in the active layer. In Fig. 4(a), the lithium concentration distribution across the electrode thickness is illustrated for three typical values of Q during the partially plastic period. The most interesting feature observed is that the lithium concentration variation in the thickness direction is much smaller in the plasticity zone of the active layer than in the elasticity zone for any given Q. This suggests that the plastic yield facilitates lithium transport towards uniform distribution. In order to find the reason for this interesting feature, the mechanical part and chemical parts of the chemical potential are plotted in Figs. 4(b) and 4(c), respectively. In the elasticity zone, the mechanical part of the chemical potential has the lowest value at the lithium inlet (free-surface), and increases with the distance from the lithium inlet. In the plastic zone, the mechanical part of the chemical potential has the highest value at the elastic/plastic boundary, and decreases with the increase in the distance from the lithium inlet. The chemical part of the chemical potential decreases monotonically with the distance from the lithium inlet. Hence, the mechanical part of the chemical potential has an opposite gradient of the chemical part in the elasticity zone. However, the gradient trend of both parts of the chemical potential is the same in the plasticity zone. Therefore, the mechanical stress in the elasticity zone impedes the lithium diffusion, while the mechanical stress in the plasticity zone facilitates lithium moving forward from the lithium inlet.

Figure 5(a) shows the through-thickness distribution of lithium concentration for three typical values of Q when the plastic yield has extended over the whole thickness of the active layer. Compared with the through-thickness variation of lithium during the partially plastic period shown in Fig. 4(a), the through-thickness variation of lithium during the fully plastic period is much smaller. In other words, the plastic flow of the active layer makes the lithium distribution over the whole thickness quite uniform. This agrees with the conclusion obtained for the partially plastic period, i.e., the plastic yielding facilitates the lithium diffusion. More interestingly, the lithium concentration near the interface (downstream the lithium) is slightly higher than that near the free-surface (lithium inlet) for Q=9.8% and Q=23%. This is an unusual feature in the sense that the amount of lithium near the inlet is less than that at the downstream. This distribution is possible because the total chemical potential has the highest value at the lithium inlet and the lowest at the interface and approaches the steady-state as Q increases (see Fig. 5(b)). It is known that the chemical part of the chemical potential drives lithium to move always towards the location of low concentration from the location of high concentration because the gradient of the chemical part of the chemical potential is proportional to the concentration gradient. Therefore, the chemical part of the chemical potential cannot produce a distribution that the concentration at the inlet is lower than the concentration in the interior. Figure 5(b) also shows that the mechanical part of the chemical potential has the highest value at the inlet and the lowest value at the interface. Therefore, the mechanical part of the chemical potential constantly provides the driving force to pump lithium ions from the inlet into the material interior, leading to abnormal lithium distribution.

4.2 Temporal-spatial evolution of deformation and stress

Here, let us see how stress and elastoplastic deformation develop in the electrode as lithiation goes on. Figures 6(a)-6(d) illustrate the through-thickness distributions of the in-plane stretch ratio, the transverse stretch ratio, the in-plane stress, and the transverse stress for typical values of Q during the elastic period, respectively. The total in-plane stretch ratio λ_θ and the total transverse stretch ratio λ_r consist of the eigen-stretch ratio due to lithiation and the elastic stretch produced by the mechanical stress. The eigen-stretch ratio λ^s ascends along the Z-direction according to the lithium concentration profile shown in Fig. 3(a). Since the in-plane stress is compressive and the transverse stress is tensile for Q=0.5% and 3.2% (see Figs. 6(c) and 6(d)), the in-plane elastic stretch ratio makes the total in-plane stretch ratio λ_θ smaller than the eigen-stretch ratio and the total transverse stretch ratio λ_r while larger than the eigen-stretch ratio. Furthermore, Fig. 6(b) shows that the transverse stretch is dominated by the chemical eigen-stretch, ascends along the Z-direction, and thereby has the biggest value at the lithium inlet and the lowest value at the interface.

The temporal-spatial evolution of the in-plane stress (see Fig. 6(c)) shows the phenomenon given in Fig. 2(b) more clearly. As can be seen in the figure, the in-plane stress at the free surface is always compressive for any changing time whereas the in-plane stress at the interface exhibits a small tensile value at the very beginning stage and a compressive value in the following time. Within an early stage charging time, a region near the free surface within the active layer of the bilayer is compressed, and the remaining part, which is separated by the neutral plane from the compressive region, is in tension due to the zero total in-plane force. The neutral plane of the bilayer starts from the free surface and moves towards the interface and the current collector as the charging starts and goes on. Therefore, the in-plane stress at the interface of the bilayer is always tensile at the beginning, changes its sign when the neutral plane arrives at the interface, and then is compressive continuously for the following charge time. Nevertheless, in most charging time, the in-plane stress within the whole active layer is compressive.

It is found that the in-plane stress gradient in the active layer changes its sign as Q increases from 0.5% to 3.2% (see Fig. 6(c)). This is the competition result of the diffusion induced stress and bending stress. At the beginning of lithiation, the curvature is small, and the overall stress is mainly determined by the lithiation induced stress that is always compressive. According to the lithium concentration profile shown in Fig. 3(a), we know that the value of the compressive stress must decrease with the distance from the lithium inlet. As more and more lithium ions enter the active layer, the bending deformation, that produces tensile stress at the free surface, will make the total stress gradient change its sign, which will result in the change in the sign of the gradient of the mechanical part of the chemical potential. Therefore, we can conclude that the mechanical stress facilitates lithium diffusion in the early stage of the elastic period but impedes the lithium diffusion at the later stage of the elastic period, which has been pointed out in Subsection 4.1.

In the partially plastic period, the through-thickness distributions of the in-plane stretch ratio, the transverse stretch ratio, the in-plane stress, and the transverse stress are illustrated for typical values of Q in Figs. 7(a)-7(d), respectively. The total in-plane stretch ratio λ_θ and the total transverse stretch ratio λ_r include three parts, i.e., the chemical eigen-stretch, the elastic stretch, and the plastic stretch. The in-plane elastic contraction and the transverse elastic expansion observed in the elastic period exhibit again here. More importantly, Figs. 7(a) and 7(b) show that the large in-plane plastic contraction and the large transverse plastic expansion are induced by the plastic flow in the plasticity zone. Therefore, it is true again that the total in-plane stretch ratio is smaller than the chemical eigen-stretch ratio whereas the total transverse stretch ratio is larger than the chemical eigen-stretch ratio. Furthermore, the total transverse stretch ratio is much bigger than the total in-plane stretch ratio for any given Q.

The transverse stretch ratio of the active layer ascends along the Z-direction in the early stage of the partially plastic period, e.g., Q=4.6%, keeping the same trend as that in the elastic period. However, the transverse stretch ratio profile becomes concave upward in the later stage of the partially plastic period (Q=5.8%), and its value at the interface is larger than the value at the free surface. The first reason is that the large transverse plastic stretch descends along the Z-direction (see Fig. 7(b)). Second, the elastic stretch descends along the Z-direction since the transverse stress decreases with Z. Third, the eigen-stretch ratio profile in the partially plastic period has become less non-uniform compared with that in the elastic period, although the eigen-stretch ratio profile still ascends along the Z-direction.

Figures 7(c) and 7(d) show that, in the plasticity zone, the compressive in-plane stress increases with Z while the tensile transverse stress decreases with Z, making the von Mises stress constant. Although the transverse stress is very small compared with the in-plane stress, the in-plane and transverse stress gradients induced by this small transverse stress are not small, which are the root causes of the stress-assisted diffusion in the plastic zone shown and discussed in the previous subsection. Therefore, the transverse stress generated by finite elastoplastic deformation plays a crucial role in the coupled diffusion in the plastic region.

In the fully plastic period, the through-thickness distributions of the in-plane stretch ratio, the transverse stretch ratio, the in-plane stress, and the transverse stress are illustrated for typical values of Q in Figs. 8(a)-8(d). The discussion and conclusions with respect to the stress and stretch given for the partial plasticity zone are also valid for the full plasticity zone. Furthermore, the figures show that the chemical eigen-stretch ratio becomes very uniform in this charging period, and the plastic deformation plays an increasingly important role in the total deformation. Both the in-plane stress and the transverse stress decrease monotonically with Z in the whole thickness of the active material for any given Q, leading to the distribution of the mechanical part of chemical potential shown in Fig. 5(b). Therefore, the in-plane stress and transverse stress gradients in Figs. 8(c) and 8(d) generated by the full plastic flow are the root causes of the abnormal lithium concentration distribution shown in Fig. 5(a).

4.3 Temporal evolution of curvature, thickness and film stress

By using Eq. (1a), we can express the dimensionless curvature of the interface as follows: κ=h₁/r(0)=θ₀/λ_θ⁽⁰⁾. The thickness of the active layer normalized by its initial value is computed by

The film stress is defined as the in-plane force resultant acting on the active layer divided by the area of its cross-section, and is given by They are plotted in Figs. 9(a)-9(c) as a function of Q. The curvature, the normalized thickness, and the film stress are linear with Q approximately before the plastic yielding initiates. Then, the change rates of the curvature, the thickness, and the film stress with respect to Q decline, increase, and vanish, respectively. There is a good consistency in the trend between the theoretical prediction and the experiment result^[22] for the film stress. The chemical eigen-stretch ratio λ^s=(1+Ωc_max × Q)^1/3 and the chemical eigen-volume ratio J^s=(1+Ωc_max × Q), which stand for the stretch and volume expansion of the active material under the free-standing condition, are also plotted in Fig. 9(b) for comparison. The change in the total thickness of the active layer is caused by the chemical eigen-stretch and mechanical deformation. For the active layer bonded on a rigid substrate, the in-plane deformation is completely restricted, and the volume expansion of the active material should be completely transformed into thickening. Therefore, the chemical eigen-volume ratio can be alternatively considered as the upper bound of the normalized thickness of the bilayer. The chemical eigen-stretch ratio under the free-standing condition can be considered as the normalized thickness of the active layer on a very soft substrate. Therefore, it is the lower bound of the normalized thickness of the bilayer. Our prediction on the normalized thickness is in-between the two bounds. It is interesting that the normalized thickness is quite larger than the chemical eigen-stretch ratio in the plastic period, and the difference between them, i.e., the excess thickening, is getting bigger as Q increases. This suggests that the finite plastic deformation is one of the main causes of the large excess thickening observed in experiment^[27].

5 Conclusions

We have developed a fully coupling model for the lithium diffusion induced large elastoplastic bending of bilayer electrodes, which accounts for the stress-assisted diffusion and the diffusion induced finite elastoplastic deformation simultaneously. The governing equations along with initial and boundary conditions are derived and solved by the finite difference method along with the Newton-Raphson iterative scheme. It is found that eight dimensionless parameters determine the evolutional chemical-mechanical behaviors of the bilayer electrodes, i.e., the electrochemical load factor , the maximum volumetric strain Ωc_max, the coupling coefficient the normalized yield stress σ_Y/E₁, the modulus ratio E₁/E_c, the thickness ratio h_c/h₁, and Poisson's ratios ν₁ and ν_c.

In the process of lithiation, the bilayer electrode undergoes three deformation stages, i.e., the elastic period, the partially plastic period, and the fully plastic period, in which different coupling behaviors are seen. In the elastic bending period, the mechanical stress facilitates lithium diffusion in the early stage of the elastic period, while impedes the lithium diffusion at the later stage of the elastic period. In the partially plastic period, during which the plasticity zone starts from the interface and extends over a partial thickness of the active layer, the finite elastoplastic bending generates the tensile transverse stress in the active layer, elevates the chemical potential at the elastic/plastic boundary with respect to that at the interface, and thereby drives lithium towards the interface. In the fully plastic period, the plasticity makes the abnormal distribution possible that the lithium concentration at the interface is higher than the concentration at the lithium inlet. The thickness of the active layer is much larger than the eigen-stretch due to lithiation when the bilayer electrode is plastically yielded.

Appendix A

(A1)

(A2)

(A3)

(A4)

(A5)

(A6)

(A7)

(A8)

(A9)

(A10)

(A11)

(A12)

(A13)

(A14)