Bruggeman’s Exponents for Effective Thermal Conductivity of Lithium-Ion Battery Electrodes

Simulations of lithium-ion battery cells are usually performed with volume averaging methods that employ effective transport properties. Bruggeman’s model, which is widely used to determine these effective properties, is primarily based on the volume fraction of porous electrodes. It does not consider actual particle shape, size or the topology of constituent phases; these play a crucial roleindeterminingeffectivetransport.Inthispaper,ageneralderivationoftheeffectivethermalconductivityofmultiphasematerials,whichcanbecorrelatedwiththesefactors,isderivedusingthevolumeaveragingtechnique.Three-dimensionalﬁnitevolumemeshesoffully-resolvedlithium-ionbatterycathodemicrostructuresarereconstructedfromscannedimages.Effectivevolumeaveragedthermalconductivityisthendeterminedfromnumericalanalysisofthermaltransportonthesemeshes.ItisshownthattheBruggemanmodelforeffectivethermalconductivitymustberecalibratedtoﬁttheeffectivethermalconductivitycomputedfromthesedetailedsimulations.Therelevanceoftheseresultstoeffectivetransportpropertiestypicallyemployedinelectrochemicalsimulationsispresented.Commonlyusedtheoriesforeffectivethermaltransportincompositesareevaluatedforcomparison.Furthermore,itisshownthatBruggeman’sexponentsyieldanimportantquantitativemeasure,theconnectivity,tocharacterizethephysicalpathfortransportthroughtheunderlyingphases.©TheAuthor(s)2015.PublishedbyECS.ThisisanopenaccessarticledistributedunderthetermsoftheCreativeCommonsAttribution4.0License(CCBY,http://creativecommons.org/licenses/by/4.0/),whichpermitsunrestrictedreuseoftheworkinanymedium,providedtheoriginalworkisproperlycited.[DOI:10.1149/2.0151602jes]Allrightsreserved.

The importance of lithium-ion batteries is now well recognized in light of the global energy crisis, global warming and the need for efficient and inexpensive energy storage options. 1,2 Battery physics encompass thermodynamics, electrochemistry, material science, transport phenomena and solid mechanics, and span multiple length and time scales. 3 Realistic modeling of batteries across these disparate physics and scales is critical for their effective and safe commercialization. 4 Research in lithium-ion batteries has been primarily driven by the need to develop cathode, anode and electrolyte materials that deliver high potential and capacity. 5,6 On the cell level, the battery consists of a porous composite anode and cathode filled with an electrolyte and separated by a separator. Typically, the thickness of such a sandwich is hundreds of microns, while lateral dimensions are of the order of centimeters. Lithium ions shuttle between cathode and anode during charge and discharge. Active cathode particles (for e.g. LiCoO 2 , LiMn 2 O 4 , LiFePO 4 ), which form the key constituent of batteries, are lithium insertion compounds and have high lithium ion conductivity. These active particles are spatially dispersed and held together by binder (such as polyvinylidene fluoride [PVDF]), while the electrolyte resides in the pores, wetting the active particles, and facilitating transport of lithium ions. The electronic conductivity of these active materials is generally low compared to lithium ion conductivity. Therefore additives like carbon particles are dispersed in the cathode to promote electronic conductivity, albeit at the cost of reduced energy density. On the cell level, experimentalists have tried to optimize the performance of batteries by changing the dimensions of active particles, volume fraction of the constituents and pores, etc. 7,8 But only a few computational and simulation studies exist to guide these experiments.
The most widely used models for simulating electrochemical and thermal behavior of batteries on a cell level are volume averaged models. 9,10 Volume averaging involves homogenization of the constituent phases of the porous electrode, eliminating the need to describe and represent the micro-structural configuration. 11 The elimination of microstructural information is reflected in the need for modeling the effective properties of the volume-averaged medium. These models must, in principle, take into account the intrinsic properties of constituent phases, their volume fraction, and the structure of the battery electrodes, among other factors.
However, until recently, it was not possible to obtain real 3D microstructure of porous electrodes. This led to the numerous simplifications in the process of volume averaging. For instance, these models assume that the active particles are spherical in shape or are dispersed far apart from each other. Thus the resulting effective transport properties take into account only the volume fraction and the intrinsic properties of constituents. In reality, the topology and geometrical shape of the constituents significantly influence the effective transport properties. Topology, loosely defined here as the 3D network formed by the constituents, defines the pathways for ionic, electronic and thermal transport. It is therefore possible to have varying effective properties for the same volume fraction because of the differences in topology. Similar effects can be attributed to the deviation of the geometry of particles from canonical shapes like spheres or ellipsoids. 12 At the same time, topology and geometry cannot be considered static for batteries. They change during the process of charging/discharging because of particle volume changes or due the fracture/cracks developed during cycling as a result of differential stresses. 13 It is therefore necessary to develop models based on the actual three-dimensional geometry, fully resolving the various phases (active material, carbon, binder, and electrolyte) of the electrode. This is being made possible with the advent of instruments like micro-computed tomography (μCT ) scanners and focused ion beam-scanning electron microscopes (FIB-SEM) in conjunction with advanced image processing techniques. 14 Unlike volume averaging, performing transport simulations on geometries explicitly partitioned into various constituent phases can be termed as fully resolved pore-scale or particle-scale simulations (some authors 4 also use the term direct numerical simulation (DNS)). A few papers have been published using pore-scale simulations, employing simple synthetic geometries and other simplifications. [15][16][17] Additionally, pore-scale electrochemical simulations on real particle bed reconstructions have recently been reported. 18,19 However, the length scales of fully-resolved reconstructed volumes attainable are orders of magnitude smaller than the global length scale of a battery cell and therefore are not amenable to simulations directly on the global scale. A feasible pathway is to use a combination of pore scale simulations and volume averaging using a representative elemental volume (REV). An REV refers to a spatial volume that is large enough compared to the particle length scale but small compared to the dimensions of the cell. Because of periodic nature of the porous structure, an ideal REV can be chosen that is able to capture the essential topological and geometrical information statistically. By performing a pore-scale simulation on these fully resolved REVs, it is possible to develop constitutive models for transport properties for use in volume averaged A120 Journal of The Electrochemical Society, 163 (2) A119-A130 (2016) transport models. Bruggeman's relation 20,21 has commonly been used to correlate and model effective transport in batteries.
Thermal issues in lithium-ion batteries remain one of the key bottlenecks to achieving efficiency, safety, and life. 22,23 One of the critical issues is thermal runaway, in which inadequate heat dissipation raises battery temperature above critical limits, and triggers exothermic reactions in a positive feedback loop. 22 Therefore it is important that effective thermal properties of battery materials like thermal conductivity and specific heat capacity be characterized accurately. Few experimental measurements of these properties have been reported in this regard. 24,25 Thus far, volume-averaged thermal simulations have been carried out using such data as are available. 26,10,27 In this paper, we develop the basic methodology for deriving effective thermal conductivity from detailed pore-scale simulations of heat transport in battery microstructures. We correlate our computed transport properties not only to the volume fractions of constituents but also to their microstructure geometry using modified Bruggeman's exponents. Furthermore, we address other effective transport properties such as mass diffusivity, ionic conductivity and electrical conductivity. We demonstrate that commonly used power law exponents in the Bruggeman model 28 must be recalibrated to fit the effective properties computed from these fully resolved pore-scale simulations. An important outcome of our simulations is the insight that these modified Bruggeman's exponents provide into the structure and connectivity of the underlying electrode material.
Section Theory and modeling describes the volume-averaging of the multiphase heat conduction equation from first principles. The procedure for obtaining the real 3D microstructure geometry and the corresponding finite volume mesh for numerical simulations is described in section Geometry reconstruction and mesh generation. Section Numerical simulation describes the process of performing pore-scale simulations on these reconstructed geometries. The results of these simulations, including the computation of modified Bruggemann's exponents, are discussed in section Results. The comparison of thermal results with other effective medium theories, quantification of the microstructure geometry with exponents, and the relevance of the exponents to electrochemical simulation are also presented in this section. We close the paper with a brief description of follow-up work.

Theory and Modeling
We present here a concise derivation of the determination of effective thermal conductivity of a two-phase material using the volumeaveraging technique proposed by Slattery. 29 Mass, momentum (flow) and heat transport equations for multiphase or porous materials have previously been derived based on this technique. [30][31][32]11 Intrinsic heat conduction equation.-The transient heat conduction equation for two phases in the Fourier limit is given by where the subscript i refers to i th material phase. Here J i is the heat flux, T i is the spatial distribution of temperature, k i , ρ i , C pi are the intrinsic thermal conductivity, density and specific heat capacity respectively, and ∇T i is the gradient of the temperature field for phase i. In addition, temperature and heat flux must be continuous at the interface S between the two phases.
Volume-averaged equations.- Figure 1 shows the schematic representation of a two-phase material, where phase 1 (henceforth referred as the particle phase) is dispersed in phase 2 (referred to as the matrix). The length scale of the particles is l p . Figure 1b shows a representative elemental volume (REV). In order to capture statistical properties adequately, 32 the length scale of the REV, l REV , must be larger than l p , but much smaller than that the global length scales of the domain under consideration, l d i.e. l p l REV l d . Assuming that ρ i , C pi and k i are spatially constant, one can apply the volume averaging theorems in Refs. 29,11 to the equations for the intrinsic media (Eq. 1) in the REV, and derive the volume averaged equations for phase i Here the local volume average temperatureT i is defined as The normal vector n i j points from phase i to phase j at the interface. The term Q i j represents the interfacial heat transfer from phase i to j through the interfacial area S inside the REV. The term i represents the thermal tortuosity vector and quantifies the deviations of the heat flux vector caused by the physical obstruction of the transport path in one phase due to the presence of the other phase. 11 The term k i ∇T i represents the global heat flux vector, while the tortuosity vector represents the effect of the local thermal gradient within the REV. According to Eq. 3, the tortuosity term i should be a positive quantity that decreases the term k i ∇T i . The effect of structure geometry enters the governing equations through the terms i and Q i j . Therefore volume-averaged models postulate approximate models for these terms to obviate the need for explicit microstructural information.
It is convenient to model the tortuosity term i as a fractional value of the global gradient term k i ∇T i 30,11 along with a geometrical factor H i so that It is common practice 11 to write the expression for the geometric factor H i in terms of a volume fraction ε i and an exponent d i given by Thus 0 ≤ H i ≤ 1.
Using the definition of intrinsic volume averaged temperature, T i = ε iTi , followed by the assumption that ε i is spatially and temporally invariant, the volume averaged equation for heat diffusion may be written as, In the absence of fast transients and heat generation terms, 32 the condition of local thermal equilibrium may be invoked i.e., ∇ T 1 = ∇ T 2 = ∇ T . Assuming steady state and adding the equations for the two constituent phases, we obtain Equation 11 implies an effective thermal conductivity for the combined medium given by Volume averaged charge and species transport equations.-The volume-averaged equations described above may be applied not only to heat transfer in the electrodes, but also to electrochemical simulations. In order to generalize the results in this paper, we consider a lithium-ion battery cathode as an example. The equations governing electrochemical transport in the cathode include those for the electrochemical potential in the electrolyte and the cathode particles, φ E and φ S , and those for the lithium ion concentration in the electrolyte, c E , and in the cathode particles, c S . The transport of electrons through the cathode particles due to the electric potential is described by The transport of lithium ions through the electrolyte across the cell due to their positive charge and the presence of an electric potential gradient is given by The lithium ions that enter the electrolyte from the cathode particles are transported across the cell through the electrolyte and their motion due to mass diffusion is modeled by Finally, the transport of lithium from the interior of a cathode particle to the surface where the electrolyte wets it is described by Volume averaging is usually not performed on the above equation (Eq. 16) since the diffusion length scale of the lithium ion is that of the particle length scale. This equation has typically been modeled in the literature through simplified single-particle models. 28,33,34 In the above equations, the subscripts S and E refer to the cathode particle phase and electrolyte phase respectively. Furthermore, c represents lithium ion concentration, c the volume averaged lithium ion concentration, D the mass diffusivity of lithium ions, φ the electric potential, σ S the electronic conductivity in the cathode particle phase, κ E the ionic conductivity in the electrolyte and κ D the diffusional conductivity taking into account the motion of the charged ions due to mass diffusion. The term −κ D ε d E E ∇(ln c E ) is typically negligible for typical discharge rates. Similar to Q i j in Eq. 10, T E S , T SE and R SE in Eqs. 13, 14 and 15 represent interactions at the particle-matrix interface. These equations are similar in form to the thermal equation, Eq. 10. One can define effective diffusivities and conductivities for the different phases as The pathways for transport of ions, electrons and heat are different in the cell. Consider an idealized battery where interconnected solid cathode particles are dispersed in an interconnected electrolyte phase (in reality binders are added to preserve the structure and carbon particles are added to enhance the electronic conductivity of cathode particles). After leaving the active cathode material, ions can only travel through the electrolyte phase (e.g. path 1 in Figure 2). Likewise, electrons travel only through the solid cathode phase (e.g. path 2 in Figure 2). However heat can travel though both phases, as depicted by path 3 in Figure 2. These differences result in significantly different Bruggeman's exponents, as will be discussed in section Relevance to exponents used in electrochemical model equations.
Bruggemann's effective conductivity.-Though there are many effective medium theories (discussed in section Comparison with other theories) for effective conductivity, we specifically discuss Bruggeman's theory here because of its extensive adoption in the lithium-ion battery community. Bruggeman extended Maxwell's formula using differential effective medium theory, to arrive at an expression for effective conductivity k e f f of a composite system composed of a conductive phase of volume fraction ε 2 having conductivity k 2 with a dispersed phase of particles occupying volume fraction ε 1 having conductivity k 1 given by, Bruggeman's derivation admits only particles of spherical shape and requires a wide range of sizes. Furthermore, the particles are assumed sufficiently far apart that their thermal fields do not interact, though there is no restriction on volume fraction of the constituents. For the particular case when the particles possess negligible conductivity (k 1 ∼ 0), Eq. 20 defaults to the following equation for effective conductivity, k ef f = k 2 (ε 2 ) 1.5 [21] It is this Bruggemann relationship (Eq. 21) that has often been used to correlate and predict effective properties in a porous structure, especially by the battery modeling community 28,35 for both electrochemical and thermal simulations. It is also noted here that Bruggeman's effective conductivity is zero (k ef f = 0), for the case when the conductivity of matrix phase is negligible (k 2 ∼ 0, k 1 = 0). The reason is because of the inherent assumption of discontinuity of the particulate phase for the non-interacting fields.
A number of alternative relations are available in the literature for the case of non-conducting dispersed particles. In rock physics, Archie's law 36,37 has been used to predict the electrical conductivity of fluid saturated in rocks. 38, 39 Archie's law is purely empirical and similar to Eq. 21 in form. The exponent in Archie's law has been found to vary between 1.3 and 4 in sedimentary rocks. It has been argued that the value of the exponent increases as the degree of connectivity of the corresponding conducting phase decreases. Sen et al. 40 developed a theory for conduction through the pore space in sedimentary rocks based on self-similarity principles. The theory addresses high particle volume fractions and ellipsoidal particles. They also assume that the particle fields do not interact, and also that the pore space is always continuous. Their theory yields an exponent of 1.5 for spherical particles, consistent with the Bruggemann theory. For non-spherical particles with axes perpendicular to the external gradient, the exponent is found to be greater than 1.5. For the case when the particle axis is parallel to the external gradient, the exponent is less than 1.5. De La Rue and Tobias 21 have found experimentally that for low-enough particle volume fractions (<0.4), Eq. 21 holds reasonably well for non-spherical suspensions, albeit with exponents somewhat different from 1.5.
The assumptions underlying the theories of Bruggeman and Sen et al. are not strictly satisfied for battery microstructures. Both the phases must be interconnected for the battery to function, and thus the particles are not sufficiently far apart to be non-interacting. In addition, particles are not always spherical in shape and the connectivity of the phases may evolve with cycling. Despite these violations, many published papers employ a value for the exponent in Eq. 21 of d = 1.5 for both electrochemical 28,34 and thermal simulations. 10 Thorat et al. 41 summarizes the results of relevant papers that use other exponents. Our objective is to obtain the exponents in Eq. 12 through porescale simulation of the heat conduction equation in fully-resolved microstructures.
Bruggeman's exponent and connectivity.-Bruggeman's exponents (d 1 or d 2 in Eq. 12) in the limiting case of either k 1 → 0 or k 2 → 0 may be interpreted as representing the degree of connectivity of the conducting phase. Following Glover in, 37 we consider a two-phase material in which one phase has zero conductivity, and the other, the i th phase, has non-zero conductivity. We define a geometry factor G i defined as: Defining the term connectivity as (which is the inverse of the tortuosity frequently used in the literature 37,42 ), we obtain, The value of connectivity χ i ranges from 0 to 1, with 1 representing maximum connectivity and 0 representing a discontinuous phase. When χ i = 1, Eq. 23 postulates that the area of cross-section available for transport is proportional to the volume fraction ε i . Fractional values of χ i represent additional resistance to transport resulting from local tortuosity, constriction and weak connectivity. That is, the model represents the effect of these imperfections as a decrease in the effective cross-sectional area for transport. A value of χ i = 0 represents a complete absence of percolative pathways, resulting in zero effective conductivity of phase i. Low values of χ i imply high values of d i , i.e., the lower the connectivity of the phase, the higher the exponents in Eqs. 12, 17, 18 or 19). When the non-conducting phase is composed of disconnected spheres of varying diameters embedded in a conducting medium, the exponent d i takes a value of 1.5 21,40 for the conducting phase. Thus, the connectivity χ i is a good indicator of the underlying microstructural pathways for transport. In section Connectivity and Bruggeman's exponents, we present χ i values for different microstructures to help understand the structure of the underlying cathode material.

Geometry Reconstruction and Mesh Generation
We turn next to the reconstruction of cathode geometry from X-ray tomography or dual beam FIB-SEM. Both techniques produce a series of 2D image slices which can be stitched together to obtain a full 3D image of the microstructure. X-ray tomography is a non-destructive technique and has been used for microstructure characterization of sintered porous materials, 43 metal foams 44 and lithium-ion batteries, 14 among others. Dual beam FIB-SEM is a destructive technique wherein the sample is milled using FIB and subsequently imaged using SEM. Alternate milling and scanning produces a series of 2D image slices. FIB-SEM is being widely used for microstructure characterization of solid oxide fuel cells (SOFC) and lithium-ion batteries (see, for example in Ref. [45][46][47].
For this work we reconstruct 3D geometries from images of battery microstructures obtained from various research groups 48,47 acquired using both techniques, the details of which are given in section Battery electrode microstructures. Figure 3a shows representative sample images of a battery electrode. The brighter regions correspond to the active material, and the darker regions to the non-active phases. The next step is to partition each individual image into sub-regions that represent the constituent phases. The process is called segmentation, which employs extensive image processing.
Segmentation of each individual image is performed by choosing a threshold value of brightness. However, as can be observed from Figure 3a, the transition of gray scale values is not crisp. Therefore image processing must be performed to eliminate spurious artifacts which would compromise the quality of the finite volume meshes ultimately generated. We employ the image processing module ScanIP of the software package SIMPLEWARE 49 for this purpose. Features of the software like image resampling, noise removal, island removal, floodfill filters etc. are applied on all the original scanned images. Parameters for applying these features can be optimized with some trial and error. We refer to our previous papers for detailed procedures. 43 The result is a set of sharpened images that can now be employed for region identification and segmentation based on a chosen threshold brightness value. Figure 3b shows some representative segmented 2D images for the battery electrodes.  The next step is to reconstruct a 3D surface/volume representation by stitching the series of 2D images thus acquired. This volume is converted to a finite volume mesh for numerical simulations using the ScanFE-Grid module of SIMPLEWARE. 49 To optimize the number of cells in the mesh, the volumetric and boundary mesh adaptation features described in Ref. 49 are used. The smoothed conformal mesh generated by the package is of very high quality and consists of approximately 10 7 cells. The mesh is then exported to our in-house finite volume simulation software package MEMOSA. 50 Figures 4c-4d show the representative reconstructed volumes of cathode phase, nonactive materials phase, the combined phases, and the corresponding finite volume mesh respectively for the battery electrode.
The intent is to treat the volume of reconstructed microstructure as an REV defined in section Volume-averaged equations. The dimensions of the volume to be reconstructed are chosen such that the assumptions underlying the use of an REV are satisfied as closely as possible. The details of dimensions of all the reconstructed volumes are given in section Battery electrode microstructures.

Numerical Simulation
The governing equations are the microscopic heat conduction equations for constituent phases given by the steady part of Eq. 1 with continuity of heat flux and temperature imposed at the interface between the phases. These governing equations are solved numerically on the REV computational domain. For thermal conductivity computation, the conjugate heat conduction problem is solved by imposing given-temperature boundary conditions on two opposing faces, with all other boundaries being held adiabatic. The effective thermal conductivity in a particular direction x i may be computed as, where q is the heat rate at the domain face, which has an area A normal to the direction of the gradient. The same procedure is applied to all three directions and the average of these three values is reported as the effective thermal conductivity.
We employ a conservative cell-centered finite volume scheme 51 for solving the governing equations (Eq. 1). The meshes are composed of arbitrary convex polyhedral cells. Temperatures and thermal conductivities are stored at cell centers. The governing equations are discretized by enforcing conservation on each cell as described in Ref. 52 which yields a system of linear equations for temperature at cell centers. A variety of linear solvers, including an algebraic multigrid solver, may be used for solving the resulting linear equations. 51 The governing equations are solved until convergence. Convergence criteria, based on the conditions of scaled residuals falling below a specified value, are employed for termination. 53 For the domain sizes considered in this work, meshes contain approximately 10 − 20 × 10 6 cells. Numerical solution of a conjugate conduction problem on these meshes consumes a maximum of 4 hours of computational time on one core of a 2.3 GHz 12-Core AMD Opteron 6176 processor. Table I summarizes the information on the battery microstructures that we procured from two research groups. Most of the results presented in this paper employ microstructural scans from the Wood group, ETH Zurich, 48 obtained using synchrotron radiation X-ray tomographic microscopy (SRXTM) (Sample A in Table I). The cathode material is a layered oxide LiNi 1/3 Mn 1/3 Co 1/3 O 2 (NMC cathode). Various batteries have been manufactured by changing the weight fraction of carbon black and PVDF binder (2%, 3%, 4% and 5% of each) and under different compression pressures (0 GPa, 0.03 GPa, 0.06 GPa and 0.2 GPa), resulting in 16 combinations of cathodes. A 3% 0.6GPa case is a battery made with 94% weight fraction of cathode material, 3% carbon and 3% binder, and then compressed with 0.6 GPa after the battery is cast and heated in the oven. The cathode particles used are spherical in shape (Figure 1a of Ref. 48) ranging in diameter from 5 to 15 microns approximately. From one of the reconstructed volumes in Figure 4, we observe that the cathode particles are generally spherical in shape. All the batteries are uncycled. As batteries are generally manufactured in a fully discharged state (maximum cathode lithiation), these cathodes can be assumed to be in maximum state of charge.

Battery electrode microstructures.-
SRXTM permits scanning of relatively large reconstructed volumes (700 μm × 700 μm × 70 μm for SAMPLE A, for example). 48 Four sub-regions (REVs) of size 55.5 μm × 55.5 μm × 55.5 μm from different random locations were cut from the image slices as depicted in Figure 5, and finite volume meshes generated for each of the  samples for numerical analysis. The procedure was repeated for each of the 16 microstructures, thus yielding a total of 64 REV samples.
The volume fractions of cathode and non-active regions as well as the corresponding porosity were determined from the geometry of the REV. As previously mentioned, the non-active phase is formed of carbon additives, PVDF binder and electrolyte; the electrolyte fills the pore regions. Porosity refers exclusively to the volume fraction of the pore region, i.e., the volume occupied by the electrolyte. For a given dry weight fraction of carbon black, given by w C B , and that of PVDF, given by w PV DF , the porosity ε P is given by, 48 [25] The volume fraction of the cathode particle phase ε C is calculated from the reconstructed geometry. The volume fractions of the carbon black additives ε C B and PVDF binder ε PV DF are then determined from their densities ρ C B and ρ PV DF respectively as, [27] Following, 48 we use ρ C B = 2000 kg/m 3 and ρ PV DF = 1780 kg/m 3 . Figure 6a shows the variation of porosity for each of the 16 electrodes with compression pressure and weight percent of constituents obtained using geometric reconstruction. The error bars represent the variation of the four realizations of each case. The porosities that we obtain using geometry reconstruction match the trends shown in Figure 3 of Ref. 48. Small differences may be attributed to the differences in the geometric reconstruction process. Figure 6a demonstrates the effect of compressive pressure on porosity. Overall, a larger compressive pressure results in a denser particle packing and therefore decreased porosity. Figure 6b re-plots the same data to show variation of porosity with the weight percentages of carbon black (CB) and PVDF (a given value of x has x weight percent of carbon black and PVDF each). Figures 6a and 6b show that compression has a more significant influence on porosity than carbon black and PVDF content, as noted in Ref. 48. Porosities range from 22%-48% across the cases considered.
We also obtained image slices from the Barnett group 47 at Northwestern University (Samples B, C, D). The group performed FIB-SEM scanning on samples procured from commercial lithium-ion batteries; further details of the resolution of the acquired images may be found in Ref. 47. Samples B and C are from uncycled batteries while D is from a highly cycled battery. Unlike the spherical shaped particles of sample A, these samples have particles that have random shape. Only limited results are presented with these cathode materials as we have only one sample for each of them. The volume fraction of the active cathode phase and non-active phase are 70% and 30% respectively for these samples. No information about the manufacturing process is available for these batteries, and therefore the volume fractions of constituents of non-active phase (PVDF/carbon black) cannot be determined. These battery materials are used to give qualitative insight into the effect of microstructure geometry on Bruggeman's exponents.
Effective thermal conductivity.-In order to determine the effective conductivity of the electrode, it is necessary to prescribe the thermal conductivity of the constituent materials. In Ref. 54 55 Since these samples are obtained from commercial electrodes, no information about the volume fraction of the constituents is known. LiCoO 2 and LiNi 1/3 Mn 1/3 Co 1/3 O 2 belong to the class of oxide compounds called layered oxides. To our knowledge, there are no published papers reporting the intrinsic conductivities of these active cathode particles. It should be noted that the thermal conductivity of layered oxide or battery materials depends on the state of charge (i.e., the amount of lithium intercalated in the crystal), as noted in Ref. 24. Some experimental measurements on other layered oxides (sodium, barium, calcium) show considerable change of thermal conductivity with the intercalation of these elements. 56,57 These experiments report the intrinsic thermal conductivities of these crystals, for instance for NaCoO 2 , to vary from 5-15 W/mK with varying intercalation of Na.
Though not exclusively obtained from experiments on lithium-ion batteries, the thermal conductivities of the constituents of the nonactive materials phase are available from other experiments. In this paper, the conductivity of the PVDF binder was taken as 0.2 W/mK, 58 that of the electrolyte as 0.16 W/mK 59 and that of carbon black as 0.2 W/mK. 60,25 In this paper, all variables with subscript C represents active cathode particle phase while the subscript N A represents the non-active material phase. The thermal conductivities and volume fraction of cathode particles and non-active phase are referred as k C , ε C and k N A , ε N A , respectively. For sample A batteries, the volume fractions of the individual non-active constituents were reported in Ref. 48 from the data used to manufacture the battery cathodes. The conductivity of the non-active material phase k N A was approximated as the volume averaged conductivities of the carbon additives, PVDF binder and electrolyte. Therefore, where subscripts C B, PV DF and E stand for carbon additives, PV DF binder and electrolyte phases. This resulted in conductivity k N A of 0.16-0.18 W/mK for various volume fractions of the cathode phase.
In the results that follow, the cathode particle conductivity k C is varied from 0-1000 W/mK. This spans a range for the k C /k N A ratio of 0 − 10 5 . With these values for conductivities and the procedure outlined in section Numerical simulation, the effective conductivity k ef f is determined in all three coordinate directions and averaged. The variation of the volume fraction of the cathode phase or nonactive material phase captures the variation of weight percent and the compression of the 16 cases of sample A. The volume fraction of the cathode phase was found to be in the 0.4-0.7 range for the cases considered.
The normalized effective conductivity k ef f /k C is plotted against the volume fraction of the cathode particle phase ε C for all 16 cases for various values of the ratio k N A /k C in Figure 7. The conductivity of the non-active phase, obtained from Eq. 28, varies slightly for the 16 different cases because of the differences in the volume fractions of the individual components in non-active phase. The error bars shown in the figure are a result of this variability. The value of k N A used for non-dimensionalization is the mean of the non-active materials phase conductivity averaged over all volume fractions, k N A,mean = 0.1691 W/mK. The variance of the normalized effective conductivity plotted in the figure corresponds to the four realizations considered for each case.
The case for k N A /k C → 0 represents the case when heat conduction is solely through the cathode particles. As the volume fraction of cathode particle phase ε C varies from 0.4 to 0.7, the normalized effective thermal conductivity varies from 0.055 to 0.338. These values are significantly lower than unity and also lower than the upper bounds predicted by a variety of effective medium theories (section Comparison with other theories), indicating a tortuous conducting pathway through a cathode phase. Effective conductivity increases with both increasing cathode phase volume fraction and increasing k N A /k C value as shown in Figure 7.
Bruggeman's exponents for effective thermal conductivity. -Equation 12 is rewritten as where d C and d N A are the exponents corresponding to cathode particle phase and non-active material phase respectively. We term these exponents Bruggemann's exponents, understanding that they may differ from the Bruggemann value of 1.5. We first consider the case k N A /k C → ∞, i.e., the case when the cathode particles are non-conducting. Figure 8a shows the variation of d N A . In this limit, the effective thermal conductivity is solely a function of geometry of the microstructure (refer Figure 4b) and the conductivity of the non-active phase. We observe that exponent d N A has a near-constant value of around 1.5 across all the cases. The mean value of d N A across all 16 samples is 1.514, i.e., the effective conductivity of the composite cathode material behaves in accordance with Bruggemann's theory 20 in the limit (Eq. 21) and also agrees with the theory of Sen et al. 40 for spherical particles. A near-constant value across volume fraction suggests that the connectivity of the noncathode phase is well maintained for the various compressions and volume fractions. Furthermore, we note that the exponent of 1.5 holds for particle volume fractions as high as 70%.
We consider next the case k N A /k C → 0, i.e., the case when the non-active phase is non-conducting. Here, the exponent and the effective conductivity are solely a function of the geometry of the microstructure (refer Figure 4a) and the conductivity of cathode particle phase. Compared to d N A , d C values are higher than the Bruggemann exponent of 1.5 across all cases and vary significantly with volume fraction ε C , as seen in Figure 8b. The mean of d C across all 16 samples is 2.849. This suggests that Bruggeman's original formula (of the form of Eq. 21), used widely in volume-averaged battery formulations, 10,28 is inadequate for the characterization of effective thermal conductivity (or other transport properties, see section Relevance to exponents used in electrochemical model equations).
The high value of d C may be interpreted as a loss of connectivity or a high degree of tortuosity of the cathode phase, as discussed in section Bruggeman's exponent and connectivity. Figure 8b suggests that the degree of connectivity changes with changing compression (and correspondingly, volume fraction). (Figure 11 plots the change of cathode connectivity with volume fraction and is useful in interpreting Fig. 8b). For the 5% wt case (w C B = 5%, w PV DF = 5%), the exponent decreases from 3.26 to 2.93 across the range of volume fractions, with the lower volume fractions of the cathode phase resulting in higher exponents and lower connectivity. At higher compressions (higher volume fractions), cathode particles start making better contact with each other, thereby decreasing the value of the exponent. At the other extreme, for 2% wt case (w C B = 2%, w PV DF = 2%), we observe that the connectivity of cathode particles decreases with increase in volume fraction. Though the reason for this opposing trend is not entirely clear, we speculate that it could be the following: The cathode particles are already well connected at 0 GPa because of higher weight percent of cathode material. Since there is not enough non-active phase to absorb the applied pressure, it is likely that the cathode particles start developing cracks at higher compression pressures, as seen in Figure 1j in Reference 48. Cracks increase the length of conduction pathways, thereby increasing the tortuosity and thus increasing d C .
In Figure 9a, we present the variation of the exponent d C with volume fraction for different values of the conductivity ratio k N A /k C , holding d N A at 1.5. Since the variation with respect to volume fraction is small, the average value of d C across all volume fractions is plotted in Figure 9b. As expected, for low k N A /k C , the exponent d C is much higher than unity, but falls to values close to unity as the cathode phase becomes more conducting. When both phases are conducting, Figure 9. (a) Variation of d C with volume fraction of the cathode particle phase ε C for various k N A /k C ratios, and (b) variation of average value of d C with k N A /k C . The exponent d N A is held assumed to be 1.5 in these graphs.
there is little variation of the d C exponent with volume fraction; the primary dependence is on the conductivity ratio. These plots make it clear that using Bruggemann's exponents of 1.5 for the cathode and non-active phases is inappropriate for typical thermal conductivity ratios of interest in battery simulations.
Comparison with other theories.-There have been various attempts to determine analytical expressions for effective thermal conductivity of multiphase materials. For pure heat conduction, the series and parallel models serve as the lower and upper bounds respectively. 61 For κ = k N A /k C , the series model is given by [30] whereas the parallel model is given by The Maxwell-Euken model 61 provides more restrictive lower and upper bounds for effective conductivity which lie between those of the series and parallel models. The lower bound, referred to as Maxwell-Eucken 1 (MW-Lower), is given by, [32] and the upper bound, referred to as Maxwell-Eucken 2 (MW-Upper), is given by, Maxwell-Eucken models are derived by assuming that the particle phase embedded in the matrix phase does not form a continuous network. Particles therefore do not come into contact with neighboring particles. The model also assumes that particle inclusions are spherical in shape and are dispersed sufficiently far apart so that the temperature distribution near the sphere is only local and do not interfere with the temperature field around its neighbor particles. Maxwell-Eucken 1 models the case when the conductivity of the particle phase (the discontinuous phase) is higher than matrix phase (the continuous phase), while Maxwell-Eucken 2 gives the conductivity for the reverse case.
Effective medium theory (EMT), as explicated in Refs. 61-63, is another model which is derived for the effective conductivity of heterogeneous materials where both the phases form continuous media. Its form is given by [34] All the above formulas take into account the volume fraction alone and not the shape or distribution of the particles.
For the range of volume fractions of the cathode phase considered in this paper, all the expressions from Equations 30-34 are plotted in Figures 10a-10d for four ratios, k N A /k C → 0, k N A /k C = 0.02, k N A /k C = 0.11 and k N A /k C → ∞. In addition, we also plot Bruggeman's formula for effective conductivity for finite conductivities of the two phases (Eq. 20) in Figures 10b and 10c. This equation is not valid for k N A /k C → 0 while it defaults to Eq. 21 for k N A /k C → ∞. The normalized thermal conductivities obtained from numerical simulations are also plotted along with the expressions. Exponents in Eq. 29 obtained from the simulations are shown in the legend. We see that EMT and Bruggeman's formula are closest to the volume-averaged fit Eq. 29 for k N A /k C = 0.02 and k N A /k C = 0.11. The former provides an upper bound, while the latter provides a lower bound. Both the curves deviate more from the volume average fit (Eq. 28) when the ratio k N A /k C gets smaller. Though the fit is not perfect, on average EMT and the Bruggeman formula perform better than other theories for the range of k N A /k C appropriate for battery thermal simulations. This is reasonable since both the cathode and non-cathode phases form a continuous network in the battery microstructure. For k N A → 0, our simulations show a k ef f value that increases with the volume fraction of the cathode particle phase. This suggests that a percolative path through the particulate phase does exist even at low volume fractions, though the high value of d c suggests that the particulate phase does not form a strongly connected medium. For k N A /k C → ∞, our computations show good agreement with MW2, and are in relatively good agreement with EMT, again suggesting that the non-active phase forms a well-connected pathway commensurate with a value of d N A = 1.5.
Connectivity and Bruggeman's exponents.-We now examine the connectivities χ N A and χ C defined in section Bruggeman's exponent and connectivity for all the cases considered in this paper. We plot the connectivities for both phases against the volume fraction of the cathode particle phase ε C in Figure 11. We note that the maximum value of the connectivity is unity. With increasing volume fraction of the cathode phase, the connectivity χ N A of the non-active material phase decreases, while that of the cathode phase χ C increases. The rate of increase in connectivity of cathode phase (1.3486) is larger in magnitude than the rate of decrease of non-active phase (−0.748) indicating that active material connectivity is more sensitive to volume fraction. Although the non-active phase connectivity decreases, it maintains a well-connected configuration as shown in Figure 8a where d N A is seen to stay constant for all volume fractions. Additionally, the non-active phase maintains a higher connectivity value than the cathode phase, even when there is more than twice the amount of cathode material.
We also present some results for batteries B, C and D (Table I) for limiting values of k N A /k C to discern their connectivities. From visual inspection, we observe that the particles for samples B, C and D are quite large and have random shapes. A representative sample is shown in Figure 12. The thicknesses of the reconstructed volumes are small and of the order of the particles size. Thus the corresponding volumes cannot strictly be considered REVs. Nevertheless we present these results for first order comparison. All samples from B-D have a cathode particle phase volume fraction of ∼0.7. For comparison, the values for sample A (case 0.2 GPA, 2% − wt) are also provided that match this volume fraction. As discussed previously, battery samples B, C are uncycled while D is highly cycled.
Bruggeman's exponents d C and d N A in the limiting conditions for battery samples B-D are presented in Table II. The cathode phase connectivities for B-D are significantly lower compared to batteries of Sample A as can be seen from the higher values of d C and therefore lower χ C values. The non-active phase connectivities of B-D are

Table II. Values of Bruggemann exponents and connectivity for batteries A-D.
Name Comparing the uncycled batteries B and C and the cycled battery D, it is observed that the connectivity of the non-active phase is similar for all samples, and the value of d N A is close to the Bruggeman value of 1.5. However, we note that the value of d C is very high for samples B-D, and particularly so for the cycled sample D. Similarly, the value of the cathode connectivity χ C is lower for sample D than for the other samples. In reality, cathode particles must be in contact with each other or through the additive particles to close the electronic path. Lithiation and delithiation causes volume expansion and contraction of particles which may cause particles to lose contact with each other, thereby changing thermal pathways (or diffusive pathways) and reducing effective conductivity. 19,4 Relevance to exponents used in electrochemical model equations.-The results presented above may be used to evaluate the correctness of effective properties widely used in volume averaged electrochemical models. As explained in section Volume averaged charge and species transport equations, the effective transport properties used in species and charge transport are given by Eqs. 17-19. These are similar in form to the effective thermal conductivities in Eq. 29 under the limiting conditions k N A /k C → 0 and k N A /k C → ∞.
Consider first the transport of species in the electrolyte phase (Eq. 14), corresponding to path 1 in Figure 2. Here, the intrinsic mass diffusivity of lithium ions in the electrolyte (or non-active region), D E , is large compared to that of lithium in the cathode particles, D S . This corresponds to the limit k N A /k C → ∞ discussed above. We see from Figure 8a that an exponent d E 1.5 is a good approximation to use in Eq. 17. This is consistent with the Bruggeman exponent of 1.5 used in the battery simulation literature. 64,42,34 This value of the exponent is expected to be approximately correct for nearly-spherical or moderately elongated particles, as suggested by Table II. Significant departures from the value of 1.5 may be expected for significantly elongated particles or other geometric arrangements, such as for LiCoO 2 (LCO) and graphite electrodes in. 65,66 Our simulations exhibit a weak dependence of d N A on volume fraction (Figure 8a). It remains to be explored whether the same trend applies to significantly non-spherical particles. The transport of lithium ions due to electric potential is modeled by (Eqs. 14 and 19) and again may be described by the k N A /k C → ∞ scenario. Thus, an exponent d E 1.5 may be used in conjunction with the ionic conductivity κ E as well.
In the absence of carbon particles, electrons flow through the cathode particle phase through path 2, as depicted in Figure 2. In this case, when considering Eq. 13, the electronic conductivity of the electrolyte phase σ E is zero, which corresponds to the case k N A /k C → 0; thus, the exponent in Figure 8b must be used. In contrast to the k N A /k C → ∞ case, this exponent varies with volume fraction and is significantly greater than 1.5 due to poor connectivity between particles. Thus, the use of an exponent d S = 1.5 in Eq. 18, as in, 28,64,34 is inappropriate. This exponent is a strong function of the volume fraction and the manufacturing process, as demonstrated in Table II. If carbon particles significantly changed the electronic conductivity, the exponents presented here would not be valid; in that event, these exponents would provide an upper bound. The methodology presented here may be used to develop more accurate exponents if more accurate scanning techniques were available which could resolve the geometry of carbon particles.

Conclusions
In this paper, we developed a procedure for determining the effective thermal conductivity of two-phase mixtures based on numerical simulations of thermal transport in reconstructed microstructures. We showed that the commonly-used exponent of 1.5 in Bruggeman's model is inappropriatefor thermal conductivity calculations, and that the exponent depends on the conductivity ratio of the non-active to the cathode phase. Furthermore, existing theories for composite thermal conductivity were evaluated against our simulations. It was found that published effective medium theories, though not entirely accurate, provide a reasonable first approximation to direct numerical simulation of the composite. We also showed that the exponents used for effective properties for electrochemical simulations must be calibrated carefully to experiments or simulations, depending on whether ion or electron transport is being considered. An important outcome of the work is the connection between Bruggeman's exponents and the physical connectivity of the medium. We showed that our conductivity calculations in the limit of large and small k N A /k C yield an important quantitative measure, the connectivity, to characterize how continuous the particulate and electrolyte pathways in the cathode are.
A number of extensions of the current work suggest themselves. First, the current work assumes that the electrode is a two-phase material. In reality, porous cathode insertion compounds are more complex multiphase materials made of four phases -active material, carbon particles, binder and the pore space where electrolyte resides. All these constituents have different intrinsic properties and it is still a matter of debate what the right proportions of these constituents are. The procedure outlined using volume averaging can easily be extended to multiphase materials to determine the right volume fraction of constituents needed for efficient and safe batteries. Furthermore, the optimal size and shape of the both active and carbon particles for efficient batteries is not well understood, and thus an area of active research. Recent advances in materials processing have enabled the fabrication of particles with controlled size and shape. The current approach can help optimize particle size and shape to maximize performance. Furthermore, the effective conductivity in this paper was correlated to two parameters -volume fraction and particle size. With access to more diverse microstructure data, it may be possible to achieve a better understanding of the true dependence of transport properties on geometry.
Thus far, we have only discussed steady state thermal problems. Fast charge and discharge rates may involve non-equilibrium thermal transport, thus necessitating multi-temperature volume-averaged models and proper characterization of interfacial heat transfer. Accurate models for interfacial transport may be obtained from unsteady simulations of thermal transport in these geometries. Finally, the modeling of batteries is essentially a coupled electrochemical -thermal problem, and effective properties for the other governing equations may similarly be derived. The procedure described in this paper can be generalized to develop constitutive models from fully resolved microstructure simulations for use in volume-averaged transport models for batteries.