A diffusion–reaction–deformation coupling model for lithiation of silicon electrodes considering plastic flow at large deformation

Abstract

Silicon (Si) is extensively used as an electrode in lithium (Li) ion batteries for its superiority in the charging capacity. However, the lithiation of Si would induce large deformation and significant stress jump across a two-phase interface in Si electrodes and might eventually lead to structure failure of batteries. To improve the performance and life of Li ion batteries, it is of great importance to exactly model the lithiation process. In the past decade, although many models have been developed for the lithiation, most of them just relate the deformation with the diffusion, without considering the effect of electrochemical reactions, so that they cannot simulate the formation of the two-phase interface. In this paper, we aim to develop a fully coupling diffusion–reaction–deformation model for large deformation cases and apply it to simulate the lithiation process, where a reaction barrier effect is proposed to describe the formation of the two-phase interface due to fast reaction and the mechanical behaviors of Si electrodes such as large plastic flow during lithiation process are predicted.

Appendix

In this Appendix we present the numerical procedure by finite difference scheme for modeling the lithiation of Si electrodes. First, Eq. (67') can be rewritten as

$$\frac{{\partial \tilde{\phi }^{{{\text{Li}}}} }}{{\partial \tilde{t}}} = k_{1} (\tilde{\phi }^{{{\text{Li}}}} ,\tilde{\ell })\frac{{\partial^{2} \tilde{\phi }^{{{\text{Li}}}} }}{{\partial \tilde{Z}^{2} }} - k_{2} (\tilde{\phi }^{{{\text{Li}}}} ,\tilde{\ell })\frac{{\partial \tilde{\phi }^{{{\text{Li}}}} }}{{\partial \tilde{Z}}} + k_{3} (\tilde{\phi }^{{{\text{Li}}}} ,\tilde{\ell })$$

(A1)

where

$$\begin{aligned} & k_{1} = \frac{1}{{\lambda^{2} }} \\ & k_{2} = \frac{3}{{\lambda^{3} }}\frac{\partial \lambda }{{\partial \tilde{Z}}} - - \frac{1}{{\lambda^{2} }}\frac{{\partial \tilde{\tau }^{{{\text{Li}}}} }}{{\partial \tilde{Z}}} \\ & k_{3} = - \left( {\frac{x}{{\lambda^{2} }}\frac{{\partial^{2} \tilde{\ell }}}{{\partial \tilde{Z}^{2} }} + \frac{{\tilde{C}^{{{\text{Li}}}} }}{{\lambda^{3} }}\frac{{\partial^{2} \lambda }}{{\partial \tilde{Z}^{2} }} - \frac{{\tilde{C}^{{{\text{Li}}}} }}{{\lambda^{2} }}\frac{{\partial^{2} \tilde{\tau }^{{{\text{Li}}}} }}{{\partial \tilde{Z}^{2} }}} \right) \\ & \quad + \left( {\frac{3x}{{\lambda^{3} }}\frac{\partial \lambda }{{\partial \tilde{Z}}}\frac{{\partial \tilde{\ell }}}{{\partial \tilde{Z}}} - \frac{x}{{\lambda^{2} }}\frac{{\partial \tilde{\tau }^{{{\text{Li}}}} }}{{\partial \tilde{Z}}}\frac{{\partial \tilde{\ell }}}{{\partial \tilde{Z}}} - \frac{{2\tilde{C}^{{{\text{Li}}}} }}{{\lambda^{3} }}\frac{\partial \lambda }{{\partial \tilde{Z}}}\frac{{\partial \tilde{\tau }^{{{\text{Li}}}} }}{{\partial \tilde{Z}}} + \frac{{3\tilde{C}^{{{\text{Li}}}} }}{{\lambda^{4} }}\left( {\frac{\partial \lambda }{{\partial \tilde{Z}}}} \right)^{2} } \right) \\ \end{aligned}$$

(A2)

We solve Eq. (A1) numerically using a forward in time, centered space finite difference scheme. A one-dimensional grid is constructed of $m + 1$ points with uniform spacing $\Delta \tilde{Z}$. A superscript $n$ and a subscript $i$ are used to denote the time increment and the spatial location, respectively, where $i = 1$ corresponds to $\tilde{Z} = 0$, and $i = m + 1$ corresponds to $\tilde{Z} = 1$. Hence, we approximate $\tilde{\phi }^{{{\text{Li}}}}$ with the finite difference equations [47]

$$\begin{aligned} & \frac{{\partial \tilde{\phi }^{{{\text{Li}}}} }}{{\partial \tilde{t}}} \approx \frac{{\tilde{\phi }_{i}^{n + 1} - \tilde{\phi }_{i}^{n} }}{{\Delta \tilde{t}}} \\ & \frac{{\partial \tilde{\phi }^{{{\text{Li}}}} }}{{\partial \tilde{Z}}} \approx \frac{{\tilde{\phi }_{i + 1}^{n} - \tilde{\phi }_{i - 1}^{n} }}{{2\left( {\Delta \tilde{Z}} \right)}} \\ & \frac{{\partial^{2} \tilde{\phi }^{{{\text{Li}}}} }}{{\partial \tilde{Z}^{2} }} \approx \frac{{\tilde{\phi }_{i + 1}^{n} - 2\tilde{\phi }_{i}^{n} + \tilde{\phi }_{i - 1}^{n} }}{{\left( {\Delta \tilde{Z}} \right)^{2} }} \\ \end{aligned}$$

(A3)

Here, for simplicity, we omit a superscript Li of the discrete value $\tilde{\phi }_{i}^{n + 1}$ and use $\Delta \tilde{t}$ to denote the time increment. In this way, the discrete form of Eq. (A1) is written as

$$\dot{\tilde{\phi }}_{i} + \left( { - \frac{{k_{1} }}{{\left( {\Delta \tilde{Z}} \right)^{2} }} - \frac{{k_{2} }}{{2\Delta \tilde{Z}}}} \right)\tilde{\phi }_{{i - {1}}} + \frac{{2k_{1} }}{{\left( {\Delta \tilde{Z}} \right)^{2} }}\tilde{\phi }_{i} + \left( { - \frac{{k_{1} }}{{\left( {\Delta \tilde{Z}} \right)^{2} }} + \frac{{k_{2} }}{{2\Delta \tilde{Z}}}} \right)\tilde{\phi }_{{i + {1}}} = k_{3}$$

(A4)

where a super dot represents the time derivative and the superscript $n$ is hidden because the time is not discrete at this moment.

Applying a forward and backward space finite difference to the boundary conditions at $\tilde{Z} = 0$ and $\tilde{Z} = 1$, respectively, Eq. (68′) yields, respectively,

$$\begin{aligned} & \tilde{\phi }_{2} = \tilde{\phi }_{1} \\ & \tilde{\phi }_{m + 1} = {\tilde{{J}}}_{0} \Delta \tilde{Z}\lambda^{2} \tilde{\phi }_{{{\text{max}}}} + \left( {1 - {\tilde{{J}}}_{0} \Delta \tilde{Z}\lambda^{2} } \right)\tilde{\phi }_{m} \\ \end{aligned}$$

(A5)

with the normalized flux rate $\tilde{J}_{0}$ being defined by $\tilde{J}_{0} = \frac{{J_{0} h_{0} }}{D}$ and the assumptions $\left. {\frac{\partial \ell }{{\partial \tilde{Z}}}} \right|_{{\tilde{Z}{ = 0}}} = 0$ and $\left. {\frac{\partial \ell }{{\partial \tilde{Z}}}} \right|_{{\tilde{Z}{ = 1}}} = 0$.

Next, we write Eq. (A4) in a vector form as

$$\left[ {0\;\;1\;\;0} \right]\left[ \begin{gathered} \dot{\tilde{\phi }}_{i - 1} \hfill \\ \dot{\tilde{\phi }}_{i} \hfill \\ \dot{\tilde{\phi }}_{i + 1} \hfill \\ \end{gathered} \right] + \left[ {\frac{{ - k_{1} }}{{\left( {\Delta \tilde{Z}} \right)^{2} }} - \frac{{k_{2} }}{{2\Delta \tilde{Z}}}\;\;\;\frac{{2k_{1} }}{{\left( {\Delta \tilde{Z}} \right)^{2} }}\;\;\;\frac{{ - k_{1} }}{{\left( {\Delta \tilde{Z}} \right)^{2} }} + \frac{{k_{2} }}{{2\Delta \tilde{Z}}}} \right]\left[ \begin{gathered} \tilde{\phi }_{i - 1} \hfill \\ \tilde{\phi }_{i} \hfill \\ \tilde{\phi }_{i + 1} \hfill \\ \end{gathered} \right] = k_{3}$$

(A6)

Applying Eq. (A5) in Eq. (A6), we obtain

$$\begin{aligned} & \left[ {\begin{array}{*{20}c} 1 & 0 \\ \end{array} } \right]\left[ \begin{gathered} \dot{\tilde{\phi }}_{2} \hfill \\ \dot{\tilde{\phi }}_{3} \hfill \\ \end{gathered} \right] + \left[ {\frac{{k_{1} }}{{\left( {\Delta \tilde{Z}} \right)^{2} }} - \frac{{k_{2} }}{{2\Delta \tilde{Z}}}\;\;\;\frac{{ - k_{1} }}{{\left( {\Delta \tilde{Z}} \right)^{2} }} + \frac{{k_{2} }}{{2\Delta \tilde{Z}}}} \right]\left[ \begin{gathered} \tilde{\phi }_{2} \hfill \\ \tilde{\phi }_{3} \hfill \\ \end{gathered} \right] = k_{3} \\ & \left[ {\begin{array}{*{20}c} 0 & 1 \\ \end{array} } \right]\left[ \begin{gathered} \dot{\tilde{\phi }}_{m - 1} \hfill \\ \dot{\tilde{\phi }}_{m} \hfill \\ \end{gathered} \right] + \left[ {\frac{{ - k_{1} }}{{\left( {\Delta \tilde{Z}} \right)^{2} }} - \frac{{k_{2} }}{{2\Delta \tilde{Z}}}\;\;\;\frac{{2k_{1} }}{{\left( {\Delta \tilde{Z}} \right)^{2} }} + \left( {\frac{{ - k_{1} }}{{\left( {\Delta \tilde{Z}} \right)^{2} }} + \frac{{k_{2} }}{{2\Delta \tilde{Z}}}} \right)\left( {1 - \tilde{J}_{0} \Delta \tilde{Z}\lambda^{2} } \right)} \right]\left[ \begin{gathered} \tilde{\phi }_{m - 1} \hfill \\ \tilde{\phi }_{m} \hfill \\ \end{gathered} \right] \\ & \quad = k_{3} - \left( {\frac{{ - k_{1} }}{{\left( {\Delta \tilde{Z}} \right)^{2} }} + \frac{{k_{2} }}{{2\Delta \tilde{Z}}}} \right)\tilde{J}_{0} \Delta \tilde{Z}\lambda^{2} \tilde{\phi }_{{{\text{max}}}} \\ \end{aligned}$$

(A7)

Then, assembling Eqs. (A6) and (A7), we get

$${\mathbf{D}}^{n} \dot{\tilde{\phi }}^{n} + {\mathbf{K}}^{n} \tilde{\phi }^{n} = {\mathbf{P}}^{n}$$

(A8)

with the matrixes ${\mathbf{D}}^{n}$, ${\mathbf{K}}^{n}$, ${\mathbf{P}}^{n}$, and the vectors $\dot{\tilde{\phi }}^{n}$, $\tilde{\phi }^{n}$ being

$$\begin{aligned} & {\mathbf{D}}^{n} = \left( {\begin{array}{*{20}c} 1 & \ldots & 0 \\ \vdots & \ddots & \vdots \\ 0 & \cdots & 1 \\ \end{array} } \right)_{m - 1} ,\;\;\;\;\dot{\tilde{\phi }}^{n} = \left[ \begin{gathered} \dot{\tilde{\phi }}_{2}^{n} \hfill \\ \; \vdots \hfill \\ \dot{\tilde{\phi }}_{m}^{n} \hfill \\ \end{gathered} \right]_{m - 1} ,\;\;\;\;\tilde{\phi }^{n} = \left[ \begin{gathered} \tilde{\phi }_{2}^{n} \hfill \\ \; \vdots \hfill \\ \tilde{\phi }_{m}^{n} \hfill \\ \end{gathered} \right]_{m - 1} , \\ & {\mathbf{P}}^{n} = \left[ \begin{gathered} \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;k_{3}^{n} \hfill \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;k_{3}^{n} \hfill \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; \vdots \hfill \\ k_{3}^{n} - \left( {\frac{{ - k_{1}^{n} }}{{\left( {\Delta \tilde{Z}} \right)^{2} }} + \frac{{k_{2}^{n} }}{{2\Delta \tilde{Z}}}} \right)\tilde{J}_{0} \Delta \tilde{Z}\lambda^{2} \tilde{\phi }_{{{\text{max}}}} \hfill \\ \end{gathered} \right]_{m - 1} , \\ & {\mathbf{K}}^{n} = {\text{diag}}\left( {{\mathbf{v}}_{1}^{n} } \right) + {\text{diag}}\left( {{\mathbf{v}}_{2}^{n} ,1} \right) + {\text{diag}}\left( {{\mathbf{v}}_{3}^{n} , - 1} \right) \\ \end{aligned}$$

(A9)

where the main diagonal ${\mathbf{v}}_{1}^{n}$ has components: ${\mathbf{v}}_{1}^{n} \left( 1 \right) = \frac{{k_{1}^{n} }}{{\left( {\Delta \tilde{Z}} \right)^{2} }} - \frac{{k_{2}^{n} }}{{2\Delta \tilde{Z}}},$${\mathbf{v}}_{1}^{n} \left( {2:m - 2} \right) = \frac{{2k_{1}^{n} }}{{\left( {\Delta \tilde{Z}} \right)^{2} }},$${\mathbf{v}}_{1}^{n} (m - 1) = \frac{{2k_{1}^{n} }}{{\left( {\Delta \tilde{Z}} \right)^{2} }} + \left( {\frac{{ - k_{1}^{n} }}{{\left( {\Delta \tilde{Z}} \right)^{2} }} + \frac{{k_{2}^{n} }}{{2\Delta \tilde{Z}}}} \right)\left( {1 - {\tilde{{J}}}_{0} \Delta \tilde{Z}\lambda^{2} } \right)$, the first upper diagonal ${\mathbf{v}}_{2}^{n}$ has components: ${\mathbf{v}}_{2} (1:m - 2) = \frac{{ - k_{1}^{n} }}{{\left( {\Delta \tilde{Z}} \right)^{2} }} + \frac{{k_{2}^{n} }}{{2\Delta \tilde{Z}}}$, and the first lower diagonal ${\mathbf{v}}_{3}^{n}$ has components: ${\mathbf{v}}_{3}^{n} \left( {1:m - 2} \right) = \frac{{ - k_{1}^{n} }}{{\left( {\Delta \tilde{Z}} \right)^{2} }} - \frac{{k_{2}^{n} }}{{2\Delta \tilde{Z}}}$.

Finally, using the forward finite difference in time for Eq. (A3)₁ in Eq. (A8), values of $\tilde{\phi }^{n}$ at the grid points are updated in each time step by

$$\tilde{\phi }^{n + 1} = (\frac{{{\mathbf{D}}^{n} }}{{\Delta \tilde{t}}})^{ - } \left( {\left( {\frac{{{\mathbf{D}}^{n} }}{{\Delta \tilde{t}}} - {\mathbf{K}}^{n} } \right)\tilde{\phi }^{n} + {\mathbf{P}}^{n} } \right)$$

(A10)

with the initial condition $\tilde{\phi }^{1} = 0$.

As for chemical kinetics Eq. (76'), we update them, respectively, by

$$\tilde{\ell }^{n + 1} = \tilde{\kappa }_{0}^{ + } \lambda^{ - 2} \exp \left( {0.01\tilde{\tau }} \right) \, \left( {\tilde{\phi }^{n} - x\tilde{\ell }^{n} } \right)^{2} \left( {1 - \tilde{\ell }^{n} } \right)\Delta \tilde{t} + \tilde{\ell }^{n}$$

(A11)

with the initial condition $\tilde{\ell }^{1} = 0$ and the parameters $\tilde{\kappa }_{0}^{ + }$ defined by $\tilde{\kappa }_{0}^{ + } = \kappa_{0}^{ + } \frac{{h_{0}^{2} }}{D}C_{0}^{{{\text{Si}}}}$.

Equation (62) is discretized to calculate the plastic stretch by the following implicit schemes:

$$\lambda^{{{\text{p}},\;n}} = \lambda^{{{\text{p}},\;n - 1}} \exp \left( { - \dot{d}_{0} \frac{{\sigma^{,\;n} }}{{\left| {\sigma^{,\;n} } \right|}}\left( {\frac{{\left| {\sigma^{,\;n} } \right|}}{{\sigma_{0} }} - 1} \right)^{{m_{0} }} H\left( {\frac{{\left| {\sigma^{,\;n} } \right|}}{{\sigma_{0} }} - 1} \right)\Delta t} \right),\;\;{\text{when}}\;\left| {\sigma^{,\;n} } \right| \le \sigma_{0}$$

(A12)

Here $n$ and $n - 1$ after a comma indicate a quantity belongs to the current and previous step.

A summary of the numerical procedure for lithiation in silicon electrode is provided in Box 1. To ensure the convergence, the reasonable time step and spatial step should be carefully chosen. For this case, $\Delta \tilde{t} = 10^{ - 6}$ and $\Delta \tilde{Z} = 0.01$ are taken for calculation.

Box 1 Summary of algorithm for lithiation of Si