Simulation of the Electrochemical Impedance in a Three-Dimensional, Complex Microstructure of Solid Oxide Fuel Cell Cathode and Its Application in the Microstructure Characterization

Electrochemical impedance spectroscopy (EIS) is a powerful technique for material characterization and diagnosis of the solid oxide fuel cells (SOFC) as it enables separation of different phenomena such as bulk diffusion and surface reaction that occur simultaneously in the SOFC. In this work, we simulate the electrochemical impedance in an experimentally determined, three-dimensional (3D) microstructure of a mixed ion-electron conducting (MIEC) SOFC cathode. We determine the impedance response by solving the mass conservation equation in the cathode under the conditions of an AC load across the cathode’s thickness and surface reaction at the pore/solid interface. Our simulation results reveal a need for modifying the Adler-Lane-Steele model, which is widely used for fitting the impedance behavior of a MIEC cathode, to account for the difference in the oscillation amplitudes of the oxygen vacancy concentration at the pore/solid interface and within the solid bulk. Moreover, our results demonstrate that the effective tortuosity is dependent on the frequency of the applied AC load as well as the material properties, and thus the prevalent practice of treating tortuosity as a constant for a given cathode should be revised. Finally, we propose a method of determining the aforementioned dependence of tortuosity on material properties and frequency by using the EIS data.

Introduction

With rising CO₂ levels in the atmosphere, low emission energy technologies for energy conversion and storage are needed to mitigate further increases in the global temperature. Such technologies in active research and development include, Li-ion batteries (Nitta et al., 2015), supercapacitors (Zhang L. et al., 2018), and several types of fuel cells like polymer electrolyte membrane fuel cells (PEMFC) (Wang et al., 2020), solid acid fuel cells (SAFC) (Lim et al., 2020), biofuel cells (BFC) (Shakeel et al., 2019), and solid oxide fuel cells (SOFC) (Mahato et al., 2015). In particular, SOFCs can be used to produce electricity or, when used in reverse (as a solid oxide electrolyzer cell), to produce fuel, depending upon their applied polarization. Specifically, anodically polarized cells act as low-emission fuel-flexible electrochemical engines, and cathodically polarized cells store energy in the form of stable chemical bonds in H₂, CO, or CH₄ through electrolysis (Chen and Jiang, 2016; Mogensen et al., 2019). However, conventional high temperature SOFCs run at 800–1000°C, which makes their operation and maintenance highly cost-ineffective (Gao et al., 2016). The operating temperature has been reduced by the use of electrocatalytically active mixed ion-electron conducting (MIEC) cathodes such as (La,Sr) (Co,Fe)O_3–δ (LSCF) (Hwang et al., 2005; Liu et al., 2013; Niania et al., 2020) and Sr(Ti,Fe)O_3–δ (STF) (Yoo and Bouwmeester, 2012; Perry et al., 2015; Nenning et al., 2017), wherein the oxygen evolution reactions (OERs) and oxygen reduction reactions (ORRs) occur over the entire surface. To take advantage of this behavior, SOFC cathodes are designed to have complex and porous microstructures with a large specific surface area; however, such designs result in a high degree of tortuosity for ion transport, which limits the performance of SOFCs. Thus, it is important to accurately determine the tortuosity of a cathode microstructure and optimize it to enable high performance at low operating temperatures.

Several methods have been reported in the literature for determining the tortuosity of a porous electrode. These methods include: a) porosity-tortuosity relations, such as Bruggeman’s relations (Bruggeman, 1935); b) calculations of the tortuous path length either through the distance propagation method (Chen-Wiegart et al., 2014) or the shortest path search method (Çeçen et al., 2012); and c) calculations of tortuosity based on the effective diffusivity (Wilson et al., 2006). Bruggeman’s relations are known to be less suitable for a domain with connected solid phases and complex porous networks (as in a MIEC cathode) (Tjaden et al., 2018). The methods involving the calculations of tortuous path length or the effective diffusivity require the three-dimensional microstructural data as input, which is obtained through tomography (Wilson et al., 2006; Cooper et al., 2016). Unfortunately, availability of the necessary equipment to obtain tomographic data is not ubiquitous, sample preparation is time consuming, and the calculations are computationally expensive. Furthermore, these approaches provide tortuosity values associated with steady-state diffusion, yielding a single value and do not generally describe the tortuosity relevant to oscillating load or when surface reaction is present. Therefore, there is a need for alternate methods for determining tortuosity of a given SOFC electrode. Yu et al. (Yu et al., 2020) reported that tortuosity can be extracted from the diffusional impedance data; in particular, they showed that in systems limited by only bulk diffusion, the diffusional impedance data can be used to characterize a given microstructure in terms of tortuosity and the area of the loading boundary. They also showed that tortuosity is a function of the AC load frequency, which they define as effective tortuosity. Such use of impedance data to characterize the microstructure of a SOFC cathode is valuable to SOFC researchers since they widely use electrochemical impedance spectroscopy (EIS) to test SOFC performance.

In this study, we show that the EIS data of a SOFC MIEC cathode can be used to determine the tortuosity of the solid phase within the cathode. In particular, we consider a statistically representative (Scott Cronin et al., 2012) portion of an experimentally determined complex three-dimensional microstructure of an unbiased SOFC cathode. In the microstructure, we solve for the amplitude of the concentration-response to an applied AC load under the influence of surface reaction and bulk diffusion, and subsequently, use the solution to determine the impedance behavior. We investigate the effect of different material properties (bulk diffusion coefficient and reaction rate constant) on the impedance behavior, and on the effective tortuosity of the microstructure. From our simulations, we find that, due to the presence of surface reaction, the amplitude of oscillations in the oxygen vacancy concentration at the pore/solid interface of the cathode is lower than the amplitude within the bulk of the solid phase. Moreover, the difference between the two concentration amplitudes increases with an increase in the ratio between the reaction rate constant and the bulk diffusion coefficient. Such a difference is not taken into account by the widely used Adler-Lane-Steele (ALS) model, a macrohomogeneous model developed by Adler, Lane and Steele to predict the impedance data of a SOFC MIEC cathode (Adler et al., 1996; Adler, 1998). Thus, we propose a modification to the ALS model to account for the aforementioned difference.

Furthermore, we develop a method of extracting the effective tortuosity from the impedance data. Specifically, we compare the macrohomogeneous (modified ALS) and 1D Finite Length Gerischer relations for the impedance to determine the effective tortuosity. Our calculations for tortuosity determination reveal that the effective tortuosity is a function of the microstructure, frequency of the applied AC load, and the material properties such as the reaction rate constant and bulk diffusion constant. Additionally, the effective tortuosity of a microstructure in the low-frequency regime (where the tortuosity approaches the DC value) decreases with an increase in the ratio between the reaction rate constant and the bulk diffusion coefficient due to a decreased penetration depth of the electrode reactions. This finding suggests that the prevalent practice of using a single tortuosity value for a given electrode for extracting the material properties by employing the ALS model (Adler, 1998; Zhang S.-L. et al., 2018; Zhang S.-L. et al., 2019) should be reviewed and revised. Additionally, the finding suggests that, both the intrinsic material properties and the microstructure should be considered concurrently in designing a cathode with enhanced performance. Our findings open a new array of applications for the EIS technique in characterizing and optimization of the microstructure of a SOFC cathode.

All the previous reports primarily consider artificially generated microstructures (Rüger et al., 2007; Haffelin et al., 2013; Pereira et al., 2014), which are designed to match given macrohomogeneous properties and thus may not be representative of the local microstructural morphologies and topologies. Only a few three-dimensional impedance calculations have been reported with a consideration of experimentally determined microstructures that are statistically representative of the cathode (Kreller et al., 2011; Lynch et al., 2013). The report by Kreller et al. (Kreller et al., 2011) only explicitly considered the three-dimensional microstructure near the electrode/electrolyte interface and assumed a one-dimensional macrohomogeneous domain beyond a certain distance from the interface. On the other hand, the investigation by Lynch et al. (Lynch et al., 2013) focused on the method of calculating impedance in a microstructure where the ORR occurs either through bulk or surface pathways. However, it did not propose the utilization of the calculated impedance data for characterizing the microstructure nor did it consider the difference between the concentration amplitudes at the pore/solid interface and within the bulk of the solid, which are the two focuses of this paper.

Methods

Model Equations for Impedance Calculation

MIEC cathodes are known to have much higher electronic conductivity as compared to the ionic conductivity (Jung and Tuller, 2009); thus, we ignore any impedance contribution due to electronic resistance. Moreover, due to the high electronic conductivity, it is reasonable to assume that the electrostatic potential of the solid phase is spatially invariant. Thus, in the absence of any gas phase transport limitation, the cathode impedance can be entirely attributed to the reduction of oxygen at the pore/solid interface and the transport of the resulting oxygen ions through the solid in the cathode, which is shown in Figure 1. The oxygen ions diffuse inside the solid phase of the cathode through a vacancy mechanism. Hence, the ionic current can be calculated from the flux of the oxygen vacancies. The mass transport equation for the oxygen vacancies inside the cathode can be written as

$\frac{\partial C_v}{\partial t}\mathrm{=}\nabla \cdot \left(A{D}_v\nabla C_v\right)\in \mathrm{\Omega },(1)$

where $C_v$ represents the concentration of oxygen vacancies, $A$ represents the thermodynamic factor, which is defined by Adler, Lane, and Steele (Adler et al., 1996) as a factor relating the oxygen vacancy concentration to the partial pressure of oxygen gas in the pore phase, $A\mathrm{=}\left(\mathrm{1}/\mathrm{2}\right)\partial \ln \left(P_{O_{\mathrm{2}}}\right)/\partial \ln \left(C_v\right)$. The diffusivity of the vacancies inside the cathode is denoted by $D_v$, whereas $t$ represents time, and $\Omega$ represents the solid phase region within the cathode where the vacancy transport takes place. The reaction at the pore/solid interface within the cathode can be modeled as a first-order chemical reaction as

$\mathrm{-}A{D}_v\nabla C_v\mathrm{=}{\kappa }^{\mathrm{3}D}\left(C_v\mathrm{-}{C}_v^0\right)\in \delta \mathrm{\Omega },(2)$

where ${\kappa }^{\mathrm{3}D}$ is the reaction rate constant, $C_v^0$ is the equilibrium vacancy concentration and $\delta \mathrm{\Omega }$ represents the pore/solid interface.

FIGURE 1

FIGURE 1. A schematic of oxygen transport withing the solid phase (gray region) of a MIEC SOFC cathode; the schematic is based on Kreller, (Kreller, 2011).

We employ the smoothed boundary method (SBM) (Yu et al., 2012) to solve Eq. 1 within an experimentally determined, complex microstructure of the MIEC cathode, along with the boundary condition set in Eq. 2. Using SBM, these equations can be reformulated as

$\psi \frac{\partial C_v}{\partial t}\mathrm{=}\nabla \cdot \left(\mathrm{\psi A}{D}_v\nabla C_v\right)\mathrm{-}\left|\nabla \psi \right|{\kappa }^{\mathrm{3}D}\left(C_v\mathrm{-}{C}_v^0\right),(3)$

where $\psi$ is a domain parameter that represents the geometry of the complex microstructure. Within the pores, $\psi$ has a value of zero, and within the solid phase, it has a value of one. At the pore/solid interface, it smoothly transitions between zero and one. Further details about the process of obtaining $\psi$ for a complex microstructure can be found in our previous work (Yu et al., 2020).

The application of an AC load causes the vacancy concentration, $C_v$, to oscillate about $C_v^0$. The oscillating perturbation is denoted by $\text{Δ}C$, and thus, the relation between $C_v\mathrm{, }{C}_v^0$, and $\mathrm{\Delta }C$ can be written as

$C_v\mathrm{=}{C}_v^0\mathrm{+\Delta }C.(4)$

Since $C_v^0$ is constant throughout space at a given temperature, the spatial variation of $C_v$ can be described entirely in terms of $\mathrm{\Delta }C$ as

$\psi \frac{\partial \mathrm{\Delta C}}{\partial t}\mathrm{=}\nabla \cdot \left(\mathrm{\psi A}{D}_v\nabla \left(\mathrm{\Delta }C\right)\right)\mathrm{-}\left|\nabla \psi \right|{\kappa }^{\mathrm{3}D}\mathrm{\Delta }C.(5)$

Furthermore, $A{D}_v$ can be assumed to be a constant for a MIEC cathode at a given temperature and for a small value of $\mathrm{\Delta }C$ (Mizusaki et al., 1993). Therefore, Eq. 5 can be written as

$\psi \frac{\partial \mathrm{\Delta }C}{\partial t}\mathrm{=}A{D}_v\nabla \cdot \left(\psi \nabla \left(\mathrm{\Delta }C\right)\right)\mathrm{-}\left|\nabla \psi \right|{\kappa }^{\mathrm{3}D}\mathrm{\Delta }C.(6)$

Eq. 6 can be further simplified by expressing the time dependent part of the oscillations as complex exponential functions, i.e.,

$\mathrm{\Delta }C\mathrm{=}\frac{\tilde{C}{e}^{\mathrm{i\omega t}}\mathrm{+}{\tilde{C}}^{\ast }{e}^{\mathrm{-}\mathrm{i\omega t}}}{\mathrm{2}},(7)$

where $\tilde{C}$ is a complex quantity that varies in space with ${\tilde{C}}^{\ast }$ as its complex conjugate, $\omega$ is the frequency of the AC load, and $i$ is the imaginary unit. Let $C_R$ and $C_I$ be the real and imaginary components of $\tilde{C}$. Upon the substitution of Eq. 7 into Eq. 6, the subsequent collection of real and imaginary terms, and the cancellation of the exponential functions, we obtain

$A{D}_v\nabla \cdot \left(\psi \nabla C_R\right)\mathrm{-}\left|\nabla \psi \right|{\kappa }^{\mathrm{3}D}{C}_R\mathrm{=-}\mathrm{\psi \omega }{C}_{\mathit{I}},(\mathrm{8a})$

$A{D}_v\nabla \cdot \left(\psi \nabla C_I\right)\mathrm{-}\left|\nabla \psi \right|{\kappa }^{\mathrm{3}D}{C}_I\mathrm{=}\mathrm{\psi \omega }{C}_R.(\mathrm{8b})$

The above pair of equations is then solved to determine the $C_R$ and $C_I$ as functions of space and $\omega$. Once $C_R$ and $C_I$ are known, the impedance, $Z\left(\omega \right)$, of the cathode can be calculated as (Jacobsen and West, 1995; Nielsen et al., 2011)

$Z\left(\omega \right)=\frac{\mathit{RTL}}{\left(-n^{\mathrm{2}}{F}^{\mathrm{2}}{D}_v{C}_v^0\right)}\frac{\tilde{C}{|}_{\mathrm{x=}0}/L}{\langle \partial \tilde{C}/\partial x{|}_{\mathrm{x=}0}\rangle },(9)$

where $R$ is the universal gas constant, $T$ is temperature, $n$ is the number of moles of electrons consumed in the reaction, $F$ is Faraday’s constant, $L$ is the thickness of the cathode, and $\langle \partial \tilde{C}/\partial x{|}_{x=0}\rangle$ is the gradient of the concentration amplitude in the primary diffusion direction, $x$, averaged over the loading boundary located at $\mathrm{x =}0$, which represents the electrolyte/cathode interface. Furthermore, the first term on the right-hand side of Eq. 9 is the product of the material resistivity and thickness of the cathode, and it can be defined as $Z_0$. Therefore, Eq. 9 becomes

$Z\left(\omega \right)=Z_0\frac{\tilde{C}{|}_{\mathrm{x=}0}/L}{\langle \partial \tilde{C}/\partial x{|}_{\mathrm{x=}0}\rangle }.(10)$

For numerically solving Eqs 8a,b, we chose to make them nondimensional by defining a length scale, $l$; then, $\widehat{\nabla }\mathrm{=}l\nabla$. Hereafter, the circumflex $\left(\widehat{\mathrm{.}}\right)$ symbol denotes that the associated operator or quantity is nondimensional. Substituting $\widehat{\nabla }/\mathrm{1}$ in place of $\nabla$, and multiplying both sides with $\left(l^{\mathrm{2}}/\left(A{D}_v\right)\right)$ gives

$\widehat{\nabla }\cdot \left(\psi \widehat{\nabla }{C}_R\right)\mathrm{-}\left|\widehat{\nabla }\psi \right|{\kappa }^{\mathrm{3}D}\frac{l}{A{D}_v}{C}_R\mathrm{=-}\mathrm{\psi \omega }\frac{l^{\mathrm{2}}}{A{D}_v}{C}_I,(\mathrm{11a})$

$\widehat{\nabla }\cdot \left(\psi \widehat{\nabla }{C}_I\right)\mathrm{-}\left|\widehat{\nabla }\psi \right|{\kappa }^{\mathrm{3}D}\frac{l}{A{D}_v}{C}_I\mathrm{=}\mathrm{\psi \omega }\frac{l^{\mathrm{2}}}{A{D}_v}{C}_R.(\mathrm{11b})$

By defining ${\widehat{\kappa }}^{\mathrm{3}D}={\kappa }^{\mathrm{3}D}l/\left(A{D}_v\right)$ and $\widehat{\mathrm{\omega }}=\omega l^{\mathrm{2}}/\left(A{D}_v\right)$ we obtain

$\widehat{\nabla }·\left(\psi \widehat{\nabla }{C}_R\right)-\left|\widehat{\nabla }\psi \right|{\widehat{\kappa }}^{\mathrm{3}D}{C}_R=-\mathit{\psi }\widehat{\mathrm{\omega }}{C}_I,(\mathrm{12a})$

$\widehat{\nabla }·\left(\mathit{\psi }\widehat{\nabla }{C}_I\right)-\left|\widehat{\nabla }\psi \right|{\widehat{\kappa }}^{\mathrm{3D}}{C}_I=\mathit{\psi }\widehat{\mathrm{\omega }}{C}_R,(\mathrm{12b})$

which are solved for two different sets of boundary conditions (BCs) as listed in Table 1.

TABLE 1

TABLE 1. Sets of boundary conditions applied to the system of equations, Eq.12a,12b.

The two sets of the BCs differ only in those at the cathode/current collector interface $\left(\widehat{x}\mathrm{=}\widehat{L}\right)$, as shown in Table 1. The blocking current collector (BCC) BC represents the case of a foil-like current collector, which blocks all the ionic current at the cathode/current collector interface. This condition is more representative of the SOFC cathodes and is also used by the ALS model (Adler et al., 1996). However, in some reports (Boukamp et al., 2003; Boukamp et al., 2006; Holtappels et al., 2006) the impedance behavior of MIEC oxides has also been fit with the relation that is applicable to the transmissive boundary condition. Under this boundary condition, the vacancy concentration is set to be zero at the cathode/current collector interface, and we will refer to this as transmissive current collector (TCC) BC, which may correspond to a porous or mixed conducting current collector. We simulate the impedance behavior for both of these sets of boundary conditions to cover all the scenarios reported in the literature. The results for TCC BC are provided in Supplementary Material. Moreover, a list of all the variables and symbols used in this article is provided in Table 2.

TABLE 2

TABLE 2. List and description of the variables and symbols used in the model equations.

Numerical Implementation

The choice of SBM enabled the use of a standard Cartesian grid and the finite difference method to solve the above pair of equations (Eqs 12a,b). We used the center difference scheme to discretize the domain, which consists of 352 × 642 × 594 grid points. A uniform grid spacing ($h = 0.0272$, non-dimensional) was selected to ensure the presence of at least 3 grid points across the pore/solid interface ($\lambda = 0.0817$). The selection of the numerical parameters is discussed in Supplementary Material in detail. The Alternating-Direction Line-Relaxation (ADLR) method (Hofhaus and Van de Velde, 1996) was employed to solve the equations. The ADLR method utilizes a tridiagonal matrix solve to obtain the values of the solutions individually along the x, y, and z directions. To simultaneously solve Eqs 12a,b, an iterative scheme was developed, in which the ADLR method was first used to calculate $C_R$ for a given RHS of Eq. 12a, and the resulting values were then used to update the RHS of Eq. 12b. Subsequently, $C_I$ was obtained using Eq. 12b with the updated value of $C_R$. This procedure was repeated until both $C_R$ and $C_I$ numerically converged. We defined the convergence metric as the ratio between the absolute value of the sum of all elements of the change matrix (obtained by taking the difference between the updated concentration and previous concentration) and the sum of all elements of the old concentration matrix. The solution was deemed converged when the metric became less than a specified threshold value. We deduced the value by progressively reducing the value and observing the change in the resulting solution. We found that solution did not change appreciably between the threshold values of $\mathrm{1}{0}^{\mathrm{-10}}$ and $\mathrm{1}{0}^{\mathrm{-11}}$. Thus, we selected a value of $\mathrm{1}{0}^{\mathrm{-10}}$ as the criteria the numerical convergence of our solution.

We chose $\mathrm{l }\mathrm{=}\mathrm{ 0}\mathrm{.46 }\mathrm{\mu m}$ as the length scale for the non-dimensionalization. The reader is referred to our previous work (Yu et al., 2020) for more details on the selection of $l$. Finally, a wide range of values for $A{D}_v$ and ${\kappa }^{\mathrm{3}D}$ are reported in the literature (Zhang L. et al., 2018) because of the use of different cathode materials and operating temperatures (600–800°C). We chose three values of ${\widehat{\kappa }}^{\mathrm{3}D}$ for our study in order to cover this wide range and simulated the impedance behavior for the frequency values, $\widehat{\omega }$, between 0 and 2.672. It should be noted that the value of ${\widehat{\kappa }}^{\mathrm{3}D}$ is affected by the values of both $A{D}_v$ and ${\kappa }^{\mathrm{3}D}$.

Derivation for Impedance Expressions and Transcendental Equation for Tortuosity

Adler, Lane, and Steele proposed a macrohomogeneous model to predict the impedance, $Z^{\mathrm{3}D}$, of a finitely thick MIEC SOFC cathode, which is commonly known as the ALS model (Adler et al., 1996; Adler, 1998). The model gives the impedance response for a symmetric cell with two identical cathodes. Since we consider a half cell with a single cathode in our investigation, the ALS impedance expression is divided by 2. Thus,

$Z^{\mathrm{3}D}\mathrm{=}\frac{R_{\mathrm{chem}}}{\sqrt{\left(\mathrm{1+}\mathrm{i\omega }{t}_{\mathrm{chem}}\right)}}\coth \left(\frac{L}{\delta }\sqrt{\left(\mathrm{1+}\mathrm{i\omega }{t}_{\mathrm{chem}}\right)}\right),(\mathrm{13a})$

$R_{\mathrm{chem}}=\frac{\mathit{RT}}{4F^2}\sqrt{\left(\frac{{\tau }^2}{\left(\mathrm{1-\varepsilon }\right)C_v^0{D}_v\mathrm{a\kappa }{C}_{\mathrm{mc}}}\right)},(\mathrm{13b})$

$t_{\mathrm{chem}}=\frac{C_v^0\left(\mathrm{1-\varepsilon }\right)}{\mathrm{Aa\kappa }{C}_{\mathrm{mc}}},(\mathrm{13c})$

$\delta \mathrm{=}\sqrt{\left(\frac{C_v^0{D}_v\left(\mathrm{1-}\varepsilon \right)}{a{\tau }^{\mathrm{2}}\kappa C_{\mathrm{mc}}}\right)},(\mathrm{13d})$

where $\tau$ is the microstructure tortuosity, $\varepsilon$ is the microstructure porosity, $\kappa$ is the surface exchange coefficient, $a$ is interfacial surface area per unit cathode volume, $C_{mc}$ is the oxygen site concentration (in the unit of mol/m³) in the mixed conductor (mc). The ALS model defines the reaction rate as $\mathrm{A\kappa }{C}_{\mathrm{mc}}\mathrm{\Delta }C/C_v^0$, whereas we define it as ${\kappa }^{\mathrm{3}D}\mathrm{\Delta }C$ (Eq. 2). Therefore, $\kappa$ and ${\kappa }^{3D}$ are related as

${\kappa }^{\mathrm{3D}}\mathrm{ = }\frac{\mathrm{A\kappa }{C}_{\mathrm{mc}}}{C_v^0}.(14)$

It should be noted that in the original ALS expression, $\tau$ is used for denoting tortuosity factor, whereas in this article we use $\tau$ for denoting tortuosity. These two quantities are related; the tortuosity factor is the square of the tortuosity.

The substitution of Eqs 13b,c,d, and Eq. 14 into Eq. 13a gives

$Z^{\mathrm{3D}}=\frac{\mathit{RTL}}{4F^2{C}_v^0{D}_v}\frac{{\tau }^2}{\mathrm{1-\varepsilon }}\frac{\coth \left(\mathit{\tau L}\sqrt{\frac{\frac{a{\kappa }^{\mathrm{3D}}}{\mathrm{1-\varepsilon }}\mathrm{+i\omega }}{A{D}_v}}\right)}{\mathit{\tau L}\sqrt{\frac{\frac{a{\kappa }^{\mathrm{3D}}}{\mathrm{1-\varepsilon }}\mathrm{+i\omega }}{A{D}_v}}}.(15)$

By defining ${\kappa }^{\mathrm{1}D}\mathrm{=}{\kappa }^{\mathrm{3}D}a/\left(\mathrm{1-}\varepsilon \right)$ as the macrohomogeneous reaction rate constant and by using the definition for $Z_0$ we obtain

$Z^{\mathrm{3D}}=Z_0\frac{{\tau }^2}{\mathrm{1-\epsilon }}\frac{\coth \left(\mathit{\tau L}\sqrt{\frac{{\kappa }^{\mathrm{1D}}\mathrm{+i\omega }}{A{D}_v}}\right)}{\mathit{\tau L}\sqrt{\frac{{\kappa }^{\mathrm{1D}}\mathrm{+i\omega }}{A{D}_v}}}.(16)$

Finally, by following the same methodology as described in the model equation section, Eq. 16 can be made nondimensional as

$\frac{Z^{\mathrm{3D}}}{Z_0}=\frac{{\tau }^2}{\mathrm{1-\epsilon }}\frac{\coth \left(\mathit{\tau }\widehat{L}\sqrt{{\widehat{\kappa }}^{\mathrm{1D}}+i\widehat{\mathrm{\omega }}}\right)}{\tau \widehat{L}\sqrt{{\widehat{\kappa }}^{\mathrm{1D}}+i\widehat{\mathrm{\omega }}}},(17)$

where the nondimensional quantities are defined as $\widehat{L}\mathrm{=}L/l$, ${\widehat{\kappa }}^{\mathrm{1}D}\mathrm{=}{\kappa }^{\mathrm{1}D}{l}^{\mathrm{2}}/\left(A{D}_v\right)$, and $\widehat{\mathrm{\omega }}\mathrm{=}\omega l^{\mathrm{2}}/\left(A{D}_v\right)$.

The conversion of Eq. 13 into Eq. 17 enables a direct comparison of the impedance expression of a complex 3D microstructure with that of the standard expression for 1D Finite Length Gerischer (FLG) impedance. The 1D FLG element represents the impedance of a 1D MIEC domain where no microstructural effects are present and where the kinetics is co-limited by both surface reaction and bulk diffusion. The FLG impedance, $Z^{\mathrm{FLG}}$ for a 1D domain with the same macrohomogeneous reaction rate constant, ${\widehat{\kappa }}^{1D}$, and thickness, $\widehat{L}$, as the SOFC cathode under consideration can be written as

$\frac{Z^{\mathrm{FLG}}}{Z_0}=\frac{\coth \left(\widehat{L}\sqrt{{\widehat{\kappa }}^{\mathrm{1D}}+i\widehat{\mathrm{\omega }}}\right)}{\widehat{L}\sqrt{{\widehat{\kappa }}^{\mathrm{1D}}+i\widehat{\mathrm{\omega }}}},(18)$

It is evident from a comparison between Eq. 17 and Eq. 18 that the microstructure affects the impedance response of a SOFC cathode in two ways. First, the microstructure increases the cathode impedance by a factor of ${\tau }^{\mathrm{2}}/\left(\mathrm{1-}\varepsilon \right)$. Second, the effective thickness of the cathode is increased by a factor of $\tau$. These effects can be exploited in the following way to determine the value of $\tau$ using the value of $Z^{3D}$, which is experimentally known.

A ratio between Eq. 17 and Eq. 18 gives

$\frac{Z^{\mathrm{3D}}}{Z^{\mathrm{FLG}}}=\frac{\tau }{\mathrm{1-\varepsilon }}\frac{\coth \left(\mathit{\tau }\widehat{L}\sqrt{{\widehat{\kappa }}^{\mathrm{1D}}+i\widehat{\mathrm{\omega }}}\right)}{\coth (\widehat{L}\sqrt{{\widehat{\kappa }}^{\mathrm{1D}}+i\widehat{\mathrm{\omega }}})}.(19)$

By taking the modulus (which is analytically unnecessary but numerically required) and after some rearrangement, Eq. 19 can be written as

$\left|\mathrm{\tau coth}\left(\mathit{\tau }\widehat{L}\sqrt{{\widehat{\kappa }}^{\mathrm{1D}}+i\widehat{\mathrm{\omega }}}\right)\right|\mathrm{ = }\left|\frac{Z^{\mathrm{3D}}}{Z^{\mathrm{FLG}}}\left(\left(\mathrm{1-\varepsilon }\right)\coth \left(\widehat{L}\sqrt{{\widehat{\kappa }}^{\mathrm{1D}}+i\widehat{\mathrm{\omega }}}\right)\right)\right|(20)$

where terms containing $\tau$ are collected on the left-hand side of the equation, yielding a transcendental equation for $\tau$. The nondimensional macrohomogeneous reaction rate constant, ${\widehat{\kappa }}^{\mathrm{1}D}\mathrm{=}{\kappa }^{\mathrm{1}D}{l}^{\mathrm{2}}/\left(A{D}_v\right)\mathrm{=}a{l}^{\mathrm{2}}{\kappa }^{\mathrm{3}D}/\left(\left(\mathrm{1-}\varepsilon \right)A{D}_v\right)$, is dependent on material properties, namely ${\kappa }^{3D}$ and $A{D}_v$, which can be determined from the literature for standard materials, as well as on microstructural characteristics, namely $a$ and $ϵ$, both of which can be determined experimentally (e.g., using Mercury intrusion porosimeter (Rootare and Prenzlow, 1967). Therefore, with both ${\widehat{\kappa }}^{\mathrm{1}D}$ and $Z^{\mathrm{3}D}$ known, Eq. 20 can be solved to determine the effective tortuosity. However, before implementing this method, we need to make a modification to the ALS model as described in the following section.

Simulation Results and Discussion

Modification to the ALS Model

The ALS model (Adler et al., 1996; Adler, 1998) is derived using the volume averaged form of the mass conservation equation (Eqs 2 and 3), which can be written as

$\left(\mathrm{1-}\epsilon \right)\frac{\partial C_v}{\partial t}\mathrm{=}\frac{A{D}_v\left(\mathrm{1-}\epsilon \right)}{{\tau }^{\mathrm{2}}}\frac{{\partial }^{\mathrm{2}}{C}_v}{\partial x^{\mathrm{2}}}\mathrm{-}a{\kappa }^{\mathrm{3}D}\left(C_v^s\mathrm{-}{C}_v^0\right).(21)$

where $C_v^s$ is the oxygen vacancy concentration at the pore/solid interface. Eq. 21 can be written in term of $\text{Δ}C$ by using the relation described by Eq. 4 as

$\left(\mathrm{1-\varepsilon }\right)\frac{\partial \mathrm{\Delta }C}{\partial t}=\frac{A{D}_v\left(\mathrm{1-\epsilon }\right)}{{\tau }^2}\frac{{\partial }^2\mathrm{\Delta }C}{\partial x^2}\mathit{-a}{\kappa }^{\mathrm{3D}}\mathrm{\Delta }{C}^s.(22)$

The ALS model assumes the concentration oscillations at the pore/solid interface, $\mathrm{\Delta }{C}^s$, and within the bulk of the solid, $\text{Δ}C$, to be equal. However, in reality, the amplitude of the concentration oscillation, $\tilde{C}$, is smaller at the interface than the amplitude within the bulk because of the surface reaction. We demonstrate this phenomenon by considering a model microstructure with solid cylindrical domain surrounded by the pore phase, as shown in Figure 2A. The dimensions of the microstructure are also shown in the figure, and the primary diffusion direction of the model microstructure is along the x-axis. The reaction occurs at the pore/solid interface, i.e., at the surface of the cylinder. Using this model geometry, we determine $C_R$ and $C_I$ within the cylinder for $\widehat{\omega } =0.5$ and for two different values of ${\widehat{\kappa }}^{\mathrm{3}D}$ by solving Eqs 12a, 12b, which are subject to the BBC BC at $\widehat{x}=1$. Figures 2B,C show the distribution of $C_R$ and $C_I$ along the radius of the cylinder at different positions along the x-axis for ${\widehat{\kappa }}^{\mathrm{3}D}\mathrm{ =0}\mathrm{.02}$ and ${\widehat{\kappa }}^{\mathrm{3}D}\mathrm{ =0}\mathrm{.2}$, respectively. The solution values outside the cylinder $\left(\left|\widehat{y}\right|\mathrm{>1}\right)$ have no physical meaning as the domain parameter, $\psi$, is zero there, and therefore are not included in the plots. As can be seen, for a small value of ${\widehat{\kappa }}^{\mathrm{3}D}$, $0.02$, the magnitude of $C_R$ at $\widehat{y}=0$ and $\widehat{y}=1$ are similar (Figure 2B), whereas for a larger value of ${\widehat{\kappa }}^{3D}$, 0.2, the magnitude of $C_R$ is significantly smaller at $\widehat{y}=1$ than at $\widehat{y}=0$, as shown in Figure 2C. The larger extent of spatial variation in $C_R$ for higher ${\widehat{\kappa }}^{\mathrm{3}D}$ is caused by the greater magnitude of $\widehat{\nabla }{C}_R$ at the interface, which is directly proportional to ${\widehat{\kappa }}^{\mathrm{3}D}$. In other words, a large ${\widehat{\kappa }}^{\mathrm{3}D}$ suggests that the rate of diffusion is slower than the reaction rate, which causes $C_R$ to be smaller in magnitude at the interface than within the bulk. This spatial variation is observed throughout the cylinder, and a similar trend is seen for $C_I$ as well.

FIGURE 2

FIGURE 2. (A) The model microstructure, which contains a cylindrical solid domain with unit radius and length surrounded by pore phase. The distribution of $C_R$ and $C_I$ along the radial axis of the cylinder at different $\widehat{x}$ positions for $\widehat{\omega }\mathrm{=0}\mathrm{.5}$ for (B)${\widehat{\kappa }}^{\mathrm{3}D}\mathrm{=0}\mathrm{.02}$, and (C)$0.2$. The solution values outside the cylinder $\left(\left|\widehat{y}\right|\mathrm{>1}\right)$ have no physical meaning as the domain parameter, $\psi$, is zero there, and therefore are not included in the plots. Due to the radial symmetry of the cylinder, results are only shown for $\widehat{y}$ between zero and one.

It is evident from the above model case that the difference between the amplitude of the oscillations at the interface and within the bulk cannot be neglected. Thus, we now consider a modification to the ALS model. Let $\alpha$ be the ratio between the interface and the bulk concentration oscillations, i.e., $\mathrm{\alpha = }\mathrm{\Delta }{C}^s/\mathrm{\Delta }C$. Therefore, Eq. 22 becomes

$\left(\mathrm{1-\varepsilon }\right)\frac{\partial \mathrm{\Delta }C}{\partial t}=\frac{A{D}_v\left(\mathrm{1-\varepsilon }\right)}{{\tau }^2}\frac{{\partial }^2\mathrm{\Delta }C}{\partial x^2}\mathit{-a}{\kappa }^{\mathrm{3D}}\left(\alpha \mathrm{\Delta }C\right).(23)$

Eq. 23 suggests that the effect of $\alpha$ can be incorporated into the ALS model by appropriately modifying ${\widehat{\kappa }}^{\mathrm{3}D}$. However, this cannot be done as straightforwardly because 〈α〉 varies in space. Nonetheless, the complexity can be reduced by considering the volume averaged $\alpha$, or 〈α〉, to modify the value of ${\kappa }^{3D}$. The quantity 〈$\alpha$〉 is defined as the ratio between average interface and average bulk concentration amplitudes, i.e.,

$\langle \alpha \rangle =\frac{{\displaystyle \mathrm{\int }}\tilde{C}\left|\widehat{\nabla }\mathit{\psi }\right|d\widehat{V}}{{\displaystyle \mathrm{\int }}\left|\widehat{\nabla }\mathit{\psi }\right|d\widehat{V}}.{\left(\frac{{\displaystyle \mathrm{\int }}\tilde{C}\mathit{\psi } d\widehat{V}}{{\displaystyle \mathrm{\int }}\mathit{\psi }d\widehat{V}}\right)}^{\mathrm{-1}},(24)$

where the integral is over the entire volume of the cathode microstructure.

To characterize 〈$\alpha$〉, we first determined its dependence on $\widehat{\omega }$ for the three values of ${\widehat{\kappa }}^{\mathrm{3}D}$. The results are summarized in Figure 3. Note that examination of the numerical accuracy is presented in Supplementary Material, and it has been shown that the error in the results is within 2%. It can be seen that, for a particular value of ${\widehat{\kappa }}^{\mathrm{3}D}$, 〈$\alpha$〉 has much higher real component than the imaginary component at all $\widehat{\omega }$ values. Moreover, the real component remains constant with increasing $\widehat{\omega }$ before transitioning into a regime where it decreases with increasing $\widehat{\omega }$. The $\widehat{\omega }$ value where this transition occurs increases with the value of ${\widehat{\kappa }}^{\mathrm{3}D}$. Furthermore, since 〈$\alpha$〉 remains constant for a fairly wide range of frequency values, we made a further simplifying assumption by approximating $\langle \alpha \rangle$ by its DC value, 〈${\alpha }_0$〉. We calculated $\langle {\alpha }_0\rangle$ for 19 values of ${\widehat{\kappa }}^{\mathrm{3}D}$ ranging from $5\times {10}^{-3}$ to $6\times {10}^{-1}$ and summarized the results in Figure 4. The value of 〈${\alpha }_0\rangle \approx \mathrm{1}$ for small values of ${\widehat{\kappa }}^{\mathrm{3}D}$, and it monotonically decreases with increasing ${\widehat{\kappa }}^{\mathrm{3}D}$. Both of these characteristics are expected because small values of ${\widehat{\kappa }}^{\mathrm{3}D}$ suggest that the reaction rate is lower than the diffusion rate, and thus there is little difference between the interface and bulk concentration amplitudes. However, for large values of ${\widehat{\kappa }}^{\mathrm{3}D}$, the reaction rate is much larger than the diffusion rate, which causes the concentration amplitude at the pore/solid interface to be lower than the bulk.

FIGURE 3

FIGURE 3. Calculated 〈$\alpha$〉 vs. $\widehat{\omega }$ (BCC BC) at (A)${\widehat{\kappa }}^{3D} = 0.0231$, (B)$0.05$, and (C)$0.1$. The real component (black solid curve) and the imaginary component (blue dashed curve) of 〈$\alpha$〉 are plotted on the left and right y-axes, respectively. The circles indicate the calculated data point. The DC values of 〈$\alpha$〉 are shown with horizontal magenta dotted lines as a reference.

FIGURE 4

FIGURE 4. DC value of 〈$\alpha$〉, 〈${\alpha }_0$〉, as a function of ${\widehat{\kappa }}^{\mathrm{3}D}$.

Using 〈${\alpha }_0$〉 as a function of ${\widehat{\kappa }}^{\mathrm{3}D}$, to capture the difference between the surface and bulk concentration amplitudes, we propose a modified macrohomogeneous reaction rate constant, ${\widehat{\kappa }}^{\mathrm{1}D}$, in the ALS expression as

$\frac{Z^{\mathrm{3D}}}{Z_0}=\frac{{\tau }^2}{\mathrm{1-\epsilon }}\frac{\coth \left(\mathit{\tau }\widehat{L}\sqrt{\langle {\alpha }_0\rangle {\widehat{\kappa }}^{\mathrm{1D}}+i\widehat{\mathrm{\omega }}}\right)}{\mathit{\tau }\widehat{L}\sqrt{\langle {\alpha }_0\rangle {\widehat{\kappa }}^{\mathrm{1D}}+i\widehat{\mathrm{\omega }}}}.(25)$

Thus, the expression for calculating the effective tortuosity can be modified as

$\left|\mathrm{\tau coth}\left(\mathit{\tau }\widehat{L}\sqrt{\langle {\alpha }_0\rangle {\widehat{\kappa }}^{\mathrm{1D}}+i\widehat{\mathrm{\omega }}}\right)\right|\mathrm{ = }\left|\frac{Z^{\mathrm{3D}}}{Z^{\mathrm{FLG}}}\left(\left(\mathrm{1-\varepsilon }\right)\coth \left(\widehat{L}\sqrt{\langle {\alpha }_0\rangle {\widehat{\kappa }}^{\mathrm{1D}}+i\widehat{\mathrm{\omega }}}\right)\right)\right|.(26)$

Simulated Impedance and Tortuosity Calculation

Figure 5 shows the impedance spectra of the experimentally obtained microstructure for three values of ${\widehat{\kappa }}^{\mathrm{3}D}$, along with the impedance spectra obtained from the analytical expression for the 1D FLG impedance in Eq. 18. As expected, the two curves for the same value of ${\widehat{\kappa }}^{\mathrm{3}D}$ deviate significantly, with the FLG curve underestimating the impedance value. To enable a better visualization of the difference a few iso-frequency points, on both the curves, are marked with magenta triangles. We note that the discussion below is for a fixed cathode thickness provided in the numerical implementation section. The effect of cathode thickness on the impedance and tortuosity is discussed qualitatively in Supplementary Material.

FIGURE 5

FIGURE 5. Nyquist plots obtained from 3D calculations for the BCC BC (black curve) and 1D FLG model (blue curve) for (A)${\widehat{\kappa }}^{3D}=0.0231$, (B)$0.05$, and (C)$0.1$. The iso-frequency points are marked with upright magenta triangles on the black curves and inverted magenta triangles on the blue curves. Three iso-frequency points are connected with arrows in a).

By solving Eq. 26, we calculated $\tau$ for the three values of ${\widehat{\kappa }}^{\mathrm{3}D}$ over the entire frequency range examined. The results are shown in Figure 6. From this comparison, three trends can be observed in the tortuosity data. First, for a given value of ${\widehat{\kappa }}^{\mathrm{3}D}$, $\tau$ remains almost a constant before it begins to decrease with an increase in $\widehat{\omega }$. The relation between $\tau$ and $\log \left(\widehat{\omega }\right)$ at high frequencies can be represented by a linear function with a slope of $-0.1033$, as shown by black dashed line in Figure 6. Second, the low frequency value of $\tau$ decreases with an increase in ${\widehat{\kappa }}^{\mathrm{3}D}$. Finally, all three $\tau$ vs. $\widehat{\omega }$ curves coincide at high frequencies.

FIGURE 6

FIGURE 6. Calculated $\tau$ vs. $\widehat{\omega }$ relations (BCC BC) for ${\widehat{\kappa }}^{3D}=0.0231$ (blue circles), $0.05$ (magenta squares), and $0.1$ (green diamonds). The black dashed line represents the fit with a linear function at high frequencies.

To provide further insights, Figure 7 shows the distribution of $C_R$ and $C_I$ in the complex microstructure at a low value of $\widehat{\omega }$, $0.018$, which is close to the DC case, for three different ${\widehat{\kappa }}^{\mathrm{3}D}$ values. For the low value of ${\widehat{\kappa }}^{\mathrm{3}D}$, $0.0231$, the gradient in the concentration amplitude is non-zero in much of the microstructure thickness, as shown in Figure 7A. However, with an increase in the value of ${\widehat{\kappa }}^{\mathrm{3}D}$, the gradient spans lesser and lesser of the microstructure at the same frequency, as qualitatively shown in Figures 7B,C. Thus, it is evident that the penetration depth of the diffusing species (oxygen vacancies) at the same frequency of the AC load decreases with an increase in ${\widehat{\kappa }}^{\mathrm{3}D}$. Furthermore, Figure 8 shows the distribution of $C_R$ and $C_I$ in the microstructure for three values of ${\widehat{\kappa }}^{\mathrm{3}D}$ at a higher frequency value, $\widehat{\omega } = 1.038$. As can be seen, the gradient in the concentration amplitude is mostly confined to a small region near the electrolyte/cathode interface for all values of ${\widehat{\kappa }}^{\mathrm{3}D}$ at the high frequency regime. Therefore, the penetration depth for each case at the high frequency value is much smaller than at the low frequency value. Moreover, at $\widehat{\omega } = 1.038$ the penetration depth is similar for all three values of ${\widehat{\kappa }}^{\mathrm{3}D}$. This shows that the penetration depth is a function of both ${\widehat{\kappa }}^{\mathrm{3}D}$ and $\widehat{\omega }$. Since the effective tortuosity of a microstructure is directly influenced by the penetration depth (Yu et al., 2020), it is evident that the effective tortuosity is also a function of both ${\widehat{\kappa }}^{\mathrm{3}D}$ and $\widehat{\omega }$. Therefore, it can be inferred that, for a porous medium, where the rate kinetics is co-limited by both the bulk diffusion and surface reaction, the tortuosity is a function not only of the microstructural characteristics but also of ${\widehat{\kappa }}^{\mathrm{3}D}$, which is a combination of material property. In addition, it also depends on the frequency, as was found in the case for diffusional impedance case (Yu et al., 2020). The observed behavior of the effective tortuosity could also be observed from the distribution of reaction-diffusion streamlines (Zhang Y. et al., 2019), which is a 3D visualization of the flux. As the penetration depth decreases, the streamlines or the trajectories of the oxygen vacancies begin to straighten. The resulting decrease in the deviation from a straight path causes the effective tortuosity to decrease.

FIGURE 7

FIGURE 7. The distribution of the real and imaginary components of the concentration amplitude calculated for the BCC BC, $\widehat{\omega } = 0.018$ and (A)${\widehat{\kappa }}^{3D} = 0.0231$, (B)$0.05$, and (C)$0.1$.

FIGURE 8

FIGURE 8. The distribution of the real and imaginary components of the concentration amplitude calculated for the BCC BC, $\widehat{\omega } = 1.038$ and (A)${\widehat{\kappa }}^{3D} = 0.0231$, (B)$0.05$, and (C)$0.1$.

We note that there exists an error of up to 9% in the calculated impedance value at high frequencies, as described in Supplementary Material. The error arises because the accuracy of numerical approximations of the gradients in the concentration amplitudes decreases at high frequencies, as our numerical implementation employs fixed grid resolution throughout the simulation domain. Although this error can be reduced by doubling the grid resolution, such calculations become computationally expensive (as discussed in Supplementary Material) without providing additional insights. In fact, the error does not affect the observed qualitative behavior of the effective tortuosity, which is the one of the two main focuses of this work. Thus, the numerical results presented here are sufficiently accurate for demonstrating the dependence of the effective tortuosity on ${\widehat{\kappa }}^{3D}$ and $\widehat{\omega }$. In the future, we plan to implement the model into a finite element framework with adaptive-mesh capability such as PRISMS-PF framework (Aagesen et al., 2018; DeWitt et al., 2020), in which SBM is already implemented, to increase the computational efficiency of the calculations at high frequencies.

Finally, to further evaluate the accuracy of approximating 〈$\alpha$〉 with 〈${\alpha }_0$〉 in the modified ALS model, we compared the tortuosity results with and without this approximation, latter of which include the frequency dependence of 〈$\alpha$〉. A comparison between the two sets of $\tau$ value obtained for the three values of ${\widehat{\kappa }}^{3D}$ is provided in Figure 9, which show good agreement. Hence the use of 〈${\alpha }_0$〉 as the correction factor is sufficiently accurate for calculating the effective tortuosity of a MIEC SOFC cathode using the EIS data.

FIGURE 9

FIGURE 9. Comparisons of the calculated $\tau$ vs. $\widehat{\omega }$ relations (BCC BC) by using 〈${\alpha }_0$〉 (blue curve) and 〈$\alpha$〉 (black curve) for (A)${\widehat{\kappa }}^{3D} = 0.0231$, (B)$0.05$, and (C)$0.1$.

Summary

In this study, we simulated the impedance behavior of a statistically representative portion of an experimentally determined complex three-dimensional microstructure of an unbiased MIEC SOFC cathode under two different boundary conditions. Our investigation generated two key insights. First, due to the presence of surface reaction, the amplitude of oscillations in the vacancy concentration is lower at the pore/solid interface than within the solid bulk of the cathode. This difference between the interface and bulk amplitude increases with an increase in the ratio between the reaction rate constant and the bulk diffusion coefficient. Therefore, to account for this difference, we proposed a modification to the ALS model in terms of the ratio of the surface and bulk concentration amplitudes at zero frequency, 〈${\alpha }_0$〉, and provided numerically evaluated 〈${\alpha }_0$〉 as a function of the reaction rate coefficient, ${\widehat{\kappa }}^{3D}$. Second, through the examination of the three-dimensional distribution of the vacancy concentration amplitude, we showed that the penetration depth of the oxygen vacancies is a function of ${\widehat{\kappa }}^{3D}$ and the frequency of the applied AC load, $\widehat{\omega }$. Due to the direct dependence of the effective tortuosity on the penetration depth, the effective tortuosity also becomes a function of ${\widehat{\kappa }}^{3D}$ and $\widehat{\omega }$ in addition to the microstructure. Furthermore, we developed a method, which utilizes the EIS data, for determining $\tau$ as a function of ${\widehat{\kappa }}^{3D}$ and $\widehat{\omega }$ for a cathode with known material properties (such as the reaction rate constant and bulk diffusion constant), pore surface area, and porosity.

Data Availability Statement

The data for all the figures in this paper can be downloaded from the following link: https://doi.org/10.13011/m3-ywxf-yw20.

Author Contributions

VG: Methodology, Software, Formal analysis, Visualization, and Writing - original draft. DC: Formal analysis and Writing–original draft. SB: Conceptualization, Formal analysis, Writing–review and editing. KT: Conceptualization, Methodology, Formal analysis, Writing–review and editing, Supervision, Project administration, and Funding acquisition.

Funding

This study is supported by the United States National Science Foundation, Division of Material Research under the grant numbers DMR-1912151 and DMR-1912530. Computational resources were provided by the Extreme Science and Engineering Discovery Environment (XSEDE (Towns et al., 2014)) (allocation No. TG- DMR110007), which is supported by the United States National Science Foundation under grant number ACI-1053575, by the National Energy Research Scientific Computing Center (NERSC), a DOE Office of Science User Facility supported by the Office of Science of the United States Department of Energy under Contract No. DE-AC02–05CH11231.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

The authors would like to acknowledge Prof. Hui-Chia Yu, who was involved in the development of the simulation code used in this work and provided associated information. The authors also thank Dr Beck Andrews for helpful discussions.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fchem.2021.627699/full#supplementary-material

References

Aagesen, L. K., Adams, J. F., Allison, J. E., Andrews, W. B., Araullo-Peters, V., Berman, T., et al. (2018). PRISMS: An Integrated, Open-Source Framework for Accelerating Predictive Structural Materials Science. JOM 70, 2298–2314. doi:10.1007/s11837-018-3079-6

CrossRef Full Text | Google Scholar

Adler, S. B., Lane, J. A., and Steele, B. C. H. (1996). Electrode Kinetics of Porous Mixed‐Conducting Oxygen Electrodes. J. Electrochem. Soc. 143, 3554–3564. doi:10.1149/1.1837252

CrossRef Full Text | Google Scholar

Adler, S. B. (1998). Mechanism and Kinetics of Oxygen Reduction on Porous La1−Sr CoO3− Electrodes. Solid State Ionics 111, 125–134. doi:10.1016/s0167-2738(98)00179-9

CrossRef Full Text | Google Scholar

Boukamp, B., Bouwmeester, H. J. M., Boukamp, B. A., and Bouwmeester, H. J. M. (2003). Interpretation of the Gerischer Impedance in Solid State Ionics. Solid State Ionics 157, 29–33. doi:10.1016/s0167-2738(02)00185-6

CrossRef Full Text | Google Scholar

Boukamp, B., Verbraeken, M., Blank, D., and Holtappels, P. (2006). SOFC-anodes, Proof for a Finite-Length Type Gerischer Impedance?. Solid State Ionics 177, 2539–2541. doi:10.1016/j.ssi.2006.03.002

CrossRef Full Text | Google Scholar

Bruggeman, D. A. G. (1935). Berechnung verschiedener physikalischer Konstanten von heterogenen Substanzen. I. Dielektrizitätskonstanten und Leitfähigkeiten der Mischkörper aus isotropen Substanzen. Ann. Phys. 416, 636–664. doi:10.1002/andp.19354160705

CrossRef Full Text | Google Scholar

Çeçen, A., Wargo, E. A., Hanna, A. C., Turner, D. M., Kalidindi, S. R., and Kumbur, E. C. (2012). 3-D Microstructure Analysis of Fuel Cell Materials: Spatial Distributions of Tortuosity, Void Size and Diffusivity. J. Electrochem. Soc. 159, B299–B307. doi:10.1149/2.068203jes

CrossRef Full Text | Google Scholar

Chen, K., and Jiang, S. P. (2016). Review-Materials Degradation of Solid Oxide Electrolysis Cells. J. Electrochem. Soc. 163, F3070–F3083. doi:10.1149/2.0101611jes

CrossRef Full Text | Google Scholar

Chen-Wiegart, Y.-c. K., Demike, R., Erdonmez, C., Thornton, K., Barnett, S. A., and Wang, J. (2014). Tortuosity Characterization of 3D Microstructure at Nano-Scale for Energy Storage and Conversion Materials. J. Power Sourc. 249, 349–356. doi:10.1016/j.jpowsour.2013.10.026

CrossRef Full Text | Google Scholar

Cooper, S. J., Bertei, A., Shearing, P. R., Kilner, J. A., and Brandon, N. P. (2016). TauFactor: An Open-Source Application for Calculating Tortuosity Factors from Tomographic Data. SoftwareX 5, 203–210. doi:10.1016/j.softx.2016.09.002

CrossRef Full Text | Google Scholar

DeWitt, S., Rudraraju, S., Montiel, D., Andrews, W. B., and Thornton, K. (2020). PRISMS-PF: A General Framework for Phase-Field Modeling with a Matrix-free Finite Element Method. Npj Comput. Mater. 6, 1–12. doi:10.1038/s41524-020-0298-5

CrossRef Full Text | Google Scholar

Gao, Z., Mogni, L. V., Miller, E. C., Railsback, J. G., and Barnett, S. A. (2016). A Perspective on Low-Temperature Solid Oxide Fuel Cells. Energy Environ. Sci. 9, 1602–1644. doi:10.1039/c5ee03858h

CrossRef Full Text | Google Scholar

Häffelin, A., Joos, J., Ender, M., Weber, A., and Ivers-Tiffée, E. (2013). Time-Dependent 3D Impedance Model of Mixed-Conducting Solid Oxide Fuel Cell Cathodes. J. Electrochem. Soc. 160, F867–F876. doi:10.1149/2.093308jes

CrossRef Full Text | Google Scholar

Hofhaus, J., and Van de Velde, E. F. (1996). Alternating-Direction Line-Relaxation Methods on Multicomputers. SIAM J. Sci. Comput. 17, 454–478. doi:10.1137/S1064827593253872

CrossRef Full Text | Google Scholar

Holtappels, P., Verbraeken, M., Vogt, U., Blank, D., and Boukamp, B. (2006). Preparation and Electrochemical Characterisation of Supporting SOFC-Ni-YZT Anodes. Solid State Ionics 177, 2029–2032. doi:10.1016/j.ssi.2006.06.018

CrossRef Full Text | Google Scholar

Hwang, H. J., Moon, J.-W., Lee, S., and Lee, E. A. (2005). Electrochemical Performance of LSCF-Based Composite Cathodes for Intermediate Temperature SOFCs. J. Power Sourc. 145, 243–248. doi:10.1016/j.jpowsour.2005.02.063

CrossRef Full Text | Google Scholar

Jacobsen, T., and West, K. (1995). Diffusion Impedance in Planar, Cylindrical and Spherical Symmetry. Electrochimica Acta 40, 255–262. doi:10.1016/0013-4686(94)E0192-3

CrossRef Full Text | Google Scholar

Jung, W., and Tuller, H. L. (2009). Impedance Study of SrTi1−xFexO3−δ (x=0.05 to 0.80) Mixed Ionic-Electronic Conducting Model Cathode. Solid State Ionics 180, 843–847. doi:10.1016/j.ssi.2009.02.008

CrossRef Full Text | Google Scholar

Kreller, C., Drake, M., Adler, S. B., Chen, H.-Y., Yu, H.-C., Thornton, K., et al. (2011). Modeling SOFC Cathodes Based on 3-D Representations of Electrode Microstructure. ECS Trans. 35, 815–822. doi:10.1149/1.3570062

CrossRef Full Text | Google Scholar

Kreller, C. R. (2011). Measurement and Modeling of Material and Microstructural Factors Governing Performance of Solid Oxide Fuel Cell Cathodes. ProQuest Diss. Theses. Available at: https://proxy.lib.umich.edu/login?url=https://www.proquest.com/dissertations-theses/measurement-modeling-material-microstructural/docview/865807958/se-2?accountid=14667.

Google Scholar

Lim, D.-K., Plymill, A. B., Paik, H., Qian, X., Zecevic, S., Chisholm, C. R. I., et al. (2020). Solid Acid Electrochemical Cell for the Production of Hydrogen from Ammonia. Joule 4, 2338–2347. doi:10.1016/j.joule.2020.10.006

CrossRef Full Text | Google Scholar

Liu, Z., Liu, M., Nie, L., and Liu, M. (2013). Fabrication and Characterization of Functionally-Graded LSCF Cathodes by Tape Casting. Int. J. Hydrogen Energ. 38, 1082–1087. doi:10.1016/j.ijhydene.2012.10.048

CrossRef Full Text | Google Scholar

Lynch, M. E., Ding, D., Harris, W. M., Lombardo, J. J., Nelson, G. J., Chiu, W. K. S., et al. (2013). Flexible Multiphysics Simulation of Porous Electrodes: Conformal to 3D Reconstructed Microstructures. Nano Energy 2, 105–115. doi:10.1016/j.nanoen.2012.08.002

CrossRef Full Text | Google Scholar

Mahato, N., Banerjee, A., Gupta, A., Omar, S., and Balani, K. (2015). Progress in Material Selection for Solid Oxide Fuel Cell Technology: A Review. Prog. Mater. Sci. 72, 141–337. doi:10.1016/j.pmatsci.2015.01.001

CrossRef Full Text | Google Scholar

Mizusaki, J., Yasuda, I., Shimoyama, J. i., Yamauchi, S., and Fueki, K. (1993). Electrical Conductivity, Defect Equilibrium and Oxygen Vacancy Diffusion Coefficient of La1 − X Ca X AlO3 − δ Single Crystals. J. Electrochem. Soc. 140, 467–471. doi:10.1149/1.2221070

CrossRef Full Text | Google Scholar

Mogensen, M. B., Chen, M., Frandsen, H. L., Graves, C., Hansen, J. B., Hansen, K. V., et al. (2019). Reversible Solid-Oxide Cells for Clean and Sustainable Energy. Clean. Energ. 3, 175–201. doi:10.1093/ce/zkz023

CrossRef Full Text | Google Scholar

Nenning, A., Volgger, L., Miller, E., Mogni, L. V., Barnett, S., and Fleig, J. (2017). The Electrochemical Properties of Sr(Ti,Fe)O3-Δfor Anodes in Solid Oxide Fuel Cells. J. Electrochem. Soc. 164, F364–F371. doi:10.1149/2.1271704jes

CrossRef Full Text | Google Scholar

Niania, M. A. R., Rossall, A. K., Van den Berg, J. A., and Kilner, J. A. (2020). The Effect of Sub-surface Strontium Depletion on Oxygen Diffusion in La0.6Sr0.4Co0.2Fe0.8O3−δ. J. Mater. Chem. A. 8, 19414–19424. doi:10.1039/d0ta06058e

CrossRef Full Text | Google Scholar

Nielsen, J., Jacobsen, T., and Wandel, M. (2011). Impedance of Porous IT-SOFC LSCF:CGO Composite Cathodes. Electrochimica Acta 56, 7963–7974. doi:10.1016/j.electacta.2011.05.042

CrossRef Full Text | Google Scholar

Nitta, N., Wu, F., Lee, J. T., and Yushin, G. (2015). Li-ion Battery Materials: Present and Future. Mater. Today 18, 252–264. doi:10.1016/j.mattod.2014.10.040

CrossRef Full Text | Google Scholar

Pereira, J. M. C., Navalho, J. E. P., Amador, A. C. G., and Pereira, J. C. F. (2014). Multi-scale Modeling of Diffusion and Reaction-Diffusion Phenomena in Catalytic Porous Layers: Comparison with the 1D Approach. Chem. Eng. Sci. 117, 364–375. doi:10.1016/j.ces.2014.06.028

CrossRef Full Text | Google Scholar

Perry, N. H., Kim, J. J., Bishop, S. R., and Tuller, H. L. (2015). Strongly Coupled Thermal and Chemical Expansion in the Perovskite Oxide System Sr(Ti,Fe)O3−α. J. Mater. Chem. A. 3, 3602–3611. doi:10.1039/c4ta05247a

CrossRef Full Text | Google Scholar

Rootare, H. M., and Prenzlow, C. F. (1967). Surface Areas from Mercury Porosimeter Measurements. J. Phys. Chem. 71, 2733–2736. doi:10.1021/j100867a057

CrossRef Full Text | Google Scholar

Rüger, B., Weber, A., and Ivers-Tiffee, E. (2007). 3D-Modelling and Performance Evaluation of Mixed Conducting (MIEC) Cathodes. ECS Trans. 7, 2065–2074. doi:10.1149/1.2729320

CrossRef Full Text | Google Scholar

Scott Cronin, J., Muangnapoh, K., Patterson, Z., Yakal-Kremski, K. J., Dravid, V. P., and Barnett, S. A. (2012). Effect of Firing Temperature on LSM-YSZ Composite Cathodes: A Combined Three-Dimensional Microstructure and Impedance Spectroscopy Study. J. Electrochem. Soc. 159, B385–B393. doi:10.1149/2.053204jes

CrossRef Full Text | Google Scholar

Shakeel, N., Ahmad, A., Ahamed, M. I., Inamuddin, , and Asiri, A. M. (2019). Kraton Based Polymeric Nanocomposite Bioanode for the Application in a Biofuel Cell. Enzyme Microb. Technol. 127, 43–49. doi:10.1016/j.enzmictec.2019.04.003

PubMed Abstract | CrossRef Full Text | Google Scholar

Tjaden, B., Brett, D. J. L., and Shearing, P. R. (2018). Tortuosity in Electrochemical Devices: a Review of Calculation Approaches. Int. Mater. Rev. 63, 47–67. doi:10.1080/09506608.2016.1249995

CrossRef Full Text | Google Scholar

Towns, J., Cockerill, T., Dahan, M., Foster, I., Gaither, K., Grimshaw, A., et al. (2014). XSEDE: Accelerating Scientific Discovery. Comput. Sci. Eng. 16, 62–74. doi:10.1109/MCSE.2014.80

CrossRef Full Text | Google Scholar

Wang, Y., Ruiz Diaz, D. F., Chen, K. S., Wang, Z., and Adroher, X. C. (2020). Materials, Technological Status, and Fundamentals of PEM Fuel Cells - A Review. Mater. Today 32, 178–203. doi:10.1016/j.mattod.2019.06.005

CrossRef Full Text | Google Scholar

Wilson, J. R., Kobsiriphat, W., Mendoza, R., Chen, H.-Y., Hiller, J. M., Miller, D. J., et al. (2006). Three-dimensional Reconstruction of a Solid-Oxide Fuel-Cell Anode. Nat. Mater 5, 541–544. doi:10.1038/nmat1668

PubMed Abstract | CrossRef Full Text | Google Scholar

Yoo, C.-Y., and Bouwmeester, H. J. M. (2012). Oxygen Surface Exchange Kinetics of SrTi1−xFexO3−δ Mixed Conducting Oxides. Phys. Chem. Chem. Phys. 14, 11759–11765. doi:10.1039/c2cp41923h

PubMed Abstract | CrossRef Full Text | Google Scholar

Yu, H.-C., Adler, S. B., Barnett, S. A., and Thornton, K. (2020). Simulation of the Diffusional Impedance and Application to the Characterization of Electrodes with Complex Microstructures. Electrochimica Acta 354, 136534. doi:10.1016/j.electacta.2020.136534

CrossRef Full Text | Google Scholar

Yu, H.-C., Chen, H.-Y., and Thornton, K. (2012). Extended Smoothed Boundary Method for Solving Partial Differential Equations with General Boundary Conditions on Complex Boundaries. Model. Simul. Mater. Sci. Eng. 20, 075008. doi:10.1088/0965-0393/20/7/075008

CrossRef Full Text | Google Scholar

Zhang, L., Hu, X., Wang, Z., Sun, F., and Dorrell, D. G. (2018). A Review of Supercapacitor Modeling, Estimation, and Applications: A Control/management Perspective. Renew. Sustain. Energ. Rev. 81, 1868–1878. doi:10.1016/j.rser.2017.05.283

CrossRef Full Text | Google Scholar

Zhang, S.-L., Cox, D., Yang, H., Park, B.-K., Li, C.-X., Li, C.-J., et al. (2019). High Stability SrTi1−xFexO3−δ Electrodes for Oxygen Reduction and Oxygen Evolution Reactions. J. Mater. Chem. A. 7, 21447–21458. doi:10.1039/c9ta07548h

CrossRef Full Text | Google Scholar

Zhang, S.-L., Wang, H., Lu, M. Y., Zhang, A.-P., Mogni, L. V., Liu, Q., et al. (2018). Cobalt-substituted SrTi0.3Fe0.7O3−δ: a Stable High-Performance Oxygen Electrode Material for Intermediate-Temperature Solid Oxide Electrochemical Cells. Energ. Environ. Sci. 11, 1870–1879. doi:10.1039/c8ee00449h

CrossRef Full Text | Google Scholar

Zhang, Y., Yan, F., Ma, J., Yan, M., Ni, M., and Xia, C. (2019). In Operando Monitoring of Reaction-Diffusion Streamlines and Uncovering of Electrochemo-Structural Interactions in Electrodes. J. Mater. Chem. A. 7, 10256–10263. doi:10.1039/c9ta01784d

CrossRef Full Text | Google Scholar