A Three-Dimensional Numerical Assessment of Heterogeneity Impact on a Solid Oxide Fuel Cell’s Anode Performance

Abstract

In this research, a fully three-dimensional, multiphase, microstructure-scale heterogeneous (non-continuous) electrode, Solid Oxide Fuel Cell (SOFC) stack model is implemented in order to assess the impact of homogeneity disturbance in an SOFC anode. The Butler–Volmer model is combined with recent empirical relations for conductivity and aspects of the Maxwell–Boltzmann kinetic theory describing the transport of mass within the porous medium. Methods for the localized quantification of electrode morphology parameters (such as triple phase boundary length) are implemented. The exchange current distribution in the electrode, the partial pressures and the electric potential fields for each phase are computed numerically. In order to simulate heterogeneity, transfer barriers of varying placement and size are added to an otherwise homogeneous, virtual microstructure based on data from FIB-SEM tomography. The results are compared to a model based on the continuous electrode theory, and the points of discrepancy are highlighted.

1. Introduction

Fuel cells are energy conversion devices in which the chemical energy of fuel is converted directly into electrical energy. The Solid Oxide Fuel Cell (SOFC) is the type of fuel cell where the electrolyte is composed of a ceramic material. Due to their high efficiency, high quality waste heat and ability to use a wide variety of fuels, SOFCs are a promising technology proposed for distributed power generation, the reduction of pollution and decreasing fossil fuel consumption [1,2]. An SOFC electrode is made of ceramic-metal composite characterized by a complex multi-phase microstructure. Each phase of the microstructure has a specific transport function. In an SOFC anode, the oxide ions are transported via an oxide phase, the gaseous reagents diffuse through a network of pores and the electronic current generated in the half-reaction flows through the metallic phase. Thus, the SOFC electrode is characterized by high internal complexity. The morphological parameters of the microstructure, such as the boundary between the three phases (triple phase boundary), connectivity, phase tortuosity and volume fraction affect both conductivity and reaction rate, to an extent to which they become a major concern in the design of SOFC electrodes, cells and stacks [3,4,5,6,7,8].

Insight into an SOFC electrode can be gained using FIB-SEM tomography. In FIB-SEM tomography, the sample is subjected to a Focused Ion Beam (FIB), which removes layers of material. The cross-sections may then be imaged using a Scanning Electron Microscope (SEM). Put together, the images may be used to produce a virtual, three-dimensional representation of the anode microstructure [9]. The virtual reconstruction is then subjected to computer graphic algorithms, which allow the quantification of the aforementioned parameters [10]. In the past decade, the method has gained prominence in the SOFC scientific community.

Since computational simulation of transport phenomena may be used to predict the performance of an electrode, it is commonly employed in research and the design of SOFC devices. The models describing the anodic reaction, as well as the transport of charge and reagents in the direction normal to the electrode surface generally employ continuous electrode theory. In those models, the electrode is treated as a continuous medium, in which the transport of all species occurs in parallel. The reaction takes place in the entire computational domain, and the morphological features of the microstructure are modeled as homogeneously distributed macroscopic parameters [3,5,11,12,13]. Attempts at creating fully-three-dimensional, non-continuous models with irregular microstructure geometry have been comparatively rare [14,15,16,17,18]. Models employing the continuous electrode theory, while inexpensive regarding computing power, may be insufficient when heterogeneity influences the process. It is, however, worth noting that when some conditions are met, there exists an analytical solution for transport phenomena in the SOFC anode, as has been demonstrated by Kulikovsky [19].

The virtual reconstruction may be modified or artificially generated. This can be advantageous in simulations aimed at investigating the influence of a specific microstructure feature, such as a mesostructure or an anisotropically-placed grain. Chen et al. have simulated the influence of meso- and macro-pores on reagent concentration and catalyst utilization by editing an artificially-generated pore network [20]. Hsu et al. have compared real and synthetic microstructures to investigate the origin of inhomogeneity in heterogeneous electrodes [21]. In our previous study, we have constructed a finite-volume method-based, three-dimensional, non-continuous model of transport phenomena in a solid oxide fuel cell anode in order to discuss the influence of heterogeneity on electrode performance [18]. The virtual reconstruction of the microstructure was modified in order to include an artificial, non-homogeneous grain-like disturbance. Several such cases were considered, and for each of them, two simulations were performed: one using the aforementioned 3D, the heterogeneous model, and one using a homogeneous, analytical model with the averaged microstructure parameters. Then, the discrepancy between the models was assessed. While the results met qualitative expectations, some questions arose regarding the relationship of the size and placement of the disturbance to the drop in the expected current density.

In the present research, we aim at a more comprehensive parametric study involving different shapes, sizes and placement of the homogeneity-disturbing anomaly. Physically, the anomaly would correspond to a periodically-occurring disturbance, such as agglomerations of impurities, detached grains or a void with no connection to the open porosity. The purpose of the study is to widen the understanding of the behavior of the potential fields and the current distribution in situations where the distribution of the phases within the electrode is inhomogeneous or anisotropic.

2. Microstructure Data

The microstructure discussed in this article has been based on samples from a commercial SOFC anode manufactured by the SOLIDpower Group. The electrode has a high volume fraction of YSZ, low open porosity and insignificant surface diffusivity. The anode microstructure was investigated using FIB-SEM tomography, a technique in which layers of material are removed using Focused Ion Beam (FIB) and then imaged by a Scanning Electron Microscopy (SEM) device (See Figure 1). The captured images, may then be “stacked” to form a 3D reconstruction of the sample. The FIB-SEM tomography of samples discussed in this study was performed at AGHUniversity of Science and Technology’s Academic Center for Materials and Nanotechnology (ACMiN), Krakow, Poland.

3. Mathematical Model

The mathematical model was based on an approach described by us in our previous paper [18]. At the heart of the model lies a set of Poisson equations with a charge transfer rate as a source term, calculated using the Butler–Volmer relationship. The transfer of the charge in the solid phases is based on the following set of equations:

$$\left\{\begin{array}{cc}\hfill 0& =-i+{\displaystyle \frac{\partial }{\partial x}}\left({\sigma }_{\mathrm{el}}{\displaystyle \frac{\partial {\varphi }_{\mathrm{el}}}{\partial x}}\right)+{\displaystyle \frac{\partial }{\partial y}}\left({\sigma }_{\mathrm{el}}{\displaystyle \frac{\partial {\varphi }_{\mathrm{el}}}{\partial y}}\right)+{\displaystyle \frac{\partial }{\partial z}}\left({\sigma }_{\mathrm{el}}{\displaystyle \frac{\partial {\varphi }_{\mathrm{el}}}{\partial z}}\right),\hfill \\ \hfill 0& =i+{\displaystyle \frac{\partial }{\partial x}}\left({\sigma }_{\mathrm{ion}}{\displaystyle \frac{\partial {\varphi }_{\mathrm{ion}}}{\partial x}}\right)+{\displaystyle \frac{\partial }{\partial y}}\left({\sigma }_{\mathrm{ion}}{\displaystyle \frac{\partial {\varphi }_{\mathrm{ion}}}{\partial y}}\right)+{\displaystyle \frac{\partial }{\partial z}}\left({\sigma }_{\mathrm{ion}}{\displaystyle \frac{\partial {\varphi }_{\mathrm{ion}}}{\partial z}}\right),\hfill \\ \hfill 0& =-i+{\displaystyle \frac{\partial }{\partial x}}\left({\displaystyle \frac{2F{D}_{{\mathrm{H}}_2}}{RT}}{\displaystyle \frac{\partial p_{{\mathrm{H}}_2}}{\partial x}}\right)+{\displaystyle \frac{\partial }{\partial y}}\left({\displaystyle \frac{2F{D}_{{\mathrm{H}}_2}}{RT}}{\displaystyle \frac{\partial p_{{\mathrm{H}}_2}}{\partial y}}\right)+{\displaystyle \frac{\partial }{\partial z}}\left({\displaystyle \frac{2F{D}_{{\mathrm{H}}_2}}{RT}}{\displaystyle \frac{\partial p_{{\mathrm{H}}_2}}{\partial z}}\right),\hfill \\ \hfill 0& =i+{\displaystyle \frac{\partial }{\partial x}}\left({\displaystyle \frac{2F{D}_{{\mathrm{H}}_2\mathrm{O}}}{RT}}{\displaystyle \frac{\partial p_{{\mathrm{H}}_2\mathrm{O}}}{\partial x}}\right)+{\displaystyle \frac{\partial }{\partial y}}\left({\displaystyle \frac{2F{D}_{{\mathrm{H}}_2\mathrm{O}}}{RT}}{\displaystyle \frac{\partial p_{{\mathrm{H}}_2\mathrm{O}}}{\partial y}}\right)+{\displaystyle \frac{\partial }{\partial z}}\left({\displaystyle \frac{2F{D}_{{\mathrm{H}}_2\mathrm{O}}}{RT}}{\displaystyle \frac{\partial p_{{\mathrm{H}}_2\mathrm{O}}}{\partial z}}\right),\hfill \end{array}\right.$$

(1)

are electric potentials at electron-conducting and ion-conducting phases respectively, $i\mathrm{A}{\mathrm{m}}^{-3}$ is the local, volumetric charge transfer rate, $x,y,z$ (m) are the Cartesian coordinates, F ($\mathrm{A}\mathrm{s}{\mathrm{mol}}^{-1}$) is the Faraday constant, $R\left(\mathrm{J}{\mathrm{mol}}^{-1}{\mathrm{K}}^{-1}\right)$ is the universal gas constant and $T\left(\mathrm{K}\right)$ is the temperature. ${\sigma }_{\mathrm{el}}$ and ${\sigma }_{\mathrm{ion}}$ (${\mathsf{\Omega }}^{-1}{\mathrm{m}}^{-1}$) are the conductivities of electronic and ionic phases. They are described by the following empirical relationships [22,23]:

$$\left\{\begin{array}{cc}{\sigma }_{\mathrm{el}}\hfill & =3.27\times {10}^6-1065.3T\hfill \\ {\sigma }_{\mathrm{ion}}\hfill & =3.4\times {10}^4\times \exp \left(-{\displaystyle \frac{10350}{T}}\right)\hfill \end{array}\right..$$

(2)

The volumetric charge transfer rate is calculated as follows:

$$i=i_0{l}_{\mathrm{tpb}}\left[exp\left({\displaystyle \frac{\alpha F{\eta }_{\mathrm{act}}}{RT}}\right)-exp\left(-{\displaystyle \frac{(1-\alpha )F{\eta }_{\mathrm{act}}}{RT}}\right)\right],$$

(3)

where $l_{\mathrm{tpb}}\left(\mathrm{m}{\mathrm{m}}^{-3}\right)$ is the TPB density, $\alpha =0.5$ is the charge transfer coefficient (it is assumed that the backward and forward reaction coefficient are equal in value) and $i_0\left(\mathrm{A}{\mathrm{m}}^{-1}\right)$ is the equilibrium exchange current density derived from the data by de Boer [24]:

$$i_0=31.4·p_{{\mathrm{H}}_2}^{-0.03}·p_{{\mathrm{H}}_2\mathrm{O}}^{0.4}·exp\left(-{\displaystyle \frac{152155}{RT}}\right)$$

(4)

where activation overpotential ${\eta }_{\mathrm{act}}\left(\mathrm{V}\right)$ is computed from the following relationship:

$${\eta }_{\mathrm{act}}={\varphi }_{\mathrm{el}}-{\varphi }_{\mathrm{ion}}-{\eta }_{\mathrm{conc}},$$

(5)

in which the concentration overpotential ${\eta }_{\mathrm{conc}}\left(\mathrm{V}\right)$ is obtained as follows:

$${\eta }_{\mathrm{conc}}={\displaystyle \frac{RT}{2F}}ln\left({\displaystyle \frac{p_{{\mathrm{H}}_2}^{\mathrm{bulk}}}{p_{{\mathrm{H}}_2}}}{\displaystyle \frac{p_{{\mathrm{H}}_2\mathrm{O}}}{p_{{\mathrm{H}}_2\mathrm{O}}^{\mathrm{bulk}}}}\right),$$

(6)

where $p_i$ is the local partial pressure of species i in the electrode pores, while $p_i^{\mathrm{bulk}}$ indicates its “bulk” (anodic channel) value. In the employed numerical scheme, it is equivalent to the value at the boundary of the computational domain. A local diffusion coefficient $D_i\left({\mathrm{m}}^2{\mathrm{s}}^{-1}\right)$ is derived from the Fuller–Schettler–Giddings correlation [25]. A correctional term is added to account for the resistance due to the transport phenomena in a transitional regime between the Fickian and the Knudsenian diffusion:

$$D_{\mathrm{i}}^{-1}={\left[{10}^2{\displaystyle \frac{T^{1.75}\left({\displaystyle \frac{1}{M_{{\mathrm{H}}_2}}}+{\displaystyle \frac{1}{M_{{\mathrm{H}}_2\mathrm{O}}}}\right)}{p\left(v_{\mathrm{H}2}^{1/3}+v_{{\mathrm{H}}_2\mathrm{O}}^{1/3}\right)}}\right]}^{-1}+{\left(\frac{d_{\mathrm{pore}}}{3}\sqrt{\frac{8RT}{\pi M_{\mathrm{i}}}}\right)}^{-1},$$

(7)

where $M_i\left(\mathrm{g}{\mathrm{mol}}^{-1}\right)$ is the molar mass of species i, $v_i$ is a specific constant in the Fuller–Schettler–Giddings correlation and $d_{\mathrm{pore}}$ (m) is the pore diameter. The set of adopted boundary conditions is presented in

Table 1. Scheme of boundary conditions.

$\mathit{z}=0$	$\mathit{z}={\mathit{z}}_{\mathbf{b}}$	$\mathit{x}\in \{0,{\mathit{x}}_{\mathbf{b}}\}\vee \mathit{y}\in \{0,{\mathit{y}}_{\mathbf{b}}\}$
$p_{{\mathrm{H}}_2}=p_{{\mathrm{H}}_2,\mathrm{b}}$	$\frac{\partial p_{{\mathrm{H}}_2}}{\partial z}=0$	$\frac{\partial p_{{\mathrm{H}}_2}}{\partial x}=0$	$\frac{\partial p_{{\mathrm{H}}_2}}{\partial y}=0$
$p_{{\mathrm{H}}_2\mathrm{O}}=p_{{\mathrm{H}}_2\mathrm{O},\mathrm{b}}$	$\frac{\partial p_{{{H}}_2\mathrm{O}}}{\partial z}=0$	$\frac{\partial p_{{\mathrm{H}}_2\mathrm{O}}}{\partial x}=0$	$\frac{\partial p_{{\mathrm{H}}_2\mathrm{O}}}{\partial y}=0$
${\varphi }_{\mathrm{el}}={\varphi }_{\mathrm{el},\mathrm{b}}$	$\frac{\partial {\varphi }_{\mathrm{el}}}{\partial z}=0$	$\frac{\partial {\varphi }_{\mathrm{el}}}{\partial x}=0$	$\frac{\partial {\varphi }_{\mathrm{el}}}{\partial y}=0$
$\frac{\partial {\varphi }_{\mathrm{ion}}}{\partial z}=0$	${\varphi }_{\mathrm{ion}}={\varphi }_{\mathrm{ion},\mathrm{b}}$	$\frac{\partial {\varphi }_{\mathrm{ion}}}{\partial x}=0$	$\frac{\partial {\varphi }_{\mathrm{ion}}}{\partial y}=0$

.

Subscript bindicates the value at the boundary of the computational domain. The general solution to the model is obtained using the finite volume method discretization and an SOR (Successive Over-Relaxation)-based solver with local linearization of the source term. The numerical mesh is directly based on microstructure reconstruction. Thus, it is non-continuous: each of the nodes has a phase assigned, corresponding to the phase in the virtual reconstruction. The assigned phase determines which ‘substance’ (electrons, ions and gas reagents) is transported through the node. The exchange current density (the source term) is added only in the nodes adjacent to all three conducting phases (pores, metal and ion-conductor). The mesh, though based on irregular geometry, is regular, uniform and uses cuboid finite volumes. The TPB length is calculated using the volume expansion method [10]. In this study, each finite volume had an edge length of 50 μm. The numerical model and the discretization method are discussed in detail in our previous publication [18].

4. Analytical Solution

Under an assumption that the transport phenomena perpendicular to the anodic plane can be approximated as one-dimensional, the microstructure is largely homogeneous, while the resistance of the metallic phase and concentration losses are negligible, there exists an analytical solution to the described model. The solution was first derived by Kulikovsky et al. [19] and later demonstrated by us [18] to be in good agreement with a 3D, non-continuous model for the aforementioned set of conditions and assumptions. The solution for the potential on the ionic phase and exchange current density is as follows:

$$\left\{\begin{array}{cc}{\varphi }_{\mathrm{ion}}\hfill & ={\displaystyle \frac{RT}{\alpha F}}·\left(1-{\displaystyle \frac{cosh\frac{\widehat{z}}{\widehat{\epsilon }}}{cosh\frac{1}{\widehat{\epsilon }}}}\right)·sinh\frac{\alpha F}{RT}{\varphi }_{\mathrm{b}}\hfill \\ j\hfill & ={\displaystyle \frac{V_{\mathrm{ion}}}{{\tau }_{\mathrm{ion}}}}{\displaystyle \frac{{\sigma }_{\mathrm{ion}}RT}{L\alpha F}}{\displaystyle \frac{1}{\widehat{\epsilon }}}·{\displaystyle \frac{sinh\frac{\widehat{z}}{\widehat{\epsilon }}}{cosh\frac{1}{\widehat{\epsilon }}}}·sinh\frac{\alpha F}{RT}{\varphi }_{\mathrm{b}}\hfill \end{array}\right.,$$

(8)

where $j\left(\mathrm{A}{\mathrm{m}}^{-2}\right)$ is the mean charge transfer rate, $V_i$ is the volume fraction of a phase i and ${\tau }_i$ is the tortuosity of a phase i. $L\left(\mathrm{m}\right)$ is anode thickness, and:

$$\widehat{\epsilon }=\sqrt{\frac{V_{\mathrm{ion}}{\sigma }_{\mathrm{ion}}RT}{2{\tau }_{\mathrm{ion}}{i}_{\mathrm{an}}\alpha F{z}_{\mathrm{b}}^2}}.$$

(9)

5. Microstructure Modification

The goal of the research is to create and test several cases in which the homogeneous microstructure is modified to include an anomaly of homogeneity. The anomaly is a sphere-shaped grain serving no transport or electrochemical function. Physically, this would correspond to an agglomeration of impurities on an unconnected grain or a void space with no connection to the open porosity. The sphere’s radius r is 0.75–$7.5 \mathsf{\mu }\mathrm{m}$, and it is placed at a distance h equal to $7.5 \mathsf{\mu }\mathrm{m}$, $3.75 \mathsf{\mu }\mathrm{m}$ or $0 \mathsf{\mu }\mathrm{m}$ (in which case, the sphere’s center is at the boundary). The disturbance is modeled in two distinct ways: the grain is either fully placed in the computational domain (“center-type” disturbance) or it is located at its edge (“corner-type” disturbance). All of the cases are visualized in Figure 2. The relevant microstructure morphological parameters are listed in

Table 2. Microstructure and disturbance parameters in the cases studied.

Case	h	r	${\mathit{l}}_{\mathbf{tpb}}$	${\mathit{V}}_{\mathbf{ion}}$	${\mathit{\tau }}_{\mathbf{ion}}$	Case	h	r	${\mathit{l}}_{\mathbf{tpb}}$	${\mathit{V}}_{\mathbf{ion}}$	${\mathit{\tau }}_{\mathbf{ion}}$
	${10}^{-6}\mathrm{m}$		$\mathrm{m}{\mathrm{m}}^{-3}$	-	-		${10}^{-6}\mathrm{m}$		$\mathrm{m}{\mathrm{m}}^{-3}$	-	-
0.0.0			5.48 $\times {10}^{12}$	0.44	2.52
center-type (non-conducting)						corner-type (non-conducting)
1.1.1	0.75	0.00	5.48 $\times {10}^{12}$	0.44	2.41	2.1.1	1.50	0.00	5.48 $\times {10}^{12}$	0.44	2.41
1.1.2	1.50	0.00	5.45 $\times {10}^{12}$	0.44	2.42	2.1.2	3.00	0.00	5.44 $\times {10}^{12}$	0.44	2.43
1.1.3	2.35	0.00	5.39 $\times {10}^{12}$	0.44	2.50	2.1.3	4.50	0.00	5.33 $\times {10}^{12}$	0.43	2.56
1.1.4	3.00	0.00	4.27 $\times {10}^{12}$	0.43	2.72	2.1.4	6.00	0.00	5.12 $\times {10}^{12}$	0.41	2.87
1.1.5	3.75	0.00	5.09 $\times {10}^{12}$	0.41	3.29	2.1.5	7.50	0.00	4.75 $\times {10}^{12}$	0.39	3.66
1.2.1	0.75	3.75	5.48 $\times {10}^{12}$	0.44	2.43	2.2.1	1.50	3.75	5.47 $\times {10}^{12}$	0.44	2.46
1.2.2	1.50	3.75	5.43 $\times {10}^{12}$	0.44	2.45	2.2.2	3.00	3.75	5.40 $\times {10}^{12}$	0.43	2.55
1.2.3	2.25	3.75	5.32 $\times {10}^{12}$	0.43	2.62	2.2.3	4.50	3.75	5.21 $\times {10}^{12}$	0.42	2.77
1.2.4	3.00	3.75	5.10 $\times {10}^{12}$	0.41	2.91	2.2.4	6.00	3.75	4.85 $\times {10}^{12}$	0.39	2.67
1.2.5	3.75	3.75	4.74 $\times {10}^{12}$	0.39	4.10	2.2.5	7.50	3.75	4.29 $\times {10}^{12}$	0.36	4.10
1.3.1	0.75	7.00	5.48 $\times {10}^{12}$	0.44	2.46	2.3.1	1.50	7.50	5.47 $\times {10}^{12}$	0.44	2.47
1.3.2	1.50	7.00	5.43 $\times {10}^{12}$	0.44	2.47	2.3.2	3.00	7.50	5.32 $\times {10}^{12}$	0.43	2.50
1.3.3	2.25	7.00	5.32 $\times {10}^{12}$	0.44	2.62	2.3.3	4.50	7.50	5.20 $\times {10}^{12}$	0.42	2.82
1.3.4	3.00	7.00	5.10 $\times {10}^{12}$	0.43	2.72	2.3.4	6.00	7.50	4.78 $\times {10}^{12}$	0.39	2.92
1.3.5	3.75	7.00	4.75 $\times {10}^{12}$	0.41	3.94	2.3.5	7.50	7.50	4.04 $\times {10}^{12}$	0.34	4.47
corner-type (ion-conducting)
2.1.3 (YSZ)	4.50	0	5.28 $\times {10}^{12}$	0.42	2.46
2.1.4 (YSZ)	6.00	0	5.07 $\times {10}^{12}$	0.39	2.50
2.1.5 (YSZ)	7.50	0	4.72 $\times {10}^{12}$	0.34	2.54
2.2.3 (YSZ)	4.50	3.75	5.15 $\times {10}^{12}$	0.42	2.49
2.2.4 (YSZ)	6.00	3.75	4.80 $\times {10}^{12}$	0.39	2.50
2.2.5 (YSZ)	7.50	3.75	4.24 $\times {10}^{12}$	0.34	2.47

.

6. Results and Discussion

The non-modified virtual electrode, as well as the modified virtual microstructures were input into the numerical simulation to solve for potential fields and current density distributions. The results were compared to the semi-analytical solution discussed in Section 4.

The results are presented in Figure 3 (unmodified microstructures), Figure 4 (“center-type” non-conducting disturbance), Figure 5 (“corner-type” non-conducting disturbance) and Figure 6 (“corner-type”, ion-conducting (YSZ) disturbance). The subscript ’an’ signifies the data series obtained from the one-dimensional model using the averaged microstructure parameters. The calculations performed are three-dimensional, and a 3D visualization for an exemplary set of cases is shown in Figure 7. The microstructure parameters input into the analytical model are presented in Table 2. It is worth noting that despite the removal of material, the volume fraction of the ion-conducting phase does not change significantly. Its tortuosity is more affected. However, except for the few cases with large anomalies, it remains within 10% of the value for the unmodified microstructure.

In the case where the anode is homogenous, the non-continuous numerical model and the analytical model with the averaged microstructure parameters are in decent agreement, suggesting that the models perform adequately. Within the common range of applicability, a good agreement between the models has been observed (see Figure 3). According to the falsificationist view on model validation, this implies that the 3D model performs adequately at least in these cases. Thus, we proceed to discuss the cases with the artificially introduced heterogeneity. Model verification has been elaborated on in our previous paper [18].

The simulated anode, being 30 $\mathsf{\mu }$$\mathrm{m}$ thick and having low tortuosity of the ion-conducting phase, is predicted to generate a relatively uniform distribution of exchange current, with very little distinction between the diffusive layer and the catalytic layer.

Both of the result sets for the modified microstructures show very little deviation from the homogeneous model prediction if the anomaly is smaller than 3 $\mathsf{\mu }$$\mathrm{m}$. However, in the cases where the spheres occupy a larger volume of the electrode, one can observe not only a significant drop in current density, but also a growing discrepancy between the numerical simulation and the analytical solution with the averaged microstructure parameters. This discrepancy increases when the sphere radius becomes larger. The homogeneous model predicts a similar impact on performance, due to the decrease of the averaged tortuosity and TPB density. However, it fails to reach an agreement with the numerical simulation regarding the current density (equivalent to the value predicted at maximum depth) distribution, underestimating the implied reaction rate at positions between the disturbance and the anodic channel boundary. This would likely affect the microstructure evolution during prolonged operation of the electrode.

In the numerical simulation, the mean current density increases at a rate comparable to what is seen in the homogeneous case, then stagnates in the vicinity of the disturbance. This suggests that the main factor contributing to the drop is a local decrease in TPB density, rather than additional resistance from throttling of the ionic current pathways. It is worth noting that the drop in current density is lower if the disturbance is located further from the electrolyte. Thus, it appears that even in a relatively thin anode, the volume in proximity to the electrolyte plays a more prominent role. In the cases where the sphere’s center is located at the boundary, the change is not as apparent, because the disturbance’s volume is effectively reduced by half.

In the cases where ion-conducting material has been added, the total current generated in the electrode (the value at the electrolyte boundary)—as computed from the three-dimensional, numerical model—does not diverge significantly from the predictions of the analytical model. Furthermore, the drop in current generation is lower than 10% of the initial value. However, the composition of the overpotential changes: the activation overpotential plays a greater role compared to the reference case (see Figure 3), and the ohmic resistance overpotential becomes smaller. This suggests that lowering the tortuosity of the ion current conduction channels may compensate for a large loss in TPB density. While the total current densities predicted by each of the models are similar, it is worth noting that the current distributions are not. According to the numerical model, the rate of change in proximity to the anodic channel is higher, only to plateau at depths where the YSZ disturbance causes depletion of the TPB.

In all of the simulations performed, the rate of change of the current density is not equal to zero at the anodic channel boundary. This implies that the oxygen ions can penetrate through the entire thickness of the electrode. This suggests that in practice, a thicker anode would be more efficient, due to better utilization of oxygen ions. Since the catalytic layer of the electrode is the focus of this study, considering the large computational complexity of the problem, we have decided not to investigate anodes thicker than 30 $\mathsf{\mu }$$\mathrm{m}$ at this time.

7. Conclusions

In this paper, we used a numerical and analytical transport model to investigate the impact of an anomaly in an otherwise homogeneous anodic microstructure of a thin electrode. An artificial modification of the virtual microstructure was used to create 36 parametric study cases. In each case, we solved for potential fields and current distributions using a three-dimensional non-continuous numerical model, as well as an analytical model with the averaged microstructural parameters. For the given set of boundary conditions, the results from these two models are in agreement when a non-conducting sphere is smaller than 3 $\mathsf{\mu }$$\mathrm{m}$, but become discrepant as the disturbance of the homogeneity becomes larger. Despite accounting for higher tortuosity, lower volume fraction and lower TPB density, the analytical model with the averaged microstructure parameters does not appear to predict the change of the current distribution caused by the added disturbance. In the cases where ion-conducting material was added, the models were in agreement on the total current density, despite predicting different distributions of generation. The results suggest that lowering tortuosity, and thus, ohmic losses, compensates for higher activation overpotential due to the removal of the TPB reaction sites. Even though a thin anode was considered, the impact of electrode volume in the vicinity of the electrolyte is once again highlighted. These points confirm that the homogeneous model with the averaged microstructure parameters is only sufficient when the electrode can be treated as homogeneous and isotropic. Overall, the presented research highlights the importance of using non-continuous, three-dimensional models in microstructure-oriented SOFC electrode design.