Modèle de potentiel chimique amélioré pour la déformation rapide dans les électrodes nanoparticules de silicium des batteries Li-ion

Résumé

Nous présentons une expression de potentiel chimique raffinée qui élimine l'hypothèse de faible taux de déformation, permettant une modélisation précise de la contrainte induite par diffusion et de la distribution de lithium dans les électrodes nanoparticules de silicium. En comparant la nouvelle expression avec la traditionnelle sous cyclage potentiostatique et galvanostatique, nous démontrons une divergence prononcée pendant l'opération potentiostatique, tandis que les résultats galvanostatiques restent largement inchangés. Un rayon critique — où le flux de Li et la contrainte hydrostatique sont essentiellement identiques pour les deux modèles — est identifié pour la première fois.

Introduction

Les électrodes à haute capacité comme le silicium, qui subissent des changements volumétriques jusqu'à 400 % lors de la lithiation, sont cruciales pour l'électronique portable, les véhicules électriques et le stockage à l'échelle du réseau [1-3]. Cependant, la contrainte induite par diffusion générée par l'expansion volumétrique homogène peut déclencher une fracture fragile, raccourcissant la durée de vie et réduisant la capacité [4]. Les électrodes composites présentent des morphologies complexes qui défient les descriptions purement théoriques ; la plupart des modèles traitent les propriétés matérielles comme spatialement variables tout en négligeant les effets d'interface. Les géométries d'électrode typiques considérées en théorie sont sphérique, cylindrique et plane, avec des modèles unidimensionnels ou bidimensionnels en conséquence.

Des investigations récentes sur les électrodes nanoparticules de silicium ont mis en lumière l'interaction entre la transformation de phase induite par lithiation, l'évolution morphologique et la génération de contrainte [5-7]. Ces études reposent sur l'équilibre thermodynamique des solides, où la diffusion atomique est influencée par la contrainte, et vice versa. Le cadre de Larche et Cahn [8] a introduit un potentiel chimique dépendant de la contrainte en assumant de petites déformations isotropes, tandis que la formulation de Wu [9] a incorporé le tenseur de momentum d'Eshelby. Cui et al. [10] ont étendu cela aux déformations finies, pourtant leur dérivation reste précise seulement pour de faibles taux de déformation relatifs à la diffusion. Donnée l'expansion volumétrique importante du silicium lors de la charge rapide, cette hypothèse peut mener à des erreurs significatives.

Dans ce travail, nous dérivons une expression de potentiel chimique valide pour la déformation rapide sans l'hypothèse de faible taux, indépendante de la morphologie de l'électrode. En appliquant le nouveau modèle aux électrodes nanoparticules de silicium sous contrôle potentiostatique et galvanostatique, nous quantifions l'impact sur la contrainte et la distribution de lithium. Nous identifions également un rayon critique où les deux modèles convergent, fournissant un aperçu des limites spatiales de l'influence de la contrainte.

Méthodes

Équations mécaniques

L'insertion de lithium induit un changement volumétrique. Nous décrivons la déformation en utilisant les coordonnées lagrangiennes et eulériennes :
$$\mathbf{U}=\mathbf{x}-\mathbf{X}\quad (1)$$ où U est le champ de déplacement. Le tenseur gradient de déformation est
$$\mathbf{F}=\frac{\partial \mathbf{x}}{\partial \mathbf{X}}=\mathbf{I}+\mathrm{Grad}\,\mathbf{U}\quad (2)$$

Pour une particule sphérique avec coordonnées (R,Θ,Φ) (lagrangiennes) et (r,θ,φ) (eulériennes), le tenseur gradient se simplifie à
$$\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)$$

La déformation est décomposée en changement de forme induit par lithiation et déformation élastique réversible :
$$\mathbf{F}=\mathbf{F}^e\mathbf{F}^c\quad (4)$$ avec F^c isotrope :
$$\mathbf{F}^c=(1+\Omega C)^{1/3}\mathbf{I}\quad (5)$$ où Ω est le volume molaire partiel. Par conséquent,
$$\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)$$

Le tenseur de déformation Green-Lagrange E est
$$\mathbf{E}=\frac{1}{2}\left(\mathbf{F}^T\mathbf{F}-\mathbf{I}\right)\quad (7)$$ donnant les composantes radiale et tangentielle
$$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)$$

En assumant l'élasticité linéaire, la densité d'énergie de déformation est
$$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)$$ où E_h et ν sont le module de Young et le coefficient de Poisson. La contrainte de Piola-Kirchhoff première suit
$$\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)$$

Par composantes,
$$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)$$

L'équilibre sans forces de volume requiert
$$\frac{\partial P_R}{\partial R}+\frac{2(P_R-P_{\Phi})}{R}=0\quad (16)$$ avec conditions aux limites
$$u(0,t)=0,\;P_R(R_0,t)=0\quad (17)$$

Équation de transport de masse

La concentration de lithium C(X,t) et le flux J(X,t) satisfont
$$\frac{\partial C}{\partial t}+\mathrm{Div}\,\mathbf{J}=0\quad (18)$$. Sous symétrie sphérique, cela se réduit à
$$\frac{\partial C(R,t)}{\partial t}+\frac{1}{R^2}\frac{\partial\big(R^2J(R,t)\big)}{\partial R}=0\quad (19)$$. Le flux suit la loi de Fick avec force motrice de potentiel chimique :
$$J=-\frac{CD}{R_gTF_{11}F_{11}}\frac{\partial \mu}{\partial R}\quad (20)$$.

Le potentiel chimique dérive de la densité d'énergie interne totale
$$\Pi=\varphi(C)+W(C,\mathbf{E}^e)\quad (22)$$, donnant
$$\mu=\frac{\partial\Pi}{\partial C}=\mu_0(C)+\tau(\mathbf{E}^e,C)\quad (23)$$ où τ est la partie dépendante de la contrainte. Pour une solution diluée,
$$\mu_0(C)=\mu_0+R_gT\ln(\gamma C)\quad (24)$$ avec γ=1.

Traditionnellement, l'approche d'énergie libre de Helmholtz donne
$$\tau_H=\frac{\partial W}{\partial C}= -\det(\mathbf{F}^e)\sigma_m\Omega\quad (25)$$ menant à
$$\mu={\mu}_0+R_gT\ln C-\det(\mathbf{F}^e)\Omega\sigma_m\quad (26)$$. Ici, σ_m est la contrainte hydrostatique de Cauchy
$$\sigma_m=\tfrac{1}{3}tr\,(\det^{-1}\mathbf{F}\mathbf{PF}^T)\quad (27)$$. Cette expression (26) est le modèle traditionnel de potentiel chimique.

Alternativement, en traitant l'énergie libre de Gibbs, on obtient
$$\tau_G=\frac{\partial W}{\partial C}\quad (28)$$ et par conséquent
$$\mu={\mu}_0+R_gT\ln C-\frac{\partial W}{\partial C}\quad (29)$$, notre modèle développé. Le système de transport de masse couple (19), (20), (26) et (29). Les conditions aux limites pour les opérations potentiostatiques (C(R₀,t)=C_max) et galvanostatiques (J(R₀,t)=j₀(1+u/R)^2) sont
$$C(R_0,t)=C_{\max},\;J(R_0,t)=j_0\big(1+u/R\big)^2\quad (30-31)$$ avec conditions initiales
$$C(R,0)=0\quad (32)$$.

Implémentation numérique

Les solutions analytiques sont intraitables ; nous avons résolu les EDPs couplées en utilisant COMSOL Multiphysics, adoptant les paramètres matériels listés Tableau 1. Des variables adimensionnelles furent utilisées pour la clarté. L'état de charge (SOC) est calculé comme
$$SOC=\frac{\int_0^{R_0}C(R,t)R^2dR}{\int_0^{R_0}C_{\max}R^2dR}\quad (33)$$. Les flux dérivés des deux modèles de potentiel chimique sont
$$J_H=\frac{\partial\tau_H}{\partial R},\;J_G=\frac{\partial\tau_G}{\partial R}\quad (34)$$.

Résultats et discussion

La Figure 1 illustre les profils de concentration de lithium, contrainte radiale et contrainte circonférentielle pour une nanoparticule sphérique de Si sous cyclage galvanostatique. Les modèles développé (lignes continues) et traditionnel (triangles) donnent des distributions de concentration presque identiques à tous les SOC. Les profils de contrainte diffèrent seulement légèrement, avec une élévation centrale modeste à 46,7 % et 65,5 % SOC.

La Figure 2 présente des données analogues pour le cyclage potentiostatique. Ici, la divergence entre les deux modèles est plus prononcée, particulièrement dans les contraintes centrales aux SOC intermédiaires. Parce que l'opération potentiostatique impose une concentration de surface constante, le taux de charge — et par conséquent le taux de déformation — est plus élevé qu'en mode galvanostatique. La déformation globale pour un SOC donné reste similaire, soulignant que le facteur clé est le taux de déformation plutôt que la déformation absolue.

La Figure 3 (de la Réf. [17]) démontre qu'une nanoparticule Si de 620 nm peut se dilater complètement en une minute à 2 V versus Li, soulignant la pertinence de la modélisation de déformation rapide. Dans de tels régimes, les deux modèles de potentiel chimique prédisent des états de contrainte markedly différents, bien que des données expérimentales de contrainte restent indisponibles pour validation quantitative.

La Figure 4 examine le rapport de flux J_H/J_G sous cyclage galvanostatique à divers courants de charge. Le rapport reste proche de l'unité à travers la particule, avec une légère augmentation près du centre et une diminution près de la surface. L'écart s'élargit avec un courant plus élevé, indiquant une sensibilité plus forte au taux de déformation. Toutes les courbes se croisent près d'un rayon critique — nommé région indépendante du potentiel chimique (CIR) — où les flux convergent indépendamment du modèle choisi. Le CIR se déplace plus près de la surface à mesure que le courant augmente.

Les Figures 5 et 6 comparent les distributions de contrainte hydrostatique de Cauchy pour les deux modèles. Dans le CIR, σ_m est presque nul pour la majeure partie de la période de lithiation, expliquant pourquoi les modèles s'accordent là. Hors CIR, σ_m diverge, amplifiant les différences de potentiel chimique et donc de flux.

Conclusions

Nous avons introduit une expression de potentiel chimique valide pour la déformation rapide, éliminant l'hypothèse de faible taux inhérente aux modèles traditionnels. La comparaison avec la formulation traditionnelle révèle des différences négligeables sous cyclage galvanostatique mais des écarts significatifs sous contrôle potentiostatique, pilotés par le taux de déformation. Un rayon critique — le CIR — émerge où les flux et contraintes sont virtuellement identiques, se rapprochant de la surface à des courants plus élevés. Ces résultats soulignent la nécessité d'utiliser le modèle développé pour une prédiction précise de l'évolution de la contrainte dans les électrodes nanoparticules de silicium lors de la charge rapide.