Shape Influence of Active Material Micro-Structure on Diffusion and Contact Stress in Lithium-Ion Batteries

Abstract

Electrochemical-mechanical modelling is a key issue to estimate the damage of active material, as direct measurements cannot be performed due to the particles nanoscale. The aim of this paper is to overcome the common assumptions of spherical and standalone particle, proposing a general approach that considers a parametrized particle shape and studying its influence on the mechanical stresses which arise in active material particles during battery operation. The shape considered is a set of ellipsoids with variable aspect ratio (elongation), which aims to approximate real active material particles. Active material particle is divided in two domains: non-contact domain and contact domain, whether contact with neighbouring particles affects stress distribution or not. Non-contact areas are affected by diffusion stress, caused by lithium concentration gradient inside particles. Contact areas are affected simultaneously by diffusion stress and contact stress, caused by contact with neighbouring particles as a result of particle expansion due to lithium insertion. A finite element model is developed in Ansys™APDL to perform the multi-physics computation in non-spherical domain. The finite element model is validated in the spherical case by analytical models of diffusion and contact available for simple geometry. Then, the shape factor is derived to describe how particle shape affects mechanical stress in non-contact and contact domains.

1. Introduction

Lithium-ion batteries (LIB) play a strategical role in the actual green energy conversion, indeed they are the most widespread storage system of electrical energy so far, with a huge field of application: from small electronics up to heavy-duty vehicles [1,2].

Nowadays, the main challenge is to increase the amount of energy stored in the battery, working with simulations [3] on new materials which guarantee a greater nominal energy density (amount of energy stored per unit of weight) and ensuring that the nominal properties, such as capacity, power and internal resistance, could be kept through the whole life cycle as well. Indeed, complex ageing mechanisms [4] take place in lithium-ion batteries, damaging the cell itself and causing a reduction in the nominal properties.

Besides the electrochemical ageing mechanisms, the mechanical ones have a crucial role in battery damage [5], and they can be identified at the active material particles scale. Electrode manufacturing [6] consists in the mixing of active material powder (major content), binder, conductive components and additives to get the so-called wet slurry, which is then deposited on current collector and dried. In this manner, particles of the active material are embedded in a solid porous matrix, which is filled with electrolyte when battery is assembled. Insertion and extraction of lithium ions in active material during battery operation cause an inhomogeneous ions concentration in the active material particles which results in swelling and differential strain, causing mechanical stress. Differential strain causes the so-called diffusion-induced stress [7], and swelling causes contact stress between neighbouring particles [8]. In some cases, stress in the active material can overcome the yield stress and lead to pulverization and breakages, causing isolation of some portions of active material and crack growth, which combined with SEI precipitation concur in capacity fade [9]. In particular, the analysis of active material micro-structure shows the presence of small flaws even in the fresh material, known as pre-cracking. The flaws in fresh material are the nucleation points of cracks [10] induced by the mechanical stress computed in this paper. Therefore, the analysis of flaws size in electrode microstructure [11] is necessary to estimate how mechanical stress affect fractures growth. Moreover, the electrode manufacturing process itself can cause flaws during drying, thus new processes are being studied [6].

The study of mechanical ageing in LIB dates back to the works of Christensen and Newman [12], which established one of the first diffusion-induced stress (DIS) models applied to LIB, followed then by other authors [13,14]. The weaknesses of this type of model are the assumptions of perfectly spherical geometry and standalone particle. A more realistic approach should take into account the influence of the particle shape on the stress caused by lithium diffusion, and the contact stress which comes from the interaction with other particles.

Some works [15,16,17,18,19,20] focused their attention in computing the stress on real particle geometry reconstructed from SEM images, but the computation on a singular particle shape cannot be representative for a comprehensive study, neither be extended as an equivalent shape of a certain active material. Other works [13,21,22,23] considered ellipsoid as a representative particle shape for electrochemical-mechanical modelling.

Recently, Lu et al. [8] established an analytical framework to study the stress caused by the contact between spherical active material particles due to the their expansion during lithiation process. They showed that contact is the predominant source of mechanical stress, when compared to diffusion induced stress. The effect of matrix confinement was studied in all solid state battery [24] as well as in traditional batteries according to the porous electrode theory [25,26], but contact stress is not explicitly computed.

The aim of this paper is to deepen the role of particle shape on stress amplification and to describe the types of stress which arise in active material.

The particle is split in two domains: contact domain, affected by both diffusion and contact stress, and non-contact domain, affected by just diffusion stress, far enough from contact point to not be affected.

A set of ellipsoids with a parametrized aspect ratio is considered as an approximation of real particle shape. Then, shape factor, which describes the stress amplification compared to the spherical case, is deduced as a function of aspect ratio to estimate how the particle shape amplifies mechanical stress.

Mechanical stress is computed with a multi-physics FE model in Ansys APDL, capable of computing the lithium concentration inside the particle and the stress due to concentration gradient and contact with neighbouring particles. The FE model is validated in spherical geometry with analytical models of diffusion stress studied in a previous paper [7] and with Hertz theory. A dilute assumption is made in this paper, so stress generation due to phase transition [27,28,29] is neglected so far.

2. Equivalent Particle Shape

Most of the works presented in literature concerning DIS analysis in LIB address the domain problem according to two approaches. Active material particles are considered spherical according to the former [9,12,13,14], which is an acceptable approximation, allowing easier calculation and exploiting even analytical solutions [7]. Instead, real particle shape reconstructed from SEM images is considered for DIS calculation according to the latter [15,16,17,18,19,20].

The problem related to the first approach is that it neglects the effects of particle shape on the stress computation, as only spherical particles are considered. The weakness of the second approach concerns the representativity of the results of an analysis conducted over a single real particle geometry: the peculiarities of a single particle are unlikely to be extended to all the particles of the active material, even more so to other active materials particles.

Then, an equivalent geometry which approximates the average particle shape of the active material of interest must be chosen for a general study, which can be later applied to specific cases.

Ellipsoid, constructed as an axisymmetric semi-ellipse, is adopted as equivalent geometry in this paper, because it is a good trade-off between the accuracy of particle approximation and the possibility to get a geometry totally described by few parameters. The geometry problem can be divided in two issues: size and shape. Consequently, ellipsoid is defined by two parameters: radius and aspect ratio. The first parameter is equal the radius of the sphere when aspect ratio is equal to 1. It controls the size issue, which was addressed in a previous paper [7], where was showed that greater particles suffer greater diffusion stress [13,20]. The shape problem is addressed in this paper, and it is described by aspect ratio ($\alpha$), defined as the ratio between maximum and minimum ellipse semi-axes, as shown in Figure 1. In this paper, a fixed radius equal to 5 μm and a variable aspect ratio are considered. The condition followed to compute axes length is that aspect ratio has to be respected and the surface of the ellipsoid has to be equal to the surface of the reference sphere ($\alpha =1$) with radius 5 μm. In this way the same amount of lithium, which is proportional to the battery state of charge (SOC), is stored in the particle with analysis carried out with the same simulation time and lithium flux (proportional to the (dis)charging rate of the cell). The particles of the active material of interest, captured with X-ray tomography [21] or using atomic force microscopy [19], are approximated choosing the proper aspect ratio and radius values [30], so that the equivalent ellipsoid circumscribes the average particle shape, as shown in Figure 1.

3. Multiphysics Modelling

Common active materials of LIBs, such as Lithium Manganese Oxide (LMO), Lithium Iron Phosphate (LFP), Lithium Nickel Cobalt Manganese Oxide (NMC) and Lithium Nickel Cobalt Aluminium Oxide (NCA) for cathode and Graphite, Silicon (with some exceptions) and Lithium titanate (LTO) for anode, are made of a matrix of micrometre particles which host and release lithium ions during battery operation [5,7,12,13,31].

Lithium ions cause a chemical strain when intercalated or deintercalated in active material particles, similar to thermal strain [7,13], which causes mechanical stress in turn. Lithium flux at the particle boundaries makes lithium ions diffuse from the surface to the core causing an inhomogeneous lithium concentration inside the particle. As active material particles expand proportionally to lithium content, the differential expansion of areas characterized by different lithium concentration causes the so called diffusion induced stress.

On the other hand, particles expansion causes the contact between neighbouring particles and gives rise to contact stress around the contact point due to the active material confinement.

The magnitude of diffusion-induced stress and contact stress have two different sources:

• Contact stress depends on the amount of inhibited displacement ($\delta$), proportional to the amount of expansion which would occur if the particle was able to expand freely. The expansion of the particle surface does not depends on the distribution of the lithium concentration inside active material, but just on the amount of lithium stored in particle (K). Qualitatively, it is reasonable that higher lithium content causes greater particle expansion, and then the contact with neighbouring particles is stronger as a consequence.
• Diffusion induced stress depends uniquely on the concentration gradient: high differences in lithium content in the particle makes great differential strain which lead to high stress.

3.1. DIS Model

The diffusion stress computation follows the model explained in a previous paper [7], whose main equations are reminded in

Table 1. Diffusion-induced stress (DIS) model equations [7].

Mechanical Equations
Constitutive	$${\sigma }_r=\frac{E\left[\left({\epsilon }_r-\frac{\mathsf{\Omega }c}{3}\right)(1-\nu )+2\nu \left({\epsilon }_c-\frac{\mathsf{\Omega }c}{3}\right)\right]}{(1+\nu )(1-2\nu )}{\sigma }_c=\frac{E\left[\left({\epsilon }_c-\frac{\mathsf{\Omega }c}{3}\right)+\nu \left({\epsilon }_r-\frac{\mathsf{\Omega }c}{3}\right)\right]}{(1+\nu )(1-2\nu )}$$ (1)
Congruence	$${\epsilon }_r=\frac{du}{dr}{\epsilon }_c=\frac{u}{r}$$ (2)
Equilibrium	$$\frac{d{\sigma }_r}{dr}+\frac{2}{r}\left({\sigma }_r-{\sigma }_c\right)=0$$ (3)
Diffusion Equation
Chemical potential	$$\mu ={\mu }_0+R_{g}Tln\left(c\right)-{\sigma }_h\mathsf{\Omega }$$ (4)
Flux	$$J=-Mc\nabla \mu =-D\nabla c+\frac{D\mathsf{\Omega }}{RT}C\nabla {\sigma }_h$$ (5)
Mass conservation	$$\frac{\partial c}{\partial t}+\nabla ·(\nabla J)=0$$ (6)
Solutions
Displacement	$$u\left(r\right)=\frac{\mathsf{\Omega }}{3(1-\nu )}\left[(1+\nu )\frac{1}{r^2}{\int }_0^{r}c\left(r\right)r^{2}dr+2(1-2\nu )\frac{r}{R^3}{\int }_0^{R}c\left(r\right)r^{2}dr\right]$$ (7)
Radial stress	$${\sigma }_r\left(r\right)=\frac{2\mathsf{\Omega }}{3}\frac{E}{1-\nu }\left[\frac{1}{R^3}{\int }_0^{R}c\left(r\right)r^{2}dr-\frac{1}{r^3}{\int }_0^{r}c\left(r\right)r^{2}dr\right]$$ (8)
Hoop stress	$${\sigma }_c\left(r\right)=\frac{\mathsf{\Omega }}{3}\frac{E}{1-\nu }\left[\frac{2}{R^3}{\int }_0^{R}c\left(r\right)r^{2}dr+\frac{1}{r^3}{\int }_0^{r}c\left(r\right)r^{2}dr-c\left(r\right)\right]$$ (9)
Concentration	$$c(r,t)=c_0+\frac{IR}{F_{a}D}\left[3\tau +\frac{1}{2}{\left(\frac{r}{R}\right)}^2-\frac{3}{10}-2\frac{R}{r}\sum _{n=1}^{\infty }\left(\frac{sin({\lambda }_{n}r/R)}{{\lambda }_n^{2}sin\left({\lambda }_n\right)}{e}^{-{\lambda }_n^2\tau }\right)\right]$$ (10)

.

3.2. Hertz Model

Contact stress is implemented in the model of Section 3.1 to take into account the influence of surroundings on particle expansion. Considering liquid electrolyte, the interactions between neighbouring particles are the ones which lead to noteworthy contact stress, as binder and conductive material are much softer (lower Young modulus) and do not contribute to contact stress appreciably.

Contact stress is computed analytically for spherical geometry following the Hertz approach, resumed in this section. Then, the analytical theory validates the FE model algorithm embedded in the complete model in Section 3.3.

Particle expansion due to lithium insertion causes contact between particles, this means that the magnitude of the contact stress is proportional to the amount of lithium stored in the particle, as the greater the amount of lithium inserted, the greater the particle expansion. Then, contact stress is driven by the particle surface displacement, which depends on the parameters of the electrochemical simulation (simulation time and lithium flux), and the active material properties reported in

Table 2. Chemical and mechanical properties of Li_xMn₂O₄ (LMO) [7].

Parameter	Symbol	Value	Unit
Diffusion coefficient	D	$7.08\times {10}^{-15}$	m${}^2$/s
Partial molar volume	$\mathsf{\Omega }$	$3.497\times {10}^{-6}$	m${}^3$/mol
Maximum concentration	$C_{max}$	$2.29\times {10}^4$	mol/m${}^3$
Young modulus	E	10	GPa
Poisson ratio	$\upsilon$	0.3	-
Particle radius (sphere)	R	$5\times {10}^{-6}$	m

.

Referring to Figure 2, contact stress is caused by the amount of displacement $\delta$ inhibited by the neighbouring particles. As a first approximation, the boundary of the active material is considered infinitely rigid; this means that the active material layer does not expand and the centre of the particles is fixed, referring to the case $\beta =1$ in the upper part of Figure 2. In this scenario, the displacement $\delta$, prevented by the surroundings, is identically equal to the displacement of the particle surface in free expansion condition ($u_{surf}$). On the contrary, if the edges of the active material layer expand, the centre of the spheres translates of a quantity $u_o$ to accommodate the particles deformation. This behaviour is described by the parameter $\beta$ in Equation (11), which is the ratio between the amount of displacement prevented by contact and the theoretical surface displacement due to lithium insertion without external constraints.

$$\beta =\frac{\delta }{u_{surf}}$$

(11)

The case $\beta =0.5$, which occurs when the displacement of the centre is equal to contact approach $\delta$, is shown in the lower part of Figure 2.

The case $\beta =1$ is considered for sake of simplicity in all the simulations in this paper, but in general the prevented displacement is $\delta =\beta u_{surf}$.

The displacement solution in free expansion of Table 1 is reported in Equation (13) as a function of dimensionless radius $\tilde{r}$.

$$u\left(\tilde{r}\right)=\frac{1+\nu }{1-\nu }\frac{\mathsf{\Omega }R{c}_{max}}{3}\frac{1}{{\tilde{r}}^2}{\int }_0^{\tilde{r}}\tilde{c}{\tilde{r}}^{2}d\tilde{r}+\frac{1-2\nu }{1-\nu }\frac{2\mathsf{\Omega }R{c}_{max}}{3}\tilde{r}{\int }_0^1\tilde{c}{\tilde{r}}^{2}d\tilde{r}$$

(12)

where $\tilde{c}=c/c_{max}$ and $\tilde{r}=r/R$. Then, the displacement of the sphere surface ($\tilde{r}=1$) is reported in Equation (13).

$$u_{surf}=\mathsf{\Omega }R{c}_{max}{\int }_0^1\tilde{c}{\tilde{r}}^{2}d\tilde{r}=\frac{\mathsf{\Omega }R{c}_{max}}{3}K$$

(13)

where $K=3{\int }_0^1\tilde{c}{\tilde{r}}^{2}d\tilde{r}$ gives the percentage of lithium content in active material particle: it is equal to zero when no lithium is in the particle, and it is equal to one when lithium concentration approaches the solubility limit. K is the result of the electrochemical simulation and depends on the lithium flux applied on the particle surface (proportional to the C-rate) and on the simulation time (linked to the (dis)charging time), as well as on the electrochemical properties in Table 2.

Assuming the contact between two particles as sketched in Figure 3a, the equivalent configuration in Figure 3b is studied, where the sphere, characterized by the equivalent mechanical properties and radius of Equation (14), comes in contact with a rigid half space surface.

$$\left\{\begin{array}{c}\frac{1}{E^{*}}=\frac{1-{\nu }_1^2}{E_1}+\frac{1-{\nu }_2^2}{E_2}\hfill \\ R^{*}=\frac{1}{R_1}+\frac{1}{R_2}=\frac{R_1{R}_2}{R_1+R_2}\hfill \end{array}\right.$$

(14)

For our purpose, the two spheres in contact are supposed to be equivalent, so $E_1=E_2=E,{\nu }_1={\nu }_2=\nu ,R_1=R_2=R$. The normal contact approach $\delta$ represents the displacement of the surface due to contact stress compared to free condition, and it is expressed in Equation (15), referring to Figure 3.

$$\delta =u_y\left(w\right)+y\left(w\right)=u_y\left(w\right)+\frac{w^2}{2R^{*}}$$

(15)

where $u_y\left(w\right)$ is the elastic displacement of the surface, and $y\left(w\right)$ is the original surface separation inside the contact area. According to Hertz theory [32], the surface displacement $u_y$ in the contact area is reported in Equation (16).

$$u_y\left(w\right)=\frac{\pi P_h}{4a{E}^{*}}(2a^2-w^2)$$

(16)

The values of the elastic displacement $u_y$ at the centre and at the surface of the sphere are replaced in Equation (17) using Equation (15).

$$\left\{\begin{array}{c}{u}_y\left(0\right)=\delta ⟶\delta =\frac{\pi a{P}_h}{2E^{*}}\hfill \\ u_y\left(a\right)=\delta -\frac{w^2}{2R^{*}}⟶a=\frac{\pi a{R}^{*}}{2E^{*}}\hfill \end{array}\right.$$

(17)

Then, the radius of contact area (a) and maximum contact pressure ($P_h$) are derived in Equations (18) and (19), respectively, solving the set of Equation (17).

$$a=\sqrt{\delta R^{*}}=\sqrt{\beta u_{surf}{R}^{*}}$$

(18)

$$P_h=\frac{2E^{*}}{\pi }\frac{a}{R^{*}}=\frac{2E^{*}}{\pi }\frac{\sqrt{\beta u_{surf}{R}^{*}}}{R^{*}}$$

(19)

Contact force is derived in Equation (20) from the maximum pressure $P_h$.

$$F=\frac{2}{3}\pi a^2{P}_h$$

(20)

In the end, contact stresses inside the elastic body are computed in Equation (21) [32], considering the sketch of Figure 3 for evaluating the Cartesian tensor. In general, the first principal stress is equal to the two circumferential stresses, and the third principal stress is the radial one.

$$\left\{\begin{array}{c}{\sigma }_1={\sigma }_2={\sigma }_x={\sigma }_z=-P_h\left[\left[1-{\zeta }_{a}t{g}^{-1}\left(\frac{1}{|{\zeta }_a|}\right)\right](1+\nu )-\frac{1}{2(1+{\zeta }_a^2)}\right]\hfill \\ {\sigma }_3={\sigma }_y=-P_h\frac{1}{1+{\zeta }_a^2}\hfill \end{array}\right.$$

(21)

where ${\zeta }_a$ is defined as ratio between the depth below the contact point and the radius of the contact area: ${\zeta }_a=\tilde{y}/a$.

3.3. FE Multiphysics Model

In this section, mechanical stress and lithium concentration are computed simultaneously in a coupled electrochemical-mechanical FE model in Ansys APDL. The model is implemented according to the migration model (tb,migr command in Ansys APDL), considering lithium concentration and hydrostatic stress as driving forces of lithium diffusion. Therefore, the model is a coupled multiphysics problem, as lithium concentration, governed by diffusion equations (Equations (4)–(6)), causes mechanical stress (Equations (1) and (21)), that in turn affects lithium diffusion and thus concentration distribution [7]. The geometry is modelled with the bidimensional plane element “Plane223” in axisymmetric mode, which have three degrees of freedom according to the structural-diffusion nature of the problem: displacements along x and y and lithium concentration. Contact elements “Conta172” and “Targe 169” are used to simulate contact between the particles surfaces.

The axisymmetric section of two half particles in contact is modelled as shown in Figure 4, thanks to the symmetry around mid-plane. Figure 4b shows the reason behind the division in contact and non-contact domains: contact alters the stress distribution in the contact domain, even away from the contact area. On the contrary, the non-contact domain is far enough from contact point and the stress distribution due to diffusion is not affected. The model is mechanically constrained so that rigid body motion is prevented and the condition of infinity stiffness of active material layer ($\beta =1$) is satisfied, as showed in Figure 4a. Finally, constant lithium flux is applied to particles surfaces as electrochemical boundary condition to simulate galvanostatic operation.

4. Results

The validation of FE multiphysics model with analytical models in spherical geometry is shown in this section, then the results in contact and non-contact domain are presented along with the shape factor which highlights the effect of particle shape on mechanical stress.

4.1. FE Sub-Models Validation

First, the computation of diffusion induced stress with the dedicated FE sub-model is validated. The results got through FE modelling in spherical geometry without contact contribution are compared with the analytical results of a previous study [7] in Figure 5 at different lithium flux. The results of Figure 5 validate the FE sub-model, which is expected to estimate correct results even in non-spherical geometry. Furthermore, the effect of boundary lithium flux is pointed out: greater flux causes steeper concentration gradient inside the particle (Figure 5a), which causes greater diffusion stress in turn (Figure 5b).

Then FE contact sub-model is validated following the approach of Figure 6a: the model is similar to the general multiphysics model of Section 3.3, but this sub-model predicts just contact stress related to a specific electrochemical condition.

The main differences with the model of Section 3.3 are the mechanical constraints and the external input. The constraints of the sub-model prevents rigid body motion and allow the displacement along y direction, where contact occurs.

The input of this model is the displacement $\delta =\beta u_{surf}$ imposed at the particle surface, which results from the electrochemical-mechanical simulation of Table 1. Contact force caused by displacement $\delta$, derived with the Hertz model of Section 3.2, is assigned both to FE sub-model and analytical model as boundary condition, instead of the lithium boundary flux. Actually, a distributed equivalent pressure $P_{eq}$ so that $P_{eq}=\frac{F}{\pi R^2}$ is assigned to the upper edge of the hemisphere which is free to move in vertical direction, as shown in Figure 6a. As pointed out in Section 3.2, contact stress depends just on the particle expansion $u_{surf}$, which depends on the amount of lithium ions stored in the particle (K), which is the result of the combination of lithium flux and simulation time.

Contact is simulated with penalty function algorithm, with a contact stiffness equal to 5 N/mm${}^3$, that, after some tries, it is found to guarantee no surfaces penetration and a reasonable computing time.

The comparison between analytical and FE sub-model results in Figure 6b are satisfactory and validate the contact algorithm of the FE sub-model.

Following this procedure, the contact FE sub-model gives two important results:

a

The contact algorithm is validated by the analytical model of Section 3.2 as shown in Figure 6b, so the effectiveness of the model is confirmed and can be incorporated in the electrochemical-mechanical FE model which aims to compute diffusion induced and contact stress simultaneously.

b

The contact stress due to a specific electrochemical operation is computed, independently from diffusion induced stress.

4.2. Coupling Effect

The multiphysics problem of Section 3.3 can be expressed according two approaches: coupled and uncoupled, whether mechanical stress influence on lithium diffusion is considered or not. The coupling resides in the expression of lithium flux (Table 1), and it is expressed in the uncoupled and coupled form in Equations (22) and (23) accordingly.

$$\begin{array}{cc}\hfill J=& -D\frac{\partial c}{\partial r}\end{array}$$

(22)

$$\begin{array}{cc}\hfill J=& -D\frac{\partial c}{\partial r}+\frac{D\mathsf{\Omega }}{RT}c\frac{\partial {\sigma }_h}{\partial r}\end{array}$$

(23)

Concentration gradient is the only driving force in the uncoupled formulation. On the contrary, lithium diffusion depends on the gradient of concentration and hydrostatic stress in the case of coupled problem:

• Lithium diffuses from areas with higher concentration to areas with lower concentration. For this reason, lithium diffuses towards the core of the particle (negative lithium flux in Equation (23)) when the derivative of concentration is positive, according to Equation (23).
• Lithium diffuses from areas with low hydrostatic stress to areas with higher hydrostatic stress, assuming negative stress as compression and positive as tensile. Namely, the diffusion of lithium ions is enhanced from compression areas to tensile (or lower compression) areas. For this reason, a negative derivative of hydrostatic stress enhances lithium diffusion inside the particle (negative lithium flux) according to Equation (23).

The hydrostatic stress is computed in non-contact areas and in contact areas together with its derivative with respect to radial coordinate in Figure 7a,b. The plot of Figure 7a highlights the stress contribution to lithium diffusion in non-contact and contact areas. The slight increase in hydrostatic stress from the surface towards the core of the particle when contact does not occur slightly enhances lithium diffusion all over the particle with respect to the uncoupled model [7], since the stress derivative is slightly negative, as shown by the comparison of blue and black dashed lines in Figure 7c.

On the other hand, the sharp hydrostatic stress gradient close to the contact point ($0.8\le \tilde{r}\le 1$) causes a huge decrease of stress derivative (Figure 7b), then hydrostatic stress gives a substantial contribution to lithium flux in this area, according to Equation (23). The greater lithium flux in the external part of the particle causes an increase of lithium concentration in this area, as shown in Figure 7c.

The coupling between mechanical stress and lithium diffusion, which was already remarkable with diffusion stress alone [7], becomes substantial in the areas where contact occurs. Indeed, the steeper concentration gradient close to contact point causes higher diffusion stress, and the higher lithium content causes a greater expansion of particle surface, which results in greater contact stress in turn.

Von Mises stress in contact areas increases steadily as lithium insertion goes on in time because it is strongly influenced by contact stress, as shown in Figure 8. Indeed, contact stress is proportional to particle expansion that depends on lithium content in turn, while diffusion stress keeps nearly constant through the whole SOC window after an initial transient. The differences from the coupled and uncoupled model are highlighted in Figure 8a, which shows that the influence of mechanical stress on diffusion causes higher Von Mises stress, as qualitatively explained above. Furthermore, the discrepancy between the two models increases with time, as the coupling becomes stronger because of the increase of mechanical stress. Moreover, the comparison between the results of multi-physics model and contact model of Section 3.2 in Figure 8b shows that diffusion stress gives a remarkable contribution to total stress in contact domain. This makes the multiphysics model meaningful, as the stress state cannot be deduced by considering contact and diffusion stress separately, as supposed in previous works [8].

4.3. Shape Factor

Shape factor is presented in this section to understand how the particle shape affects mechanical stress. Shape factor is computed as the ratio between Von Mises stress in ellipsoids and in the reference sphere. A set of aspect ratio values from $\alpha =0.3$ to $\alpha =3.5$ are considered. Prolate ellipsoids are characterized by aspect ratio $0.3\le \alpha <1$, oblate ellipsoids by $1<\alpha \le 3.5$, and $\alpha =1$ refers to sphere. Particles in active material are assumed to have not so flattened shapes to be described by aspect ratios lower than 0.3 or greater than 3.5.

Firstly, non-contact domain is considered, then all the following considerations hold for the areas far enough from contact point not to be affected. Concentration distribution in sphere and ellipsoids are reported in Figure 9a–c, respectively. An inhomogeneous lithium distribution at the edge of the particle is observed in ellipsoids, which increases with aspect ratio, as pointed out also by other works [13,22,23]. The solubility limit is reached earlier in the more convex areas, as confirmed by the comparison of Figure 9a–c: the concentration value is up 33% higher in the more convex area of ellipsoids with the same simulation time.

As a consequence, even displacements are not symmetric any more: convex areas suffer greater displacements compared to flat areas due to higher lithium content.

Particle shape influences Von Mises stress as well, as highlighted in Figure 10. Ellipsoid breaks the symmetry of stress distribution and makes some areas characterized by higher stress and other by lower stress compared to spherical geometry. Moreover, as the aspect ratio increases over 1, the maximum stress moves from less convex areas to more convex ones. On the contrary, peak stress is observed in flat areas in prolate ellipsoids ($\alpha <1$).

DIS shape factors are displayed in Figure 11 as the ratio between the maximum Von Mises stress computed in ellipsoid and in sphere in non-contact domain as a function of the aspect ratio. DIS shape factor is computed for x axis and y axis referring to Figure 10, which means that the maximum Von Mises stress along x direction and y direction is picked for the computation.

Flat areas in oblate ellipsoids with aspect ratio between 1 and 2 and in the corresponding prolate ellipsoids with aspect ratio between 0.5 and 1 are the most critical conditions, since a stress amplification over 20% is detected by the DIS shape factor along y and x directions, respectively, as shown in Figure 11a,b. The DIS shape factor can be fitted by a parable in these regions, then it decreases linearly outside. This trend can be deduced also by previous works [13,23], which detects greater Von Mises stress for aspect ratios between 1 and 2 in oblate ellipsoids.

On the other hand, DIS shape factors are lower than 1 in all the other regions in more convex areas, because the stress induced by diffusion is mitigated by the longer diffusion path.

The simulations of Figure 11 are conducted until a SOC of about 30%, because diffusion stress, after an initial transient, does not change beyond this value, and accordingly neither does the DIS shape factor.

As noticed before, with increasing aspect ratio over 1.5, the peak stress moves from less convex areas to more convex areas, then shape factor along y decreases and the one along x increases.

DIS shape factors have a weak dependence on lithium boundary flux, which controls the absolute magnitude of the stress induced by diffusion, but the trends remain unchanged.

Finally, shape factor is computed in contact domain in Figure 12.

Overall shape factor in contact domain is highly influenced by contact stress as showed by the comparison of red line and dotted purple line in Figure 12, and thus it depends mainly on the size of contact area: prolate ellipsoid suffers greater contact stress because of the smaller contact area, oblate ellipsoid vice versa. The dotted purple line is obtained from the FE sub-model of Section 3.2, it is time-independent and satisfactory fitted by the hyperbolic function in Equation (24), with the constants $a=0.9$ and $b=0.14$ for the case reported in Figure 6.

$$S{F}_c=\frac{j}{\alpha }+k$$

(24)

However, some differences are highlighted between the contact and complete model due to the asymmetry of displacement in ellipsoid respect to sphere, already pointed out in Section 2: convex areas expand more than flat areas due to lithium diffusion; thus, the greater the particle expansion, the greater the contact stress. This fact causes higher contact stress in the areas where the convexity is greater than in the spherical case of reference (aspect ratio lower than one), which are already characterized by higher contact stress because of the smaller contact area. On the contrary, the lower expansion of flat areas, compared to the spherical case of reference, reduces the contact stress, which, however, is already mitigated by the larger contact surfaces.

These differences are highlighted by the comparison of red line and purple line in Figure 12: the overall shape factor is greater than the contact one when aspect ratio is lower than one, vice versa for aspect ratio greater than one. This trend makes the shape of the particle emphasize the total stress, intensifying the stress in convex areas in particular. However, overall shape factor is still well fitted by a hyperbolic curve, as shown by the comparison of red line and green dashed line in Figure 12. Interestingly, shape factor basically does not depend on lithium flux and lithium content.

5. Conclusions

The aim of this paper is to overcome the assumptions of spherical and standalone active material particle in the electrochemical-mechanical modelling of lithium-ion battery.

The geometry adopted is a parametrized ellipsoid: the parameter “radius” controls the size and “aspect ratio” controls the shape (elongation). The two parameters can be fitted to approximate real active material particles from X-ray tomography or atomic force microscopy.

Shape factor describes how particle shape affects mechanical stress in active material. It is computed as the ratio between the maximum Von Mises stress in ellipsoid and in the spherical geometry of reference, as a function of aspect ratio.

Mechanical stresses in active material particle are computed with a multiphysics model in Ansys APDL which solves the electrochemical and structural problem simultaneously. The model consistency is guaranteed by the validation of diffusion and contact sub-model with analytical solution available in spherical geometry reported in [7,32] respectively.

The active material particle is split in two domains: non-contact domain and contact domain, whether the area is affected by contact with neighbouring particles or not. Diffusion stress affects the non-contact domain, and it is the result of the inhomogeneous distribution of lithium ions inside active material particle due to diffusion. The combination of diffusion stress and contact stress affects the areas in contact domain because of the expansion of the particles in a confined active material layer during lithiation.

The multi-physics model is derived according to a coupled approach to capture the influence of mechanical stress on electro-chemistry. Namely, lithium diffusion is driven by concentration gradient and hydrostatic stress gradient. The coupling affects the results considerably in the non-contact areas, as already pointed out in a previous work [7], and has a substantial influence in contact domain, where contact causes a steep hydrostatic stress gradient close to the particle surface.

Shape factor in non-contact domain highlights that an increment of diffusion stress up to 20% occurs in oblate ellipsoids with aspect ratio between 1 and 2 and in prolate ellipsoids with aspect ratio between 0.5 and 1, which are the critical shapes that cause stress intensification. All the other shapes cause a considerable stress relaxation up to 80%.

Shape factor in contact domain follows an hyperbolic trend and highlights that convex areas (aspect ratio lower than 1) are characterized by a serious stress amplification, up to five time compared to the spherical case of reference because of the smaller contact areas and greater swelling due to lithium inhomogeneity. On the contrary, flat areas show a stress relaxation because of the larger contact area. Shape factor in the contact domain is basically independent on lithium flux and lithium content.

It can thus be concluded that shape factor analysis suggests the shape characteristics which guarantee a mitigation of mechanical stress and ageing consequently, providing a guidance for electrode and active material manufacturing.

From the results of the simulation displayed in this paper is possible to predict that particle size effect has an opposite influence on contact and non-contact domains. Greater particle size increases diffusion stress in non-contact domain because of the greater diffusion length. On the other hand, greater particle size makes the larger contact area mitigate the contact stress in contact domain.