Modeling Anisotropic Transport in Polycrystalline Battery Materials

Abstract

Hierarchical structures of many agglomerated primary crystals are often employed as cathode materials, especially for layered-oxide compounds. The anisotropic nature of these materials results in a strong correlation between particle morphology and ion transport. In this work, we present a multiphase-field framework that is able to account for strongly anisotropic diffusion in polycrystalline materials. Various secondary particle structures with random grain orientation as well as strongly textured samples are investigated. The observed ion distributions match well with the experimental observations. Furthermore, we show how these simulations can be used to mimic potentiostatic intermittent titration technique (PITT) measurements and compute effective diffusion coefficients for secondary particles. The results unravel the intrinsic relation between particle microstructure and the apparent diffusivity. Consequently, the modeling framework can be employed to guide the microstructure design of secondary battery particles. Furthermore, the phase-field method closes the gap between computation of diffusivities on the atomistic scale and the effective properties of secondary particles, which are a necessary input for Newman-type cell models.

1. Introduction

Layered oxides are an important material class for lithium as well as sodium ion intercalation batteries [1]. Their stoichiometry is described by the generalized formula A${}_X$MO${}_2$, where A is the intercalated ion and M stands for one or various transition metals such as Ni, Mn, or Co [2]. The crystal structure consists of MO${}_2$ sheets with A-ions between these layers, thus resulting in strongly anisotropic material properties on the single-crystalline scale. Ion diffusion is limited to the layers in a defect-free structure due to the higher energetic barrier for cross-layer hopping. This effect is mediated through crystal defects; however, diffusivities still vary by orders of magnitude between the directions normal and parallel to the layers in a real system, e.g., in LiCoO${}_2$ [3].

Apart from the crystal structure, ion transport in electrode particles is strongly influenced by the morphology (i.e., structural appearance on the microscale). Layered-oxide compounds based on nickel, manganese, and cobalt (NMC) are state-of-the-art cathode materials for commercial lithium-ion cells in high-capacity applications such as electric vehicles [4]. Typically, NMCs form hierarchical structures of hundreds of primary crystals agglomerated in a spherical secondary particle with a diameter in the range of 5–10 µm, sometimes referred to as a meatball structure [5,6]. Similar secondary morphologies can be obtained for a wide range of layered-oxide cathode materials based on lithium and sodium as the inserting ion [4,7,8,9]. Ionic transport can be altered in polycrystalline materials through microstructure design such as tailoring crystallographic orientations [7,8] or introducing nanoporosity [10]. Particles consisting of radially aligned grains show excellent rate performance and improved cyclic stability [7,8]. This can be attributed to less tortuous pathways for ion transport and the increase of small-angle grain boundaries, which lead to an overall reduction of accumulated stress in grain boundaries. The same applies for the nano-porous particle morphology, as the open-pored structure significantly reduces diffusion paths within the active material compared to close-packed spherical structures.

Secondary particles tend to form intergranular cracks if the primary crystals undergo significant expansion and contraction [5]. This effect is strongly interconnected with the anisotropic nature of these materials because layered oxides typically exhibit much stronger changes in the c-axis (perpendicular to the transition metal oxide sheets). As the primary crystals are typically randomly oriented, stresses arise due to the misorientation of neighboring grains that can lead to fracture and eventually disintegration of the secondary structure [5,6,11]. Accordingly, small-angle grain boundaries are not strained as much because neighboring grains expand in a similar way.

NMC agglomerates have been modeled by Xu et al. [6] to study the formation of intergranular cracks. Even though the polycrystalline nature is accounted for, isotropic input parameters for diffusivity, Young’s modulus and volumetric expansion, which do not reflect the physical properties on the crystalline scale, were used. Additionally, a detailed investigation of the chemo-mechanical bulk behavior of NMC811 single crystals has been conducted by Lim et al. [12] in which anisotropic stiffness tensors and lattice expansions were computed from density functional theory (DFT) and, subsequently, employed in phase-field simulations. The mechanical properties from DFT [12] were, furthermore, used to study simplified secondary agglomerates consisting of 67 primary particles by Taghikhani et al. [13]. Works on other battery materials underline the importance to account for material anisotropy in terms of ion diffusion and elastic deformation [14,15,16,17].

In this work we propose the multiphase-field (MPF) method as a suitable modeling framework for the investigation of polycrystalline battery materials. Originating from the modeling of solidification [18], the MPF method was successfully applied to study microstructure evolution in polycrystalline systems such as martensitic transformations [19]. Large amounts of phases can be handled efficiently via a locally reduced order parameter scheme [20]. Computational efficiency for large structures is, furthermore, ensured through parallel computing based on message passing interface (MPI) [21]. The application to intercalation in polycrystalline battery materials has been demonstrated in our previous work [11]. A new development from previous studies is the inclusion of strongly anisotropic diffusion coefficients in microstructures composed of randomly oriented grains. While the formulation of anisotropic diffusion in terms of a partial differential equation is rather straightforward, an accurate numerical implementation, especially in the non-grid-aligned case, is surprisingly difficult. Strongly anisotropic fluxes also play a role in the simulation of plasma physics [22,23]. The heat flux in magnetized plasmas exhibits a ratio of heat conductivity coefficients in the order of ${10}^9$. Günter et al. [22] highlight the accuracy of a symmetric finite difference scheme in comparison with other numerical implementations especially at high anisotropy ratios. The discretization has been successfully used in subsequent works [24] and has been applied for the computation of anisotropic fluxes [25]. Based on these works, we extend the simulation framework presented in Daubner et al. [11] to account for anisotropic fluxes in hierarchical battery particles.

2. Materials and Methods

2.1. Multi-Phase Diffusion

The bulk diffusion of lithium ions inside intercalation compounds can be modeled in a similar manner as interstitial diffusion on a sublattice. In the following, all equations are expressed in terms of the local site filling fractions c. The general form of the diffusion equation that accounts for mass conservation can be expressed as

$$\begin{array}{c}\hfill \frac{\partial c}{\partial t}+\mathbf{\nabla }·\mathit{J}=0\end{array}$$

(1)

where $\mathbf{\nabla }·(·)$ denotes the divergence operator. According to linear irreversible thermodynamics, the mass flux is taken as $\mathit{J}=-\mathit{M}\left(c\right)\mathbf{\nabla }\mu$. Generally, the diffusion potential $\mu$ can include electrical, chemical, as well as mechanical effects. The gradient $\mathbf{\nabla }\mu$ is the thermodynamic driving force for diffusion in the system. The ion mobility $\mathit{M}$ depends on the site filling fraction c, and furthermore, on the crystal structure of the material.

In the context of the multiphase-field method, various phases $\alpha \dots N$ are described by their local volume fractions $\mathit{\varphi }=\{{\varphi }_{\alpha },{\varphi }_{\beta },\dots ,{\varphi }_N\}$. In this work, $\mathit{\varphi }$ is used to differentiate between grains with different orientations but with equal chemical properties. Interfaces between neighboring grains are represented by diffuse transition regions where the order parameter ${\varphi }_{\alpha }$ varies between 0 and 1. The initial polycrystalline microstructures are created by solving the corresponding multiphase-field evolution equations [11,26,27] including only interfacial forces to create close-to-equilibrium interfaces. This work is based on a multi-obstacle potential and a gradient formulation based on dual interactions [26]. More details can be found in our previous work [11]. For the remainder of this study, we assume that no phase transformations occur during (dis-)charge. Future works could encompass the effects of cation disorder and surface densification, which has been observed in high-nickel NMCs [4] and LiNiO${}_2$ [28] The ion diffusion in a multi-phase system is given by [29]

$$\frac{\partial c}{\partial t}=\mathbf{\nabla }·\left(\mathit{M}(c,\mathit{\varphi })\mathbf{\nabla }\mu \right)$$

(2)

where $\mathit{M}={\sum }_{\alpha }^N{\mathit{M}}^{\alpha }{\varphi }_{\alpha }$ denotes the average mobility. The lithium content in interfaces is described by an interpolation of phase-specific concentrations $c={\sum }_{\alpha }^N{c}^{\alpha }{\varphi }_{\alpha }$ [29]. In the bulk of a grain, Equation (2) reduces to

$$\frac{\partial c^{\alpha }}{\partial t}=\mathbf{\nabla }·\left({\mathit{M}}^{\alpha }\left(c^{\alpha }\right)\mathbf{\nabla }\mu \right)=\mathbf{\nabla }·\left({\mathit{M}}^{\alpha }\left(c^{\alpha }\right)\frac{\partial \mu }{\partial c^{\alpha }}\mathbf{\nabla }{c}^{\alpha }\right)=\mathbf{\nabla }·\left({\mathit{D}}^{\alpha }\left(c^{\alpha }\right)\mathbf{\nabla }{c}^{\alpha }\right)$$

(3)

which yields the relation ${\mathit{M}}^{\alpha }={\mathit{D}}^{\alpha }/(\partial \mu /\partial c^{\alpha })$ between mobility and chemical diffusivity ${\mathit{D}}^{\alpha }$. We employ this relation to correlate the ion mobility in simulations with measured chemical diffusivities as $(\partial \mu /\partial c^{\alpha })$, which is given from the free energy expression Equation (5). For most applications, diffusion is assumed to be ideal and isotropic, which results in a scalar diffusivity $D^{\alpha }$ that is phase-wise constant. More generally, anisotropic diffusion is described by a second order tensor ${\mathit{D}}^{\alpha }$. Under the assumption that the diffusion potentials of co-existing phases are equal at any point within the diffuse interface (quasi-equilibrium condition [26]), Equation (2) can be re-written in terms of an evolution equation for the diffusion potential [29]

$$\begin{array}{c}\hfill \frac{\partial \mu }{\partial t}={\left[\sum _{\alpha }^N\frac{\partial c^{\alpha }}{\partial \mu }{\varphi }_{\alpha }\right]}^{-1}\left(\mathbf{\nabla }·\left(\mathit{M}(c,\mathit{\varphi })\mathbf{\nabla }\mu \right)-\sum _{\alpha }^N{c}^{\alpha }\frac{\partial {\varphi }_{\alpha }}{\partial t}\right)\end{array}$$

(4)

which is computationally beneficial. Equation (4) is compatible with the simulation of phase evolution as the last term accounts for mass conservation at moving interfaces. This has been used in our previous study [11] while in this context, all phases and grain boundaries are assumed to be stationary ($\partial {\varphi }_{\alpha }/\partial t=0,\forall \alpha$).

2.2. Anisotropic Diffusion Tensor

As discussed in the introduction, layered-oxides exhibit strongly anisotropic properties in terms of diffusivity and lattice expansion. For orthorhombic materials with strong anisotropy, the diffusion tensor in the lattice coordinate system can be described by a diagonal tensor with three independent entries. The layered structure of NMC results in transverse isotropy, i.e., in-plane diffusion is isotropic with a diffusion coefficient of $D_{\left|\right|}$, while the diffusivity along the c-axis direction $D_{\perp }$ exhibits a much lower value. The diffusion tensor in material coordinates is given by

$$\begin{array}{c}\hfill {\mathit{D}}_{123}=\left[\begin{array}{ccc}{D}_{\left|\right|}& 0& 0\\ 0& D_{\left|\right|}& 0\\ 0& 0& D_{\perp }\end{array}\right].\end{array}$$

Every grain $\alpha$ has an individual orientation that can be described by Euler angles or quaternions with respect to the global coordinate system. Rotation of the diffusion tensor from the material to the reference coordinate system can be achieved via rotation matrices, e.g.,

$${\mathit{R}}_x^{\alpha }=\left[\begin{array}{ccc}1& 0& 0\\ 0& cos{\theta }^{\alpha }& -sin{\theta }^{\alpha }\\ 0& sin{\theta }^{\alpha }& cos{\theta }^{\alpha }\end{array}\right],{\mathit{R}}_y^{\alpha }=\left[\begin{array}{ccc}cos{\psi }^{\alpha }& 0& sin{\psi }^{\alpha }\\ 0& 1& 0\\ -sin{\psi }^{\alpha }& 0& cos{\psi }^{\alpha }\end{array}\right],{\mathit{R}}_z^{\alpha }=\left[\begin{array}{ccc}cos{\phi }^{\alpha }& -sin{\phi }^{\alpha }& 0\\ sin{\phi }^{\alpha }& cos{\phi }^{\alpha }& 0\\ 0& 0& 1\end{array}\right]$$

which, in the XYZ-convention, yields the rotation transformation

$$\begin{array}{c}\hfill {\mathit{D}}_{xyz}^{\alpha }={\mathit{R}}_z^{\alpha }{\mathit{R}}_y^{\alpha }{\mathit{R}}_x^{\alpha }{\mathit{D}}_{123}{\left({\mathit{R}}_x^{\alpha }\right)}^T{\left({\mathit{R}}_y^{\alpha }\right)}^T{\left({\mathit{R}}_z^{\alpha }\right)}^T.\end{array}$$

Alternatively, the rotation can be computed based on quaternions that are more robust for arbitrary rotations [30,31] and, moreover, allow for the straight-forward interpolation of intermediate configurations. Effective material properties are computed by a mixture rule, i.e., interpolation of the phase-specific properties. Throughout this work, we employ a linear interpolation of diffusion tensors $\mathit{D}={\sum }_{\alpha }^N{\mathit{D}}^{\alpha }{\varphi }_{\alpha }$ weighted by the local volume fractions ${\varphi }_{\alpha }$. The discretization is based on the finite difference scheme described in more detail in Appendix A.

2.3. Simulation Parameters for High-Nickel NMC

Lithium layered-oxides with high nickel content, such as NMC811, deliver high capacities but also face stronger capacity fading due to mechanical degradation and surface densification [4]. The cycling of LiNiO${}_2$ is accompanied by multiple phase transitions [28,32] and indications for similar transformations have been observed for NMC811 as well [4], although the voltage curve at low (dis)charge rates does not feature distinct steps and plateaus (see Figure 1). For the purpose of this study, we assume an ideal solution behavior in the voltage range 3.6–4.2 V. As a consequence, the chemical free energy and the resulting chemical potential are expressed as

$$\begin{array}{c}\hfill f_{\mathrm{NMC}}={\mu }_{\mathrm{NMC}}^{\ominus }c+Kln\left(c\right)+K(1-c)ln(1-c),\end{array}$$

(5)

$$\begin{array}{c}\hfill {\mu }_{\mathrm{NMC}}=\frac{\partial f_{\mathrm{NMC}}}{\partial c}={\mu }_{\mathrm{NMC}}^{\ominus }+Kln\left(\frac{c}{1-c}\right).\end{array}$$

(6)

The measurable open-circuit voltage (OCV) is related to the difference in chemical potentials between the cathode and the anode [33], which in this case yields

$$\begin{array}{c}\hfill V_{\mathrm{OCV}}=-\frac{{\mu }_{\mathrm{Cathode}}-{\mu }_{\mathrm{Anode}}}{e{N}_{\mathrm{A}}}=-\frac{{\mu }_{\mathrm{NMC}}-{\mu }_{\mathrm{Li}}^{\ominus }}{e{N}_{\mathrm{A}}}=V^{\ominus }-\frac{K}{e{N}_{\mathrm{A}}}ln\left(\frac{c}{1-c}\right)\end{array}$$

(7)

where the chemical potential defined in Equation (6) is used for the NMC cathode, and we assume a lithium metal counter electrode with constant chemical potential ${\mu }_{\mathrm{Li}}^{\ominus }$. The reference voltage $V^{\ominus }=({\mu }_{\mathrm{Li}}^{\ominus }-{\mu }_{\mathrm{NMC}}^{\ominus })/\left(e{N}_{\mathrm{A}}\right)$ is set to $3.89$ V, which results in the blue curve in Figure 1 for the equilibrium electro-chemical potential of the cathode with respect to a lithium metal anode. e denotes the elementary charge and $N_{\mathrm{A}}$ the Avogadro constant. The ion mobility is assumed to increase with concentration due to more random hops, while the crowding at high concentrations leads to a decrease such that $M\propto c(1-c)$. The thermodynamic factor $\partial \mu /\partial c^{\alpha }$ calculated based on Equation (6) yields $K{c}^{-1}{(1-c)}^{-1}$. Inserting the mobility and the thermodynamic factor derived from the ideal solution into Equation (3) leads to a chemical diffusivity that is independent of c.

Lithium diffusion is confined to the plane between transition metal oxide sheets in a defect-free structure, and thus, is highly anisotropic. Reported values of the apparent diffusion coefficient in NMC differ by orders of magnitude [4,34], which stems from the multiple challenges involved in aquiring reliable data. Computation of the apparent diffusion from galvanostatic intermittent titration technique (GITT) relies on the assumption of semi-infinite relaxation that is hard to ensure for small primary crystals and high diffusivities. The fracture of agglomerates and subsequent penetration of electrolytes lead to an increase in surface, and consequently, an overestimation of diffusivities obtained from GITT measurements [34]. Within this study, the in-plane diffusion coefficient is set to $D_{\left|\right|}={10}^{-10}$ cm${}^2$/s, which is one order of magnitude faster than the apparent diffusion measured in [34]. In the direction perpendicular to the planes, diffusivity is assumed to be negligible ($D_{\perp }=0$). At grain boundaries, a mixture rule is applied to interpolate diffusion tensors of neighboring grains. To the best of our knowledge, the influence of grain boundary diffusion on ion transport within polycrystals is not sufficiently understood. Due to a lack of reliable data, no enhanced grain boundary diffusion is modeled although it is technically straightforward in the framework of the multiphase-field method [35].

3. Results

Lithium diffusion in NMC811 is simulated in two-dimensional domains for varying microstructures, as shown in Figure 2. All secondary particles in this study are assumed to be spherical with an outer diameter of 10 µm [4] (i.e., $R=5$ µm). Inspired by various works with textured, rod-like microstructures [7,8,36], we assume a compact core with randomly oriented grains surrounded by a layer of radially aligned grains.

The radius r of the inner core is systematically varied, starting from $r/R=1$ that corresponds to a sphere of randomly oriented grains (“gravel-NMC” [8]), to a ratio of $r/R=1/4$ where the microstructure is dominated by aligned grains. Randomized grain structures for the inner core are created by Voronoi tesselation. All structures are pre-conditioned using the multiphase-field method, in the sense that diffuse transition regions between order parameters of indiviual grains are established by solving the corresponding phase evolution equations [11,26,27]. The grain structures are assumed to be stable at room temperature such that no microstructure evolution occurs during cycling.

All structures are subjected to the same electro-chemical simulation setup where we apply voltage steps of $\Delta V=0.1$ V, followed by a relaxation time of 1 h. The voltage steps are applied via a Dirichlet boundary condition for c at the surface of the particles employing the relation between the voltage and c in Figure 1. The characteristic relaxation time of the model setup is approximately $\tau =R^2/D_{\left|\right|}=0.7$ h. Due to the different microstructures, diffusional relaxation within the secondary particles varies strongly, as can be seen in Figure 3.

While 1 h of relaxation time is sufficiently long for the strongly textured samples ($r/R=1/3$ and $r/R=1/4$) to reach an equilibrium state ($\Delta c=0$), the other samples are out of equilibrium. Figure 3a,d clearly show that the gravel-NMC sample builds up larger concentration gradients over time as a result of slower ion diffusion. This effect is also seen in the shifting of currents for the $r/R=1$ sample in Figure 4a. The spatial ion distribution is more inhomogeneous as some grains block the direct diffusion path towards the center of the sphere, resulting in a more tortuous pathway. The strongly aligned radial grains, moreover, provide fast ion transport paths towards the inner of the particle. In both cases of aligned grains (Figure 3b,c), we observe an accumulation of ions at the grain boundary that is a result of the tapered shape of these grains. All ions moving through transition metal oxide sheets that end in a grain boundary need to diffuse further via the grain boundary. The particle with a ratio of $r/R=1/2$ exhibits inhomogeneities in the ion diffusion that appear as rays in Figure 3b. The effective ion transport through the radially aligned grains is influenced by the microstructure of the inner core. Strong misalignment at the transition between core and outer ring leads to a blocking of transport routes and as a result some radial grains exhibit stronger concentration gradients. This effect becomes negligible when the textured grains dominate the overall morphology of the particle (Figure 3c).

The composition gradients can, furthermore, be evaluated in terms of a virtual potentiostatic intermittent titration technique (PITT) experiment. The non-dimensional current resulting from the applied voltage step is calculated from the change in average concentration in the particle

$$I=\frac{1}{\int {\varphi }_{\mathrm{NMC}}\mathrm{d}V}\int \frac{\partial c}{\partial t}{\varphi }_{\mathrm{NMC}}\mathrm{d}V$$

(8)

where ${\varphi }_{\mathrm{NMC}}$ denotes the local volume fraction of all NMC grains, which in this case is identical to $1-{\varphi }_{\mathrm{electrolyte}}$. Following a transient regime, the current follows an exponential decay that can be identified as linear regions in the logarithmic plot Figure 4a. The results in Figure 4a underline that the microstructure influences the global current (i.e., influx of ions in this case). The linear regions of the various samples exhibit different slopes, which hints at varying effective diffusion coefficients. Employing a standard equation to compute the apparent diffusion coefficient from PITT experiments (Equation (2) [37])

$$D_{\mathrm{apparent}}=-\frac{\mathrm{d}ln\left(I\right)}{\mathrm{d}t}\frac{16R^2}{{\pi }^2}$$

(9)

yields the results shown in Figure 4b. Only regions where the current follows an exponential decay fulfill the assumptions that have been used to derive Equation (9). Over the course of ten potential steps, the computed apparent diffusivities are very reliable and the average values have been marked by dashed lines. It should be noted that the noise in our evaluation becomes larger the faster the average diffusivity is. This is an artefact of the small changes of ion concentration close to equilibrium, where the computed diffusivity values are dominated by trunctuation errors.

4. Discussion

The results highlight how crucially ion transport is affected by the microstructure of secondary particles for anisotropic materials. The spatial distribution of composition that is observed in the simulations (Figure 5b,d) closely resembles the experimentally obtained ion distribution in the work of Xu et al. [8]. Especially, the ray-like ion distribution in Figure 3a [8] can now be explained by the interplay of misaligned grains in the randomly distributed core with respect to the outer aligned grains.

The evaluation of apparent diffusion coefficients sheds some light on the vast range of values that can be found in the literature. In this study, the obtained values for $D_{\mathrm{apparent}}$ vary by a factor of 4, even though the input values for $D_{\left|\right|}$ and $D_{\perp }$ are identical. It should be noted that the apparent diffusion coefficients calculated from Equation (9) are too high in this case. Although $D_{\left|\right|}=1\times {10}^{-10}$ cm${}^2$/s was chosen, the apparent diffusion coefficient for rod-NMC is $D_{\mathrm{apparent}}^{\mathrm{rod-NMC}}=8\times {10}^{-10}$ cm${}^2$/s, which is physically not possible. This results from the fact that Equation (9) has been derived for spherical particles while the simulations are only a two-dimensional representation of spherical microstructures. This does not change the qualitative differences between the obtained values as all could be corrected by a geometry factor. However, this discussion underlines the importance to question the assumptions that have been used for the derivation of formulas that are typically used for the evaluation of GITT or PITT experiments. The detailed investigation unravels the influence of the microstructure on the effective transport properties, which is only one of many parameters influencing the experimental results. It should be acknowledged that any experimental measurement of diffusivity can only capture the apparent diffusion in the specific microstructure, and thus, diffusivity values should always be published in combination with a characterization of particle morphology.

Furthermore, the simulation results suggest that grain boundary diffusion could play an important role for the ion transport in dense secondary particles. In the randomly oriented sample, an enhanced grain boundary diffusion would mediate the blocking effect observed at interfaces between strongly misaligned grains. Furthermore, in the rod-NMC with a strong texture, grain boundaries play a vital role in the overall transport, as sketched in Figure 6.

Due to the tapered shape of radially aligned grains, many layers are cut by grain boundaries, and thus, the migrating ions need to diffuse further into the particle via grain boundary diffusion. An enhanced grain boundary diffusion would lead to the faster transport in the whole particle, while sluggish grain boundary diffusion would result in many trapped ions and strong concentration gradients at grain boundaries. To the best of our knowledge, the actual influence of grain boundaries on ion transport is not well understood.

5. Conclusions

The presented framework is well-suited to study the influence of particle morphology in terms of size, distribution, and alignment of primary particles on key battery performance measures such as the effective diffusivity. Through a virtual PITT test, we are able to compute an effective diffusivity depending on the actual microstructure.

Within this work, the phase-field method is employed to parametrize the polycrystalline structure of NMC agglomerates. While this study is based on virtually generated structures, the regular grid employed to solve the partial differential equations also enables the straightforward inclusion of image data. Furthermore, the hierarchical particle morphologies investigated in this work are relevant for other layered-oxide cathode materials, such as the P2-type ${\mathrm{Na}}_X{\mathrm{Ni}}_{1/4}{\mathrm{Mn}}_{3/4}{\mathrm{O}}_2$, which is a promising material for sodium intercalation batteries [9]. This material exhibits multiple plateaus in the OCV due to phase transitions. The multiphase-field framework employed throughout this work naturally allows for the computation of phase evolution within the grains such that future investigations will encompass the coupled effects of phase transformations, ion diffusion, and mechanical degradation. A systematic study of electrolyte infiltration after crack formation and the resulting effective diffusion coefficients becomes possible. By explicitly modeling the morphology dependence of ion diffusion, this work contributes to bridging the gap between diffusivity values obtained at different scales, e.g., from first principle methods such as DFT [38,39], as well as from experimental methods such as GITT [4,34] and PITT [37]. The proposed framework enables the systematic study of the effective transport properties of secondary particles, which is relevant for Newman-type models in which secondary particles are typically homogenized as isotropic spheres.