Effective Transport Properties of Porous Electrochemical Materials — A Homogenization Approach

We now develop the mathematics for homogenizing the electric conductivity of a three-phase porous electrode to obtain an effective electric conductivity using the analytic continuation method,^2,17,18,26 keeping in mind that the derivation for the effective diffusivity, ionic conductivity, and thermal conductivity are all mathematically identical³⁹ to the one for the effective electric conductivity. Let J and E be the current and electric fields, where E is the negative gradient of the electrical potential ϕ, i.e., . Then, locally, we have following constitutive law

where σ(x) is given in (3). We note that, since the electric field is the negative gradient of the electrical potential and there is no free charge, we can write

For brevity, hereafter we will not indicate the dependence on the independent variable x. Let

where e_k is a unit vector in the k-th direction for k = 1, 2, 3. The effective conductivity, which is a diagonal tensor, is then defined to be

From this we can write

and then define . This allows us to strictly examine the kk-th component of the effective conductivity tensor and simplifies the notation. Thus, using (3) we can rewrite equation (9) as σ* = ⟨e^T_k(σ₁χ₁ + σ₂χ₂ + σ₃χ₃)E⟩. Due to the homogeneity of the effective transport coefficients, i.e., the property that σ*(λσ₁, λσ₂, λσ₃) = λσ*(σ₁, σ₂, σ₃), σ* only depends on the ratios h₁ := σ₁/σ₃ and h₂ := σ₂/σ₃, and we define the function

Although strictly real values σ_j, j = 1, 2, 3, will be used to calculate the effective conductivity of the porous electrode, mathematically, we need not make such restrictions and the σ_j values could in principle be complex numbers. Thus, the two main properties of m(h₁, h₂) are that it is analytic off the interval ( − ∞, 0] in the h₁-plane and h₂-plane, and that it maps the upper half plane to the upper half plane,^1,15,17 so that it is an example of a Herglotz or Stieltjes function. Therefore, to be mathematically guaranteed here that the effective conductivity is well defined and unique for real conductivity values, each material phase must be characterized by a positive, non-zero conductivity σ_j. While in many consistent models, the electrolyte is often assumed to have an electrical conductivity of zero, within the proposed framework an effective conductivity cannot be defined if the conductivity of one of the phases is exactly zero, and therefore we must consider a very small, yet strictly positive conductivity value for the electrolyte. Certain trivial microstructures exist, such as, e.g., parallel slabs, where an effective conductivity tensor can still be uniquely determined if one material has a conductivity identically equal to zero; however, this may not be universally true for all microstructures.

For convenience, let s₁ := 1/(1 − h₁) and s₂ := 1/(1 − h₂) and then define

which is now analytic off [0, 1] in the s₁–plane and s₂–plane. From here, we must now obtain a resolvent representation for E corresponding to (7) and obeying conditions (6). First, it is observed that implies that . Then, we let W be a vector field representing the fluctuations in the electric field, such that ⟨W⟩ = 0, and we suppose that E = e_k + W. Using (3) we thus obtain

After multiple manipulations involving algebraic transformations, application of the operator and the invoking identity χ₁(x) + χ₂(x) + χ₃(x) = 1 true for all x ∈ Ω₀, we obtain

where the inverse Laplacian Δ^{− 1} involves periodic boundary conditions. Following 16–19, we introduce the operator , defined for an arbitrary vector field z as , which allows us to express the resolvent of the electric field E as

in which I is the identity operator. Then, using this representation to express F(s₁, s₂) we find

This formulation for F(s₁, s₂) now has two purposes. The first is that this provides a way to directly calculate the effective conductivity. Recalling that F(s₁, s₂) = 1 − σ*/σ₃ and that is the kk-th component of the effective conductivity tensor, we see that

A quick check indicates that formula (16) collapses to the correct solution if material three and material two occupying, respectively, Ω₃ and Ω₂ have identical properties, thus making it a composite of two materials, and in the trivial case when all three materials have the same properties. A computationally efficient approach to evaluate this expression for experimentally obtained microstructures (described in terms of functions χ₁, χ₂ and χ₃) will be introduced in the next section. The second purpose behind relation (15) is that, as shown below, it can be used to derive rigorous upper and lower bounds on the effective conductivity applicable when only partial information about the microstructure is available.

To obtain the bounds for a three-phase material,¹⁵ a suitable integral representation and an expansion about a homogeneous material must be obtained for equation (15). Here we simply state the rigorous bounds and defer a more complete derivation to Appendix A.

The Wiener or elementary bounds for a three phase material are found using only information about the volume fractions of the three phases, i.e., p₁, p₂, and p₃, along with the individual conductivities of each phase, which are σ₁, σ₂, and σ₃. Assuming real values for σ₁, σ₂, and σ₃, this results in bounds on σ*, where for a three-phase material we have

The microstructure configurations achieving these bounds are discussed in the Analysis and Results section. Since we will also consider two-phase materials in that section, we state here the corresponding rigorous Wiener bounds for the two-phase case, namely

where σ₁ and σ₂ still have real values and p₁ + p₂ = 1.

If more information about the microstructure geometry is available, namely, it can be further assumed to be statistically isotropic, the Hashin-Shtrikman or isotropic bounds can be found. Assuming σ₁ ⩽ σ₂ ⩽ σ₃, the Hashin-Shtrikman bounds for a three-phase material are

which are optimal for certain volume fractions.²⁸ Essentially, the volume fractions of any one phase cannot be extremely small. For a battery microstructure, all phases have volume fractions well within the appropriate range and these bounds (19) are optimal. Although we do not use the rigorous Hashin-Shtrikman bounds for a two-phase material in any subsequent section, we list them here for completeness. For a two-phase material where σ₁ ⩽ σ₂ and p₁ + p₂ = 1, the Hashin-Shtrikman bounds are

which are optimal for all volume fractions.

In order to compute the effective transport coefficients for a battery microstructure, the sequential set of FIB images must be digitally segmented into the three components of a Li-ion battery. A typical FIB image in a cross-sectional plane is shown in Figure 3(A), and the red square (with dimensions 512 × 512 pixels) marks the region on which we will focus. This region is the same for each individual image in the sequential set and truncation of the images to this region ensures that only battery microstructure is captured. Before the images are segmented, we must also consider the boundary conditions imposed by Assumption 1(a) in the previous section that require the battery microstructure to be periodic in all three directions. In principle, a microstructure can always be assumed to be periodic; however, when images are stacked side-by-side, sharp jump discontinuities are observed along the edges (Figure 3(B)). Therefore, as an optional step before segmentation, periodized microstructure images can be smoothed near the edges, reducing the jump discontinuities (point (i) below). Processing of the FIB images is done using an in-house MATLAB code that performs the following sequence of steps:

Zoom In Zoom Out Reset image size

Step (i)— smoothing (optional): a 3D image is reconstituted from the sequential set of 2D SEM images and then smoothed near the edges in all three directions; this is accomplished by treating the gray scale of the 3D image as a function of x₁, x₂, and x₃, and then applying the smooth image periodization algorithm of Moisan;³⁰ it is based on Fourier filtering and produces a 3D image with smoothed edges that is closest to the original 3D image in a suitable sense (Figure 3(B(i))); the smoothed 3D image is then again recast as a sequential set of 2D images; the effect of this optional smoothing step on the computed effective transport coefficients will be studied in the Analysis and Results section,

Step (ii)— thresholding: thresholding combined with gradient analysis is performed using the image processing toolbox available in MATLAB,²³ and allows one to segment the sequential set of periodized 2D images into three distinct phases, see Figure 3(C); more specifically, the in-house algorithm identifies several key aspects of an electrode microstructure taking advantage of the following features: (a) the active particles are relatively smooth, lighter on a gray-scale threshold, and have distinct vertical groves as a result of the FIB process, (b) the pore space (or electrolyte) has a very distinct low value on the gray-scale, and (c) the carbon-doped binder typically has a medium gray-scale value with a "cobweb" structure that has many sharp edges; although the algorithm includes features such as an analysis of gradient values, it is still critical to choose appropriate gray-scale threshold values and any error produced during segmentation is typically a result of using inappropriate threshold values; in addition, if the FIB process has artificially produced large regions where the carbon-doped binder phase has larger (or brighter) gray-scale value than the active particles, the algorithm still performs quite well, but has to be used with greater care — controlling variances in the FIB imaging process is especially important in this regard,

Step (iii)— compression: the sequential set of 2D segmented periodized images is reconstituted into a 3D image and then compressed into a manageable size ensuring that the resulting computational problem is tractable; this downsampling process is done using a nearest neighbor 3D stencil with the selection based on the component that occurs most frequently (an average value stencil cannot be applied here, because three components are being considered); in the present study the final size is 51 × 51 × 51 voxels and the corresponding slice of the 3D image is shown in Figure 3(D).

Although it does not significantly alter the volume fractions of the different phases, the downsampling introduced in step (iii) above may lead to some loss of information, especially as regards transport in the narrowest pathways. While this can be rather difficult to quantify precisely, our tests indicate that downsampling errors are in fact significantly larger than the errors related to FIB imaging and image segmentation. On the other hand, because it involves spatial averaging, downsampling may in the mean sense cancel errors introduced during image segmentation. Providing an accurate characterization of this effect is an interesting open problem for future research.

We now present a computational algorithm which will allow us to conveniently evaluate expression (16) for the effective conductivity by employing standard numerical techniques. The idea is to break the process into steps that can be easily implemented using tools of finite-element analysis, such as available in COMSOL.⁹ Information about the microstructure in the form of the compressed 3D image (Figure 3(D)) is imported as functions χ₁, χ₂ and χ₃ using interpolation routines. Writing the operator out explicitly, equation (14) becomes

Introducing an auxiliary variable (potential) , equation (21) can be recast as the following system

which essentially reduces the problem to the solution of an elliptic boundary-value problem (22a). We reiterate that equation (22a) is subject to periodic boundary conditions (cf. Assumption 1(a)). Then, the kk-th component of the effective conductivity tensor can be easily evaluated as

(an explicit representation of equations (21)–(23) in x, y, and z coordinates is given in Appendix B). Solution of system (22a) is obtained by employing the interactive mode in COMSOL which discretizes the problem using finite elements on a 51 × 51 × 51 grid with the order of interpolation chosen adaptively. The resulting algebraic system is then solved using the solver MUMPS available in COMSOL which can take advantage of shared-memory multicore architectures. Due to the high variability of the coefficients in system (22a), convergence of iterations needs to be closely monitored. In all results reported in the next section, convergence needed to determine the effective transport coefficients with the required number of the significant digits was typically achieved within 25 iterations. An average run takes approximately two hours using 32 2.4 GHz AMD Opteron processors and requires about 35 Gb of RAM memory.

It is in this fashion that the effective conductivity tensor is found for any geometric configuration that has three phases, each with its own individual conductivity values. This solution is found without having to specifically track the electric and current fields within the composite material; however, these fields can be studied a posteriori via suitable post-processing (cf. Figure 3(E)). We note that while the effective transport coefficients can be computed in a variety of ways, the proposed approach has the advantage of separating the geometric information about the microstructure from the material properties. The expressions used to derive the lower and upper bounds in the section on homogenization also share this property.

Although one of the key novelties of this study is the characterization of three-phase composites, some of the essential features of the proposed approach are already evident in the two-phase case which is easier to illustrate. For clarity, we will thus use the two-phase setting to make some of the key points.

In order to demonstrate the relationship between the microstructure geometry and the corresponding effective transport coefficients, Figure 4 analyzes three hypothetical composite materials in 2D, where each composite has identical volume fractions, but a different microstructure geometry. We see that the corresponding effective conductivities are in fact quite different and span a broad range of values between the rigorous lower and upper bounds, which are the corresponding Wiener bounds (18) for a two-phase composite. These bounds on the achievable range of effective transport coefficients are, in turn, quite different from the conductivities of the two materials. We also emphasize that two of the microstructures yield effective transport coefficients quite different from the constant value predicted by Bruggeman's formula (1).

Zoom In Zoom Out Reset image size

While the results presented in Figure 4 correspond to a fixed contrast h = σ₁/σ₂ = 10, it is also interesting to evaluate how the allowed range of effective transport coefficients varies with contrast h. To illustrate this, Bruggeman's formula (1) together with the upper (17) and lower (18) two-phase Wiener bounds can be rewritten as functions of the contrast h in the general form σ* = σ₂ f(h) as follows

These relations are shown in Figure 5(A) for arbitrary volume fractions of p₁ = p₂ = 1/2 and conductivity σ₂ = 1, and in Figure 5(B) corresponding to volume fractions from our sampled electrode (Table I) and an approximate electrolyte conductivity from the literature (Table II). Here, phase 1 is assumed to be the combined volume fractions of the active particles and carbon-doped binder, whereas phase 2 is the volume fraction of the electrolyte. Thus, p₁ = 0.5953 and p₂ = 0.4047, respectively, and σ₂ = 1.0 × 10^{− 8} S · cm^{− 1}. In Figures 5(A) and 5(B) we see that the gap between the upper and lower Wiener bounds increases in proportion to h and, as expected, the two bounds collapse for h = 1, which corresponds to the homogeneous material (σ₁ = σ₂). We also note that for large contrasts, the predictions of Bruggeman's formula are quite close to the upper bound. On the other hand for low contrasts, these predictions deviate from the upper bound and, for certain values of σ₁ > σ₂, become incorrect as they drop below the lower bound. Predictions of Bruggeman's formula remain unfeasible when σ₁ ⩽ σ₂, which is a consequence of the assumption underlying Bruggeman's formula that phase 1 is dominating.^4,39

Zoom In Zoom Out Reset image size

In this Section we use the framework developed above to compute the effective electric conductivity and diffusivity of a porous electrode of an actual Li-ion battery. This battery electrode is assumed to consist of three phases, namely, the electrolyte, active particles and carbon-doped binder. The methodology employed to obtain the microstructure images used in our study was described in the previous section. The volume fractions of the three phases obtained after image processing with and without optional smoothing are listed in Table I. As regards the transport properties of the individual phases, even significant uncertainties are not uncommon in the literature,³³ and one of the goals of this study is to quantify the effect of this uncertainty on the effective transport coefficients computed with the homogenization approach. Another objective is to assess the importance of smoothing the periodized microstructure (cf. step (i) in the previous section). To address these issues in a systematic manner, we assume that the transport coefficients of the dominating phase (the electrolyte for the Li-ion diffusion and the binder with carbon for the electric conduction) are known exactly, and a range of different values is considered for the transport coefficients of the remaining phases. The calculations are performed with smoothed and nonsmoothed microstructure which leads to six distinct cases (A,D), (B,E) and

(C,F) corresponding to decreasing values of the transport coefficients in the nondominating phases, in the non-smoothed and smoothed setting. The data characterizing the different cases is collected in Table II, whereas the computed effective transport coefficients are listed in Table III together with the associated lower and upper bounds and the predictions based on Bruggeman's formula (1). As regards the effective electrical conductivity, we also report the results of applying Bruggeman's formula with the assumption that the binder and active particle are combined into one phase, so that

in relation (1). Homogenization calculations reported in Table III were performed assuming transport in the direction x₁ (i.e., k = 1 in equation (16)). To assess to what extent the microstructure considered here is actually isotropic, in CASE B we performed the calculations assuming independent transport in the three different directions. The obtained standard deviation is less than 6% of the mean, indicating that the microstructure anisotropy does not in fact play a significant role. The results for cases A, B and C are plotted in Figure 6 employing the logarithmic scale, useful to assess how the effective properties compare with the rigorous bounds, and also the linear scale, to emphasize the differences between the different approaches (since the results obtained with the smoothed periodized microstructure differ very little, we do not present separate plots for cases D, E and F).

Zoom In Zoom Out Reset image size

In Table III and Figure 6 we observe that the computed effective diffusivities and conductivities are in the different cases within 20%–30% of saturating the isotropic upper bounds, meaning that the microstructure sample considered in our study is close to having optimal transport properties. As regards the effect of the uncertainty characterizing the properties of the individual phases, we note that even though the diffusivities of the binder and the active particle are allowed to vary by several orders of magnitude, the resulting effective diffusivity changes only by less than 5%. On the other hand, this effect is more pronounced in regard to the effective electrical conductivity which changes by more than 40% when the conductivities of the electrolyte and the binder are decreased (cf. Table II). Smoothing of the periodized microstructure has a small effect decreasing the effective Li-ion diffusivity (by about 5%) and increasing the effective electrical conductivity (by about 6%). Finally, concerning comparison with Bruggeman's formula (1), we observe that this formula tends to overestimate the effective Li-ion diffusivity by up to 14% and underestimate the effective electrical conductivity by up to 50%, with lumped formula (27) yielding even worse results.

To obtain the bounds for a three phase material,¹⁵ let us first consider the function analogous to (11) for a two-phase material,¹ namely G(s) := 1 − m(h), where s := 1/(1 − h) and h := σ₁/σ₂, which is analytic off [0, 1] in the s–plane. Here we briefly consider a two-phase material, because the integral representation is easier to define and can then be easily extrapolated to a three-phase material. It was proven^1,17 that G(s) has the following representation

where μ is a positive measure on [0, 1]. Formula (A1) separates the information about the material properties encoded in s from information about the microstructure geometry contained in μ, a spectral measure of the operator . Information about the geometry, expressed though a number of statistical assumptions, is incorporated into μ via its moments μ_n = ∫¹₀zⁿdμ(z) which can be calculated from the correlation functions of the composite medium as . This integral representation will now be adapted below to find bounds for a three-phase material.

For simplicity, we now assume that the σ_j values are strictly real (derivation of the general complex-valued case can be found in Ref. 15). When |s₁| > 1 and |s₂| > 1, F(s₁, s₂) has the following expansion about the point s₁ = ∞ and s₂ = ∞, which characterizes a homogeneous medium,

in which the coefficients are given in terms of the correlation functions and where q = 0, 1, 2, 3 and w = 0, 1, 2, 3 such that q + w = 3. Rewriting this expression to separate terms containing only s₁ and χ₁, or s₂ and χ₂, we get

where

Thus, F₁(s₁) and F₂(s₂) have the integral representations, similar to equation (A1),

where μ₁ and μ₂ are positive measures on [0, 1]. Similar to the two-component material version (A1), these integral representations (A5) separate the information about the material properties contained in s₁ and s₂ from information about the composite geometry contained in μ₁ and μ₂, which are spectral measures of the operators and , respectively. Again, statistical assumptions about the geometry are incorporated into μ₁ and μ₂ via their moments μ_{1, n} = ∫¹₀zⁿ₁μ₁(dz₁) and μ_{2, n} = ∫¹₀zⁿ₂μ₂(dz₂), which can be calculated based on the correlation functions of the composite medium, with and .

Wiener Bounds

To obtain the elementary bounds for a three-phase material, only the volume fractions of the three components need to be known along with the material properties incorporated into s₁ and s₂. Letting p₁, p₂, and p₃ = 1 − p₁ − p₂ be the volume fractions of the three materials, with the properties ⟨χ₁⟩ = p₁ and ⟨χ₂⟩ = p₂, we observe that F(s₁, s₂) is known only to first order, and the second order terms, including the cross-term, in equations (A2) and (A3) cannot be evaluated. Thus, F(s₁, s₂) to first order is strictly expressed by the first-order truncations of F₁(s₁) and F₂(s₂). To find the upper bound on σ*, we now minimize F(s₁, s₂) with respect to all admissible geometries, cf. (10)–(11), and since for fixed p₁ and p₂ it is known to first order, then F(s₁, s₂) is minimized by separately minimizing the first-order truncations of F₁(s₁) and F₂(s₂) with respect to all admissible geometries represented by measures dμ₁ and dμ₂

The minimums for F₁(s₁) and F₂(s₂) are then obtained by letting μ₁ = p₁δ₀ and μ₂ = p₂δ₀¹⁵ in relations (A5). Thus,

To find the lower bound, consider the auxiliary function¹

which has the same analytical properties as F(s₁, s₂), where t₁ := 1 − s₁ and t₂ := 1 − s₂. Similar to F(s₁, s₂), has the first-order truncated expansion of

and thus using the same argument as before

Putting these bounds together results in the elementary, or Wiener, bounds on σ* for a three component material

The microstructure configurations achieving these bounds are discussed in the Analysis and Results section.

Hashin-Shtrikman Bounds

If more information about the microstructure geometry is available, namely, it can be further assumed to be statistically isotropic, additional correlation terms can be found in expansion (A2) for F(s₁, s₂) which then becomes¹⁵

where d is the space dimension (for a battery microstructure d = 3) and again where q = 0, 1, 2, 3 and w = 0, 1, 2, 3 such that q + w = 3. Since the expansion truncated at the second order now has the non-zero cross-term (2p₁p₂)/(ds₁s₂), it is more convenient to consider the relation

where H(s₁, s₂) again has the same analytic properties as F(s₁, s₂), but now has the expansion

in which the cross-term with the denominator (ds₁s₂) has been eliminated under this transformation. Thus, similar to before, H(s₁, s₂) = H₁(s₁) + H₂(s₂) + ... and the minimum of H(s₁, s₂) is achieved when both H₁(s₁) and H₂(s₂) separately achieve minimums. A convenient way of including the second-order expansion terms so that minimums can be found³ is to use the transformations J₁(s₁) := 1/p₁ − 1/(s₁H₁(s₁)) and J₂(s₂) := 1/p₂ − 1/(s₂H₂(s₂)) which again have the same analytic properties as H₁(s₁) and H₂(s₂). Therefore, under these transformations, the second-order terms reduce to first-order terms only, and

Then, using the analogous integral representations as in equation (A5), the minimum for J₁(s₁) and J₂(s₂) with respect to all admissible geometries are then obtained with the measures δ₀/(dp₁) and δ₀/(dp₂).¹⁵ Finally, mapping back to H₁(s₁) and H₂(s₂), we observe that

To obtain the other bound, we use another auxiliary function ,¹ where y₂ := 1/(1 − σ₂/σ₁) and y₃ := 1/(1 − σ₃/σ₁). Then, in analogy to the argument above, the following relation is found

In this way, when σ₁ ⩽ σ₂ ⩽ σ₃, the Hashin-Shtrikman bounds for a three component material are found

which are optimal for certain volume fractions.²⁸ Essentially, the volume fractions of any one phase cannot be extremely small. For a battery microstructure, all phases have volume fractions well within the appropriate range and these bounds (A18) are optimal.

For completeness, an explicit representation of equations (21)–(23) is given below in terms of their x, y, and z components. Rewriting operator in terms of its components, letting C := (s^{− 1}₁χ₁ + s^{− 1}₂χ₂), and denoting δ_{α, k} = 1 when α = k and zero otherwise, vector equation (21) can be recast as the following system of three scalar equations

where E = [E_x, E_y, E_z]. Introducing an auxiliary variable (potential) ψ := ( − Δ)^{− 1}[∂_x(CE_x) + ∂_y(CE_y) + ∂_z(CE_z)], equations (B1) can be rewritten as the following system of four equations

Finally, the diagonal components of the effective conductivity tensor can be evaluated as