Stochastic 3D Modeling of Three-Phase Microstructures for Predicting Transport Properties: A Case Study

Abstract

We compare two conceptually different stochastic microstructure models, i.e., a graph-based model and a pluri-Gaussian model, that have been introduced to model the transport properties of three-phase microstructures occurring, e.g., in solid oxide fuel cell electrodes. Besides comparing both models, we present new results regarding the relationship between model parameters and certain microstructure characteristics. In particular, an analytical expression is obtained for the expected length of triple phase boundary per unit volume in the pluri-Gaussian model. As a case study, we consider 3D image data which show a representative cutout of a solid oxide fuel cell anode obtained by FIB-SEM tomography. The two models are fitted to image data and compared in terms of morphological characteristics (like mean geodesic tortuosity and constrictivity) as well as in terms of effective transport properties. The Stokes flow in the pore phase and effective conductivities in the solid phases are computed numerically for realizations of the two models as well as for the 3D image data using Fourier methods. The local and effective physical responses of the model realizations are compared to those obtained from 3D image data. Finally, we assess the accuracy of the two methods to predict permeability as well as electronic and ionic conductivities of the anode.

Appendices

Appendix

A Proof of Proposition 1

To prove Proposition 1, we introduce a further stationary random set $\widetilde{\varXi }_2$ defined by $\widetilde{\varXi }_2 = \lbrace t \in \mathbb {R}^3: Y(t) \ge \lambda _Y \rbrace .$ Note that the specific surface areas $S_1$ and $\widetilde{S}_2$ of $\varXi _1$ and $\widetilde{\varXi }_2$, respectively, can be computed by

$$\begin{aligned} S_1 = \frac{2}{\pi } \, e^{-u_{Z}^{2} / 2} \sqrt{-\rho ''_{Z}(0+)} \end{aligned}$$

(20)

and

$$\begin{aligned} \widetilde{S}_2 = \frac{2}{\pi } \, e^{-u_{Y}^{2} / 2} \sqrt{-\rho ''_{Y}(0+)}, \end{aligned}$$

(21)

if Z and Y are mean-square differentiable, see Eq. (6.165) in Chiu et al. (2013) and Remark 7 in Ballani et al. (2012). In that case, $\rho _Z$ and $\rho _Y$ are twice differentiable from the right and $\rho ''_{Z}(0+), \rho ''_{Y}(0+)<0$, see Adler (1981). At first, we derive a formula which allows us to express $L_{\mathrm {TPB}}$ in terms of $S_1$ and $\widetilde{S}_2.$ We show that

$$\begin{aligned} L_{\mathrm {TPB}} =\pi S_1 \widetilde{S}_2 / 4. \end{aligned}$$

(22)

Then, the assertion follows directly, when plugging Eqs. (20) and (21) into Eq. (22). To prove Eq. (22), note that the intersection $\varXi _0 = \varXi _1 \cap \varXi _2 \cap \varXi _3$ is a motion-invariant random closed set and can be considered as a spatial fiber process in the sense of Sect. 8.4 in Chiu et al. (2013). Then, $\varXi _0 \cap [0,1]^2 \times \lbrace o \rbrace$ forms a motion-invariant point process with intensity $\vartheta _0 \ge 0.$ According to Eq. (8.63) in Chiu et al. (2013) we obtain $L_{\mathrm {TPB}} = 2 \vartheta _0.$ Furthermore, note that the intersection of $\partial \varXi _1$ with an arbitrary one-dimensional subset of $\mathbb {R}^3$ forms a motion-invariant point process with intensity $\vartheta _1 \ge 0$ satisfying $S_1 = 2 \vartheta _1,$ see Eq. (8.84) in Chiu et al. (2013). In order to compute $\vartheta _0$, we use the independence of Z and Y. This gives

$$\begin{aligned} \vartheta _0&= \mathbb {E}\mathcal {H}_0(\varXi _0 \cap [0,1]^2 \times \lbrace o \rbrace ) \nonumber \\&= \int \mathbb {E}\left( \mathcal {H}_0(\partial \varXi _1 \cap \partial \widetilde{\varXi }_2 \cap ([0,1]^2 \times \lbrace o \rbrace )) \mid Y\right) \, \mathrm {d}P_Y \nonumber \\&= \vartheta _1 \int \mathcal {H}_1(\partial \widetilde{\varXi }_2 \cap ([0,1]^2 \times \lbrace o \rbrace ) \, \mathrm {d}P_Y \nonumber \\&= \frac{\pi }{8} S_1 \widetilde{S}_2, \end{aligned}$$

(23)

where the last equality is obtained by Eq. (8.83) in Chiu et al. (2013).$\square$

Fig. 6

Current field $J_1$ (b, d) in 2D microstructures (a, c, flow occurs in the complement of the white phase) obtained as the union of the YSZ and pore phases: FIB-SEM image (a, b) and PGM model (c, d). The applied electrical field is $\langle E_1\rangle =1$, $\langle E_2\rangle =0$ and axis 1 is oriented left to right on the maps. The color maps (b, d) indicate the current flow in the direction of axis 1 (lowest value in black, highest in white, values in-between in red and yellow). The color bar is restricted to current values between 0 and 2. Zones of low current values $|J_1|<0.026$ are shown in red in maps (a) and (c)

B Visualization of Dead-End Volume in 2D

The dead-end volume turned out to serve as a reasonable interpretation of the difference observed between models (e.g., the PGM model) and the tomographic image data regarding the M-factor of the nickel phase. In order to visualize the dead-end volume, we have computed the current flow in a hypothetical structure. For this purpose, we considered YSZ and the pore space in a random 2D slice of tomographic image data and of realizations of the PGM model as conducting phases, while nickel is insulating. An electrical field $\langle E_1\rangle =1$ is applied for these 2D structures. This approach is useful for several reasons. First, it is easier to compare the field patterns between different structures in a purely 2D problem. This is however possible as long as discrepancies between the M-factors of model realizations and FIB-SEM images observed in 3D are also reflected in the results of 2D computations. Second, the nickel phase does not percolate in a 2D cut; hence, we consider the flow in the complementary phase of the nickel (YSZ and pores). Indeed if the nickel phase of the models was representative of that observed in the FIB-SEM image, this would be the case also for its complementary. We observe a very significant difference between the two effective conductivities in 2D, equal to $\sigma _{{\mathrm {YSZ,pores}}}=0.21$ for the realization of the PGM model and $\sigma _{{\mathrm {YSZ,pores}}}=0.31$ for the FIB-SEM image. This suggests that the reason for the discrepancy is also present in the simpler 2D problem. In Fig. 6, regions of low current values $|J_1|<0.026$ are highlighted in red. These regions are considered as dead-end volume of the union of pores and YSZ. We observe that there is a significantly larger amount of such regions in realization of the PGM model than in tomographic image data. This can be related to much larger clusters in the model realization for the union of YSZ and of the pores than in tomographic image data (Fig. 6a and b), which act as barriers. The presence of such barriers is consistent with a higher value of effective conductivity $\sigma _{{\mathrm {YSZ,pores}}}=0.31$ of tomographic image compared to the model realization ($\sigma _{{\mathrm {YSZ,pores}}}=0.21$).

C Description of Computer Implementation of the Calculations

In the following, we provide some technical details of the implementations used to simulate the virtual microstructures and their effective transport properties. The simulation of virtual microstructures by the GBM is implemented using Java in the framework of the software library Geostoch (Mayer et al. 2004). Drawing one model realization with the parameters given in Table 1 takes about 25 min on a desktop computer. The code for generating virtual microstructures with the PGM is written in MATLAB (2015). One model realization with the parameters given in Table 2 takes about 5 min on a desktop computer. For the simulation of effective conductivity and permeability as described in Sect. 4.1, a Fortran code parallelized on a 24-cores machine is used. Calculations take about 2 h and 40 min for effective conductivity and permeability, respectively.