Improved Chemical‑Potential Model for Fast Deformation in Silicon Nanoparticle Electrodes of Li‑Ion Batteries

Introduction

High‑capacity electrodes such as silicon, which undergo volumetric changes of up to 400 % during lithiation, are pivotal for portable electronics, electric vehicles, and grid‑scale storage [1–3]. However, the diffusion‑induced stress generated by homogeneous volumetric expansion can trigger brittle fracture, shortening cycle life and reducing capacity [4]. Composite electrodes exhibit complex morphologies that challenge purely theoretical descriptions; most models treat material properties as spatially varying while neglecting interfacial effects. Typical electrode geometries considered in theory are spherical, cylindrical, and planar, with one‑ or two‑dimensional models accordingly.

Recent investigations on silicon nanoparticle electrodes have highlighted the interplay between lithiation‑induced phase transformation, morphological evolution, and stress generation [5–7]. These studies rely on the thermodynamic equilibrium of solids, where atomic diffusion is influenced by stress, and vice versa. Larche and Cahn’s framework [8] introduced a stress‑dependent chemical potential assuming small, isotropic deformations, while Wu’s formulation [9] incorporated the Eshelby momentum tensor. Cui et al. [10] extended this to finite deformations, yet their derivation remains accurate only for low deformation rates relative to diffusion. Given silicon’s large volumetric expansion during rapid charging, this assumption can lead to significant errors.

In this work, we derive a chemical‑potential expression valid for fast deformation without the low‑rate assumption, independent of electrode morphology. By applying the new model to silicon nanoparticle electrodes under both potentiostatic and galvanostatic control, we quantify the impact on stress and lithium distribution. We also identify a critical radius where the two models converge, providing insight into the spatial limits of stress influence.

Methods

Mechanical Equations

Lithium insertion induces volumetric change. We describe deformation using both Lagrangian and Eulerian coordinates:
$$\mathbf{U}=\mathbf{x}-\mathbf{X}\quad (1)$$ where U is the displacement field. The deformation gradient tensor is
$$\mathbf{F}=\frac{\partial \mathbf{x}}{\partial \mathbf{X}}=\mathbf{I}+\mathrm{Grad}\,\mathbf{U}\quad (2)$$

For a spherical particle with coordinates (R,Θ,Φ) (Lagrangian) and (r,θ,φ) (Eulerian), the gradient tensor simplifies to
$$\mathbf{F}=\begin{bmatrix}{F}_R&0&0\0&{F}_{\Theta}&0\0&0&{F}_{\Phi}\end{bmatrix}=\begin{bmatrix}1+\partial u/\partial R&0&0\0&1+u/R&0\0&0&1+u/R\end{bmatrix}\quad (3)$$

Deformation is decomposed into lithiation‑induced shape change and reversible elastic deformation:
$$\mathbf{F}=\mathbf{F}^e\mathbf{F}^c\quad (4)$$ with F^c isotropic:
$$\mathbf{F}^c=(1+\Omega C)^{1/3}\mathbf{I}\quad (5)$$ where Ω is the partial molar volume. Consequently,
$$\mathbf{F}^e=(1+\Omega(R)C)^{-1/3}\begin{bmatrix}1+\partial u/\partial R&0&0\0&1+u/R&0\0&0&1+u/R\end{bmatrix}\quad (6)$$

The Green–Lagrange strain tensor E is
$$\mathbf{E}=\frac{1}{2}\left(\mathbf{F}^T\mathbf{F}-\mathbf{I}\right)\quad (7)$$ yielding radial and tangential components
$$E_R^e=\tfrac{1}{2}\Big[\tfrac{(1+\partial u/\partial R)^2}{(1+\Omega(R)C)^{2/3}}-1\Big]\quad (9)$$
$$E_{\Theta}^e=E_{\Phi}^e=\tfrac{1}{2}\Big[\tfrac{(1+u/R)^2}{(1+\Omega(R)C)^{2/3}}-1\Big]\quad (10)$$

Assuming linear elasticity, the strain energy density is
$$W=\tfrac{1}{2}\mathbf{E}^e:\mathbf{C}:\mathbf{E}^e=\det(\mathbf{F}^c)\frac{E_h}{2(1+\nu)}\Big[\tfrac{\nu}{1-2\nu}\big(tr\,\mathbf{E}^e\big)^2+tr\,\mathbf{E}^e\mathbf{E}^e\Big]\quad (12)$$ where E_h and ν are Young’s modulus and Poisson’s ratio. The first Piola–Kirchhoff stress follows
$$\mathbf{P}=\det(\mathbf{F}^c)\frac{E_h}{1+\nu}\Big[\tfrac{\nu}{1-2\nu}tr\,\mathbf{E}^e+\mathbf{E}^e\Big]\mathbf{F}^e(\mathbf{F}^c)^{-1}\quad (13)$$

Componentwise,
$$P_R=\left(1+\Omega C\right)^{1/3}\frac{E_h}{2(1+\nu)(1-2\nu)}\left(1+\tfrac{\partial u}{\partial R}\right)\Big[(1-\nu)E_R^e+2\nu E_{\Theta}^e\Big]\quad (14)$$
$$P_{\Theta}=P_{\Phi}=\left(1+\Omega C\right)^{1/3}\frac{E_h}{2(1+\nu)(1-2\nu)}\left(1+\tfrac{u}{R}\right)\Big[\nu E_R^e+E_{\Theta}^e\Big]\quad (15)$$

Equilibrium without body forces requires
$$\frac{\partial P_R}{\partial R}+\frac{2(P_R-P_{\Phi})}{R}=0\quad (16)$$ with boundary conditions
$$u(0,t)=0,\;P_R(R_0,t)=0\quad (17)$$

Mass Transport Equation

Lithium concentration C(X,t) and flux J(X,t) satisfy
$$\frac{\partial C}{\partial t}+\mathrm{Div}\,\mathbf{J}=0\quad (18)$$. Under spherical symmetry, this reduces to
$$\frac{\partial C(R,t)}{\partial t}+\frac{1}{R^2}\frac{\partial\big(R^2J(R,t)\big)}{\partial R}=0\quad (19)$$. The flux follows Fick’s law with chemical‑potential driving force:
$$J=-\frac{CD}{R_gTF_{11}F_{11}}\frac{\partial \mu}{\partial R}\quad (20)$$.

The chemical potential derives from total internal energy density
$$\Pi=\varphi(C)+W(C,\mathbf{E}^e)\quad (22)$$, giving
$$\mu=\frac{\partial\Pi}{\partial C}=\mu_0(C)+\tau(\mathbf{E}^e,C)\quad (23)$$ where τ is the stress‑dependent part. For a dilute solution,
$$\mu_0(C)=\mu_0+R_gT\ln(\gamma C)\quad (24)$$ with γ=1.

Traditionally, the Helmholtz‑free‑energy approach yields
$$\tau_H=\frac{\partial W}{\partial C}= -\det(\mathbf{F}^e)\sigma_m\Omega\quad (25)$$ leading to
$$\mu={\mu}_0+R_gT\ln C-\det(\mathbf{F}^e)\Omega\sigma_m\quad (26)$$. Here, σ_m is the Cauchy hydrostatic stress
$$\sigma_m=\tfrac{1}{3}tr\,(\det^{-1}\mathbf{F}\mathbf{PF}^T)\quad (27)$$. This expression (26) is the traditional chemical‑potential model.

Alternatively, treating the Gibbs free energy yields
$$\tau_G=\frac{\partial W}{\partial C}\quad (28)$$ and consequently
$$\mu={\mu}_0+R_gT\ln C-\frac{\partial W}{\partial C}\quad (29)$$, our developed model. The mass‑transport system couples (19), (20), (26), and (29). Boundary conditions for potentiostatic (C(R₀,t)=C_max) and galvanostatic (J(R₀,t)=j₀(1+u/R)^2) operations are
$$C(R_0,t)=C_{\max},\;J(R_0,t)=j_0\big(1+u/R\big)^2\quad (30-31)$$ with initial
$$C(R,0)=0\quad (32)$$.

Numerical Implementation

Analytical solutions are intractable; we solved the coupled PDEs using COMSOL Multiphysics, adopting material parameters listed in Table 1. Dimensionless variables were used for clarity. State of charge (SOC) is computed as
$$SOC=\frac{\int_0^{R_0}C(R,t)R^2dR}{\int_0^{R_0}C_{\max}R^2dR}\quad (33)$$. Fluxes derived from the two chemical‑potential models are
$$J_H=\frac{\partial\tau_H}{\partial R},\;J_G=\frac{\partial\tau_G}{\partial R}\quad (34)$$.

Results and Discussion

Figure 1 illustrates lithium concentration, radial stress, and hoop stress profiles for a spherical Si nanoparticle under galvanostatic cycling. The developed (solid lines) and traditional (triangles) models yield nearly identical concentration distributions across all SOCs. Stress profiles differ only slightly, with modest central elevation at 46.7 % and 65.5 % SOC.

Figure 2 presents analogous data for potentiostatic cycling. Here, divergence between the two models is more pronounced, particularly in central stresses at intermediate SOCs. Because potentiostatic operation imposes a constant surface concentration, the charge rate—and consequently the deformation rate—is higher than in galvanostatic mode. The overall deformation for a given SOC remains similar, underscoring that the key factor is the deformation rate rather than the absolute deformation.

Figure 3 (from Ref. [17]) demonstrates that a 620‑nm Si nanoparticle can fully expand within one minute at 2 V versus Li, highlighting the relevance of fast‑deformation modeling. In such regimes, the two chemical‑potential models predict markedly different stress states, though experimental stress data remain unavailable for quantitative validation.

Figure 4 examines the flux ratio J_H/J_G under galvanostatic cycling at various charge currents. The ratio remains close to unity across the particle, with a slight increase near the center and a decrease near the surface. The spread widens with higher current, indicating stronger sensitivity to deformation rate. All curves intersect near a critical radius—termed the chemical‑potential‑independent region (CIR)—where fluxes converge regardless of the chosen model. CIR shifts closer to the surface as the current increases.

Figures 5 and 6 compare Cauchy hydrostatic stress distributions for the two models. In the CIR, σ_m is nearly zero for most of the lithiation period, explaining why the models agree there. Outside CIR, σ_m diverges, amplifying differences in chemical potential and thus flux.

Conclusions

We have introduced a chemical‑potential expression valid for rapid deformation, eliminating the low‑rate assumption inherent in traditional models. Comparison with the traditional formulation reveals negligible differences under galvanostatic cycling but significant deviations under potentiostatic control, driven by deformation rate. A critical radius—the CIR—emerges where fluxes and stresses are virtually identical, moving closer to the surface at higher currents. These findings highlight the necessity of using the developed model for accurate prediction of stress evolution in silicon nanoparticle electrodes during fast charging.

Availability of Data and Materials

The datasets analyzed in this study are available from the corresponding author upon reasonable request.

Abbreviations

CIR

A region where the diffusion fluxes predicted by the two chemical‑potential expressions are essentially identical.