Image based modelling of microstructural heterogeneity in LiFePO4 electrodes for Li-ion batteries

3D structure of LiFePO4 electrode has been imaged using X-ray nano-CT for the first time. CFD simulation has been used to compute tortuosity which is compared with the Bruggeman correlation. Local heterogeneity in the microstructure is described using vectorial tortuosity. Correlation between transport and geometrical tortuosity established for this electrode. a r t i c l e i n f o

h i g h l i g h t s g r a p h i c a l a b s t r a c t 3D structure of LiFePO 4 electrode has been imaged using X-ray nano-CT for the first time. CFD simulation has been used to compute tortuosity which is compared with the Bruggeman correlation. Local heterogeneity in the microstructure is described using vectorial tortuosity. Correlation between transport and geometrical tortuosity established for this electrode.

Introduction
Many simulation techniques treat Li battery electrodes as macro-homogeneous and isotropic, using a few scalar parameters, such as porosity and conductivity, to represent the internal structure [1]. Recent advances in tomography techniques have allowed the microstructure of nano-scale systems, such as battery and fuel cell electrodes [2] to be captured in 3D and, therefore, used to establish the validity of these assumptions for working battery electrode systems. These experimental developments are summarised in a number of recent review papers [2e4] and have found increasing application to the exploration of the complex porous media characteristic of Li battery electrodes.
The most widely used 3D imaging techniques for applications in Li battery research are FIBeSEM tomography [5e7] and X-ray micro-and nano-computed tomography [8e12]. Of these techniques, we have focussed on the application of laboratory and synchrotron X-ray techniques, which demonstrate multi-length scale capabilities [8,13] and afford the non-destructive capability to explore processes of microstructural evolution over time.
LiFePO 4 (LFP) has attracted significant attention in recent years as a safe and economical alternative to LiCoO 2 . First proposed by Huang and Goodenough [11], the bulk material does not exhibit sufficient electrical conductivity to function as a viable cathode material. However, recent nano-structuring of LiFePO 4 has provided conductivity improvements and increased its viability for commercial development. It is therefore by nature, a challenging material to image, owing to characteristically small feature size commensurate with the required conductivity improvements. The high resolution (typically <100 nm) required to accurately capture this structure has meant that, to date, 3D imaging has been conducted using FIBeSEM techniques, where sufficient resolution and contrast have been achieved to discriminate between the active material, binder and pore phases. Here we present, to the authors' knowledge, the first study of the 3D microstructure of a LiFePO 4 electrode conducted using synchrotron X-ray microscopy. We have subsequently utilised the 3D microstructure as a basis for image based modelling of transport within the complex porous structure, providing an insight into the microscopic properties of these devices as well as a tool through which we can establish and explore the effects of local microstructural heterogeneity.
In this paper, we first, briefly, introduce the necessary theory associated with diffusion in complex porous media before applying computational fluid dynamics (CFD) tools to characterise the tortuosity of these working electrodes, based on the analysis of a commercially available electrode.

Background theory
In complex porous media, two phenomena can be considered to govern the rate of diffusion. The continuum description provided by Fick's law, assumes that the primary mechanism of diffusion is associated with moleculeemolecule collisions as species are transported down to a concentration gradient. With decreasing pore size, the likelihood of moleculeewall collisions is increased e this is indicated by the Knudsen number, which is defined as the ratio of a characteristic pore length scale to the mean free path of the diffusing species. As the Knudsen number (Kn) approaches unity, the continuum Fickian description becomes limited in its ability to accurately describe the diffusion phenomena and we must instead consider a statistical description of the diffusion process [12].
Ionic diffusion in Li-ion batteries is typically considered to behave according to Fick's law with negligible contributions from Knudsen effects e in this case, the form of Fick's law most widely used, accounting for the porosity and tortuosity of the porous media, is given by Eq. (1), which relates the flux (F) to the diffusivity (D) and the concentration gradient ðV4Þ.
In this equation, both porosity, ε, and tortuosity, s, are dimensionless parameters that account for obstruction to diffusion by porous networks e geometrically, the tortuosity is considered as the ratio of the distance a molecule travels through the twisted pore network (Dl) to the overall thickness (Dx). A two-dimensional simplification is presented in Fig. 1; however, translating this definition to a three-dimensional (3D) structure is non-trivial. This definition of tortuosity can be extended to define an effective diffusion parameter, D eff , via Eq. (2), where D is the intrinsic diffusivity, usually found empirically [14].
In this form, the tortuosity is used to quantify the resistance to diffusion caused by the convolution of the pore network. In general, the values used for this purpose are taken from one of many porosityetortuosity correlations, rather than being measured from the specific sample being considered; the most widely used of which is the Bruggeman correlation [5,14e16]: The value of m that is usually assumed for battery models is 1.5 [14]; however, it is important to note that Bruggeman only used this value to describe the tortuosity of a specific structure containing dispersed spherical obstructions to diffusion [17,18]. Clearly, this single class of structure cannot accurately describe the diverse reality of porous electrode structures. A detailed discussion of this is beyond the scope of this paper (see e.g. Shen et al. [19]); however, there are two requirements of tortuosity which must be satisfied for any arbitrary structure [20]. Firstly, tortuosity must always be greater than or equal to one, as otherwise it would represent an augmentation of the diffusivity, Eq. (4). Secondly, as porosity tends to 100%, the tortuosity tends to one, because the properties of the system must be entirely described by the bulk properties of the material Eq. (5).
It is also worth noting that based on the definition of tortuosity in Eq. (1), for a system with ε ! 0 that is unconnected (i.e. does not conduct at all), s ¼ N.
The analogy between mass, heat, and momentum transport is well established [21]. Eq. (6) is a modified form of the heat transfer equation (analogous to Eq. (1)) where Q pore is the heat transfer; k is the thermal conductivity; ε is the porosity; A is the area of the face normal to the direction of heat transfer; DT is the imposed temperature difference between two parallel faces; L is the length of the sample in the direction of heat transfer; and s is the tortuosity of the conduction path. The ratio of ε to s is called the "obstruction factor" [22] and it represents the resistance to heat conduction of an otherwise uniform and fully dense material.
Modelling the transport of mass through porous media is a commonly encountered problem in engineering and, as shown by Zalc et al. [23], the tortuosity of a pore network is independent of the diffusion mechanism. Whilst Knudsen diffusion is not considered relevant in the case of ionic transport in Li-ion batteries, in order to maintain generality of the model, a heat transfer analogue has been used to calculate the pore network tortuosity here. Within the heat transfer analogy, the dimensions of the system (i.e. the sample volume) are arbitrary and so the continuous nature of the hypothetical heat conducting medium at any scale can be assumed.
The heat transfer analogy treats the pores as a solid material with a thermal conductivity of unity and adiabatic surfaces. The temperatures at two parallel faces of the control volume considered are held constant with a known temperature difference. The remaining external faces of the control volume are treated as adiabatic surfaces.
The thermal conductivity, k, need not bear any relation to the actual material being studied and was set to unity for all models described for simplicity of interpretation. This is made clear in Eq. (8) where the k terms have cancelled. The cross-sectional area of the control volume, A, is calculated for a face normal to the direction of heat transfer and the conduction length, L, is the thickness in the direction of heat transfer; it is worth noting that all control volumes described in this report are cubic. The imposed temperature difference, DT, is arbitrary and for the sake of convenience was set to 100 K. The porosity, ε, is defined as the volume of the pore network divided by the control volume. It is important to note that the value of the porosity includes those pores that are unconnected to the system and, therefore, cannot conduct. Porosity is often measured experimentally by high pressure fluid intrusion, which cannot detect unconnected pores [24,25]; however, by utilising the direct imaging method of XRM, these pores are detected. The tortuosity, s, quantifies the degree of resistance to heat transfer caused by the convolution of the pore network and, by this definition, accounts for the apparent loss in conduction caused by unconnected pore volume.
In order to calculate the tortuosity, the heat transfer measured through the pore network must be compared to that conducted through a uniform, fully dense sample with the same dimensions as the control volume. This can be found using Eq. (7), where Q CV is the heat transfer through the control volume and the other variables have the same meaning as Eq. (6).
Dividing Eq. (6) by Eq. (7) yields Eq. (8), which can be rearranged to the form shown is Eq. (9) and used to calculate the tortuosity from the porosity and the heat transfer rates.
The definition of tortuosity employed in this work is consistent with that adopted by Kehrwald et al. [18] as well as many others [23,26e28]. However, care must be taken to avoid confusion with the definition used by Shen and Chen [19] and others [29e31], who use s 2 in place of s. It should also be noted that the final equation used for the determination of tortuosity, Eq. (9), does not contain parameters relating to the heat transfer analogy, such as the temperature difference or thermal conductivity other  than the ratio of heat transfer rates, which is in line with expectations.

X-ray characterisation
A commercially available LiFePO 4 cathode, manufactured by Hydro Quebec, was used as the basis for the study. The electrode was obtained in a fresh state (i.e. before incorporation into a battery) and a small section of the electrode sheet was cut using a sharp blade, to form a point which matched the field of view of the X-ray microscope (see Refs. [32,33]).
The sample was characterised using the Xradia nanoXCT-S100 system installed on beamline 6-2 at the Stanford Synchrotron Radiation Lightsource [34]. The Xradia XMController software was used to collect 721 projection radiographs over a 180 rotation angle range, using an individual exposure time of 1 s per radiograph. A photon energy of 7120 eV was used for the exposures and images were collected on a 1024 Â 1024 pixel CCD detector with an effective pixel size of 20 nm. After image acquisition, the resulting rotational series was reconstructed using a filtered back projection algorithm integrated into the Xradia XMReconstructor software package, which subsequently produced a stack of 907 virtual slices with 20 nm cubic voxels. Imaging details are summarised in Table 1.
The projection images were reconstructed using a standard parallel beam back projection algorithm (Xradia Reconstructor) to obtain the grey-scale image stack, where each image is a 2D slice through the sample. Analysis and binarisation of the images was conducted using Avizo (VSG Bordeaux), where an Edge Preserving Smoothing Filter [35] was applied followed by grey-scale threshold segmentation. The resulting binarised reconstruction, Fig. 2(a), was of non-cubic geometry and for ease of manipulation was cropped to the largest possible cuboid volume of 900 Â 600 Â 600 voxels corresponding to a total volume of w3000 mm 3 . Finally, a cubic region 440 voxels across, Fig. 2(b), was selected for use in the simulations; although this region constitutes just 9% of the original image set by volume, it represents the largest reconstruction of an LFP electrode known to the authors (ca. 80 vol% larger than Ender et al. [5] at comparable voxel resolution).

CFD simulation methodology
The reconstructed, segmented tomography data was converted from an image stack into a surface mesh (ASCII *.stl) file using Avizo Fire and subsequently imported into the STAR-CCMþ Ò CFD package (CD-adapco, London). The surfaces generated in Avizo often required some manual repair operations to be performed in STAR, such as patching holes, to guarantee that the mesh was closed and manifold, as a polyhedral mesh can only be generated in a well defined volume. The repaired regions generally represent less than 0.001% of the total volume and care was taken to ensure that the repairs would not significantly affect the porosity of the sample. Once repaired, a cube was created which determined the control volume to be investigated and then a Boolean subtraction operation was used to enclose the pore space as opposed to the solid volume.
Before generating the volume mesh, it was necessary to use the built-in remesher tool to re-tessellate the surface, as the *.stl generation in Avizo is designed for dimensional fidelity rather than cell quality, as illustrated by Fig. 3(b). The remesher tool ensures that the aspect ratio of all the surface mesh triangles is below a specified value, which leads to an improved quality of the final volume mesh. A polyhedral volume mesh was then generated using the surface mesh as a template.
A two stage mesh refinement study was conducted in order to ensure that computationally efficient and grid independent solutions were obtained. This resulted in a 13% reduction in the number of computational cells required, which significantly reduced the processing time required for the heat transfer simulation.
Having established optimal mesh conditions, the various surfaces of the volume were labelled in order to allow the appropriate boundary conditions to be specified. To measure the heat transfer  in, for example, the x-direction, the two faces parallel to the yez plane would be assigned fixed temperatures with an arbitrary temperature difference (in this case 100 K) and all other faces were considered adiabatic. The value of thermal conductivity set for the material was also arbitrary and, as such, was consistently set as unity. Only the steady state solution was required and so the physical continuum was defined accordingly. The heat transfer through the pore network, Q pore , was calculated at the high temperature plane for convenience, although this location is arbitrary as the heat transfer will be equal through any section normal to the direction of flow at steady state. The simulation was run until the residual reached an acceptably low value (below 10 À5 ), where the solution was considered to have converged. In the context of a heat transfer simulation, the residual refers to the degree to which the conservation of energy equation is satisfied for an iteration. This process was then repeated to calculate the heat transfer in the y-and z-directions. For each of the three measured directions, a plot could be generated, Fig. 4, which showed the temperature distribution within the pore network, useful for checking that the simulation was performing in line with expectations.

Results & discussion
Compared to its nano-structured average particle size, the sample of LiFePO 4 was relatively large and, as such, facilitated a detailed investigation into the validity of the isotropy and heterogeneity assumptions implicit in many battery models. For the largest bounded cubic region of the sample, which constituted around 9% of the volume of the original data set, the overall tortuosity and porosity is given in Table 2.  The porosity can be calculated for each slice through the sample, where a slice is a plane of single voxel thickness. The graphs in Fig. 5 show this planar porosity of the sample in the x-, y-and z-directions, where the slices are normal to the axis along which they are being compared.
Following the formalism set out by Kehrwald et al. [18], to investigate the degree of anisotropy of the LFP sample, the volume was subdivided into 8 cubic, non-overlapping regions as illustrated by Fig. 6(a). Table 2 contains the porosity and the directional tortuosities of these 8 control volumes.
The tortuosities of the eight cubes were re-combined using a series/parallel resistor network analogue. The example resistor network in Fig. 6(b) would be suitable for calculating the tortuosity in the x-direction as defined in Fig. 6(a). Comparing the values from the single large LFP control volume to the combined values of the 8 smaller cubes yields an error of less than 1% for both the measured porosity and characteristic tortuosity. These small errors in the calculated results suggest that the combination method was valid for the highly percolating LFP sample, the small (<0.2 vol%) difference in porosity observed is attributed to mesh smoothing during sub-division of the main mesh.
In an attempt to compare the calculated directional tortuosities to the Bruggeman, a quantity called the characteristic tortuosity, s c , was defined in Eq. (10). Additionally, a measure of the anisotropy, s, was also defined, Eq. (11), as the ratio of the range of the directional tortuosity to s c . Table 2 shows the calculated values of the characteristic tortuosity, s c , and degree of anisotropy, s, for all regions.
Plotting the characteristic tortuosity of each sample against the Bruggeman correlation, as in Fig. 7(a), allows a clearer comparison to be made than with the range of directional tortuosities measured. However, it is concluded that the Bruggeman correlation remains a poor predictor for tortuosity.
Avizo Fire contains an algorithm for measuring a type of geometric tortuosity, which was also employed in this investigation. This method works by first calculating the pore-network centroid within each slice, and then tracking its movement between each slice. The tortuosity is subsequently calculated by dividing this centroid motion path length by the distance from the first slice to the last [36]. Table 3 gives the values measured in each direction for this geometrical tortuosity, s geo . This table also contains values for the volume specific surface area, a, which was calculated using the "area3d" function in Avizo Fire. Fig. 7(b) is a graph comparing the results for geometric tortuosity measured in Avizo against the values obtained in STAR using the heat transfer analogy. The two methods appear to be quite well correlated and a curve defined by Eq. (12) fits the points with an R-Squared value of 0.91. This curve also has two major advantages over a linear fit. Firstly, it goes through the point (1,1), which is in line with the definition of tortuosity described in Section 2 of this work, and, secondly, the maximum possible values of s CFD are expected to be much higher than those of s geo , which is well represented by a logarithmic fit.

Concluding remarks
Here the steady-state heat transfer analogy has been used to quantify the effect of tortuosity of porous electrode microstructures on the resistance to diffusion through a material. The microstructural data utilized in this investigation was obtained for an LiFePO 4 battery electrode using synchrotron nano-scale X-ray microscopy.
The tortuosities measured by this CFD simulation did not correspond well with the widely used "Bruggeman correlation", which is perhaps unsurprising as the latter is based on the assumption that tortuosity is isotropic and homogeneous throughout a material. In order to compare the highly anisotropic pore networks observed to  the structures described by the Bruggeman correlation, a term called the characteristic tortuosity was defined, which combined the three directional tortuosities into a single value. A second, simplified, measure of tortuosity based purely on the pore geometry was also obtained from the CT volume using a sliceby-slice pore centroid calculation implemented in Avizo Fire. A logarithmic relationship between the geometrical and transport tortuosities was established with reasonable confidence. Notably, the relationship between the geometrical tortuosity and the Bruggeman prediction showed a marked correlation, which is a rather unexpected result as the Bruggeman relationship is derived from a transport (not geometrical perspective). It is possible that for the porous samples studied (ε ¼ 34%e48%) the geometrical and Bruggemann tortuosities fall within similarly limited possible ranges and further work is underway to examine this.
Computational models of batteries and fuel cells tend to treat tortuosity as a single scalar parameter characterising an entire electrode. The simulations undertaken in this report demonstrate that tortuosity is not a simple scalar quantity, but instead both geometrical and transport tortuosities show a marked dependence on direction. This suggests a directional (vectorial) representation of tortuosity should be developed and an example methodology for doing this is proposed.
Finally, all simulations were run under steady-state conditions; as such, it can be conceived that pore structures with highly different geometries can exhibit the same value of tortuosity e for example the contributions to mass transport resistance from pore twistedness and pore constriction cannot be decoupled. Whilst this is not important for steady state scenarios, for electrochemical mass transport limiting scenarios, where transient mass transport must be considered, it may become important to include these characteristics. Under such conditions it is suggested that a complex analysis of mass transfer impedance (mass resistance and mass capacitance) must be considered in order to successfully decouple the effects of pore constriction and pore twistedness that will impact the transient mass transfer behaviour.