Simulation of the Deformation for Cycling Chemo-Mechanically Coupled Battery Active Particles with Mechanical Constraints

TL;DR

The paper develops a thermodynamically consistent chemo-mechanical model for lithiation/delithiation of silicon active particles under constrained swelling, using a finite-deformation framework with a multiplicative decomposition $\mathbf{F} = \mathbf{F}_{\text{ch}} \mathbf{F}_{\text{el}}$ and isotropic chemical expansion $\mathbf{F}_{\text{ch}} = \lambda_{\text{ch}} \mathbf{I}$ with $\lambda_{\text{ch}} = \sqrt[3]{1 + v_{\text{pmv}} c}$. It couples chemical and elastic effects through a free-energy density $\psi(c, \boldsymbol{\nabla}_0 \boldsymbol{u}) = \psi_{\text{ch}}(c) + \psi_{\text{el}}(c, \boldsymbol{\nabla}_0 \boldsymbol{u})$, where $\psi_{\text{ch}}$ is tied to the OCV curve and $\psi_{\text{el}}$ uses a linear elastic $(\mathds{C}, \mathbf{E}_{\text{el}})$ response; diffusion is governed by $\partial_t c = -\boldsymbol{\nabla}_0 \cdot \mathbf{N}$ with $\mathbf{N} = -m(c, \boldsymbol{\nabla}_0 u) \boldsymbol{\nabla}_0 \mu$ and $\mu = \partial_c \psi$. Obstacle contact is modeled by Signorini-type boundary conditions, leading to a variational inequality that is solved via a primal-dual active-set method interpreted as a semismooth Newton iteration, integrated with space-time adaptive finite elements. The numerical experiments on 1D and 2D geometries with silicon demonstrate substantial stress amplification and the emergence of a lithium-poor region near the obstacle in 2D, while the adaptive solver distributes computational effort efficiently across evolving contact zones and cycling. The work provides a scalable framework for studying mechanical degradation and aging in constrained battery particles, with potential extensions to more complex 3D geometries and fracture mechanisms.

Abstract

Next-generation lithium-ion batteries with silicon anodes have positive characteristics due to higher energy densities compared to state-of-the-art graphite anodes. However, the large volume expansion of silicon anodes can cause high mechanical stresses, especially if the battery active particle cannot expand freely. In this article, a thermodynamically consistent continuum model for coupling chemical and mechanical effects of electrode particles is extended by a change in the boundary condition for the displacement via a variational inequality. This switch represents a limited enlargement of the particle swelling or shrinking due to lithium intercalation or deintercalation in the host material, respectively. For inequality constraints as boundary condition a smaller time step size is need as well as a locally finer mesh. The combination of a primal-dual active set algorithm, interpreted as semismooth Newton method, and a spatial and temporal adaptive algorithm allows the efficient numerical investigation based on a finite element method. Using the example of silicon, the chemical and mechanical behavior of one- and two-dimensional representative geometries for a charge-discharge cycle is investigated. Furthermore, the efficiency of the adaptive algorithm is demonstrated. It turns out that the size of the gap has a significant influence on the maximal stress value and the slope of the increase. Especially in two dimension, the obstacle can cause an additional region with a lithium-poor phase.

Figures (13)

• Figure 1: The total deformation $\mathbf{F}{}$ can be multiplicatively decomposed into an elastic part $\mathbf{F}{}_\text{el}$ and a chemical part $\mathbf{F}{}_\text{ch}$, compare schoof2022parallelization.
• Figure 2: Schematic sketch of lithiation and following delithiation of a representative battery active particle with volume change, getting in contact with the obstacle and detaching from the obstacle again.
• Figure 3: All boundary parts for the example of a two-dimensional quarter disk domain $\Omega_0$ split up into two artificial boundaries $\Gamma_{0,y}$ and $\Gamma_{0,x}$ with additional Dirichlet constraints for the displacement $\boldsymbol{u}{}$ and the potential contact boundary $\Gamma_\mathcal{P}$ subdivided into the active contact boundary $\Gamma_\mathcal{A}$ and the inactive boundary $\Gamma_\mathcal{I}$.
• Figure 4: Dimension reduction of a three-dimensional unit sphere with surrounded obstacle to the one-dimensional unit interval with spherical symmetry and the gap function $g$, based on castelli2021numerical.
• Figure 5: Dimension reduction of a three-dimensional nanowire with surrounded rectangular obstacle to the two-dimensional quarter disk and the (time-independent) quadratic obstacle with the gap function $(g_x, g_y)^\textsf{T} = \boldsymbol{g}{} = \hat{\boldsymbol{g}{}} - \boldsymbol{X}{}_0 = (\hat{g}_x, \hat{g}_y)^\textsf{T} - (X_{0,x}, X_{0,y})^\textsf{T}$ with the components $g_x$ and $g_y$ in $x$- and $y$-directions, respectively.