Characterization of gas diffusion electrodes for metal-air batteries

TL;DR

This paper develops advanced 2D and 3D multiphase Lattice-Boltzmann models to analyze the wetting behavior and phase distribution in gas diffusion electrodes, aiding the design of more efficient metal-air batteries.

Contribution

It introduces heterogeneous wetting models for gas diffusion electrodes using FIB-SEM reconstructions and multiphase Lattice-Boltzmann simulations, enhancing understanding of phase behavior.

Findings

Pressure-saturation characteristics derived from simulations

Transport correlations for heterogeneous wetting behavior

Improved parameters for continuum-scale modeling

Abstract

Gas diffusion electrodes are commonly used in high energy density metal-air batteries for the supply of oxygen. Hydrophobic binder materials ensure the coexistence of gas and liquid phase in the pore network. The phase distribution has a strong influence on transport processes and electrochemical reactions. In this article we present 2D and 3D Rothman-Keller type multiphase Lattice-Boltzmann models which take into account the heterogeneous wetting behavior of gas diffusion electrodes. The simulations are performed on FIB-SEM 3D reconstructions of an Ag model electrode for predefined saturation of the pore space with the liquid phase. The resulting pressure-saturation characteristics and transport correlations are important input parameters for modeling approaches on the continuum scale and allow for an efficient development of improved gas diffusion electrodes.

Tables

Table 1. Table 1: Dimensions and resolution of the electrode sample in x , y , z 𝑥 𝑦 𝑧 x,y,z direction. The reconstructed electrodes are based on cubic voxels. For the LBM simulations the structure is mirrored in z 𝑧 z direction to increase the computational domain.

	Dimensions / µm	Voxels / -	Voxel size / nm
	x/y/z	x/y/z	x/y/z
FIB-SEM	27.8/23.9/10.0	504/346/84	55.25/69.1l/119
Reconstruction	27.8/23.9/10.0	504/432/180	55.25
LBM	27.8/23.9/20.0	126/108/90	221

Table 2. Table 2: Physical parameters of the system (SI units) and corresponding input parameters of the LBM simulations (lattice units).

	$\rho$		$\nu$		$\sigma$	$\Theta$
	Gas	Liquid	Gas	Liquid	Gas-Liquid	Ag	Binder
Physical	1.18	997	1.58$\cdot {10}^{-5}$	8.93$\cdot {10}^{-7}$	7.28$\cdot {10}^{-2}$	67 °[57]	140 °[45]
LBM	1	1	1/6	1/6	0.1	1.195/0.805†	0.357/1.643†
†  Density on wall nodes (gas/liquid) mu lu-3

Table 3. Table 3: Parameters of the Leverett function resulting from a fit to Eq. ( 26 ).

	$A$	$B$	$C$	$D$	$E$
configuration I	-0.813	0.233	5.079	2.323$\cdot$10-3	13.47
configuration II	4.557$\cdot$10-2	3.797$\cdot$10-2	8.138	0.01	10.98

Equations

$$\rho =\sum _i^{nv}{f}_i.$$

$$\mathbf{v}=\frac{1}{\rho }\sum _i^{nv}{\mathbf{e}}_i{f}_i.$$

$$p_k$$

$$f_i^k(\mathbf{x}+{\mathbf{e}}_i,t+\Delta t)=f_i^k(\mathbf{x},t)+{\Omega }_i^k\left(f_i^k(\mathbf{x},t)\right),$$

$${\Omega }_i^k={\Omega }_i^{k,1}\left({\Omega }_i^{k,3}+{\Omega }_i^{k,2}\right).$$

$$f_i^k(\mathbf{x},t^{\ast })=f_i^k(\mathbf{x},t)+{\Omega }_i^{k,1}\left(f_i^k(\mathbf{x},t)\right)$$

$${\Omega }_i^{k,1}=-\frac{f_i^k(\mathbf{x},t)-f_i^{k,\text{eq}}(\mathbf{x},t)}{{\tau }_k},$$

$${\tau }_k=\frac{3{\nu }_k}{c^2}+0.5\Delta t.$$

$$\bar{\nu }={\left(\sum _k\frac{{\rho }_k}{\rho {\nu }_k}\right)}^{-1}.$$

$$f_i^{k,\text{eq}}=\rho \left({\varphi }_i^k+w_i\left[3\frac{{\mathbf{e}}_i\mathbf{v}}{c^2}+\frac{9}{2}\frac{({\mathbf{e}}_i\mathbf{v}{)}^2}{c^4}-\frac{3}{2}\frac{{\mathbf{v}}^2}{c^2}\right]\right),$$

$$f_i^k(\mathbf{x},t^{\ast \ast })=f_i^k(\mathbf{x},t^{\ast })+{\Omega }_i^{k,2}\left(f_i^k(\mathbf{x},t^{\ast })\right)$$

$$\psi (\mathbf{x},t)=\frac{{\rho }_{\text{gas}}-{\rho }_{\text{liq}}}{{\rho }_{\text{gas}}+{\rho }_{\text{liq}}}.$$

$${\Omega }_i^{k,2}\left(f_i^k(\mathbf{x},t^{\ast })\right)=\frac{A_k}{2}\left|▽\psi \right|\left[w_i\frac{({\mathbf{e}}_i▽\psi {)}^2}{{\left|▽\psi \right|}^2}-B_i\right].$$

$$f_i^k(\mathbf{x},t^{\ast \ast \ast })={\Omega }_i^{k,3}\left(f_i^k(\mathbf{x},t^{\ast \ast })\right)$$

$${\Omega }_i^{\text{gas},3}\left(f_i^{\text{gas}}(\mathbf{x},t^{\ast \ast })\right)=\frac{{\rho }_{\text{gas}}}{\rho }{f}_i+\beta \frac{{\rho }_{\text{gas}}{\rho }_{\text{liq}}}{{\rho }^2}\text{cos}\left({\varphi }_i\right)f_i^{\text{eq}}{|}_{\mathbf{v}=0}$$

$${\Omega }_i^{\text{liq},3}\left(f_i^{\text{liq}}(\mathbf{x},t^{\ast \ast })\right)=\frac{{\rho }_{\text{liq}}}{\rho }{f}_i-\beta \frac{{\rho }_{\text{gas}}{\rho }_{\text{liq}}}{{\rho }^2}\text{cos}\left({\varphi }_i\right)f_i^{\text{eq}}{|}_{\mathbf{v}=0},$$

$$\text{cos}\left({\varphi }_i\right)=\frac{{\mathbf{e}}_i▽\psi }{\left|{\mathbf{e}}_i\right|\left|▽\psi \right|}.$$

$$f_i^k(\mathbf{x}+{\mathbf{e}}_i,t+1)=f_i^k(\mathbf{x},t^{\ast \ast \ast })$$

$$Bo=\frac{\text{gravitational forces}}{\text{surface forces}}=\frac{\rho g{l}^2}{\sigma }\sim 4\cdot 10^{-7}$$

$$Ca=\frac{\text{viscous forces}}{\text{surface forces}}=\frac{\mu \mathbf{v}}{\sigma }\sim 1\cdot 10^{-8}$$

$$Re=\frac{\text{inertial forces}}{\text{viscous forces}}=\frac{\rho \mathbf{v}l}{\mu }\sim 2\cdot 10^{-6}$$

$$\Delta p_c=\frac{\sigma }{r}\text{(2D)}\text{and}\Delta p_c=\frac{2\sigma }{r}\text{(3D)}$$

$$\Delta p_c=\frac{\sigma }{{\sigma }^{\text{LBM}}}\frac{1(lu)}{\Delta x}\Delta p_c^{\text{LBM}},$$

$$J(s)=\frac{\Delta p_c}{\sigma \left|cos(\Theta )\right|}\sqrt{\frac{B_0}{{\epsilon }_0}},$$

$$J(s)=\frac{\Delta p_c}{\sigma }\sqrt{\frac{B_0}{{\epsilon }_0}}.$$

$$D_k^{\text{eff}}/D_k^0={\left({\epsilon }_k\right)}^{\beta },$$

$$J(s)=A+B{e}^{C(s-0.5)}-D{e}^{-E(s-0.5)}$$

$$\mathbf{e}=c\scalebox{0.75}{\mbox{$\displaystyle\begin{pmatrix}0&1&0&-1&0&1&-1&1&-1\\ 0&0&1&0&-1&1&1&-1&-1\end{pmatrix}$}}$$

$$\mathbf{e}=c\scalebox{0.75}{\mbox{$\displaystyle\setcounter{MaxMatrixCols}{19}\begin{pmatrix}0&1&-1&0&0&0&0&1&-1&1&-1&1&-1&1&-1&0&0&0&0\\ 0&0&0&1&-1&0&0&1&-1&-1&1&0&0&0&0&1&-1&1&-1\\ 0&0&0&0&0&1&-1&0&0&0&0&1&-1&-1&1&1&-1&-1&1\end{pmatrix}$}}\;,$$

$$w_{i}=\left\{\begin{array}[]{l l}4/9&\quad i=0\\ 1/9&\quad i=1\dots 4\\ 1/36&\quad i=5\dots 8\end{array}\right.\;\;$$

Full text

††journal: Journal of Power Sources

url]http://www.dlr.de/tt/en/

1 Introduction

Metal-air batteries possess a very high theoretical energy density which makes them interesting for both mobile and stationary applications [1, 2]. At the negative electrode, metals like Al [3], Li [4], Mg [5], Na [6, 7], and Zn [8, 9] were suggested in the literature [10]. In recent years Li-air batteries received the most attention in the battery community [11, 12]. However, the only system which successfully reached the stage of mass production is the primary Zn-air battery. At the positive electrode oxygen is reduced and evolved during discharge and charge, respectively. Sufficient supply of O2 during discharge is accomplished by the concept of porous gas diffusion electrodes (GDEs). The electrodes are commonly made of carbon materials. However, in aqueous systems carbon is known to dissolve (’carbon corrosion’) and carbon-free GDEs were proposed [13, 14, 15]. Hydrophobic binder materials ensure the coexistence of gas and liquid phase in the porous structure of the GDE. The saturation behavior is characteristic for the porous material and can be described by capillary pressure saturation ($p_{c}-s$) curves. The amount and distribution of the liquid phase has a strong influence on transport processes. The transport in the gas phase ensures a good supply of O2 and, thus, allows to draw high current densities. Moreover, the binder improves the mechanical stability of the electrode.

Due to their application in alkaline fuel cells the concept of gas diffusion electrodes was studied already in the 1960s. Design optimizations of the electrodes were mainly done by intensive experimental studies. In recent years the improvements in computational efficiency made computer simulations a common design tool in the engineering disciplines. However, traditional CFD (computational fluid dynamics) tools like the volume of fluid (VOF) method [16, 17] have their limitations in the simulation of multiphase flow in complex geometries. In recent years the lattice Boltzmann method (LBM) [18, 19, 20] became increasingly popular for this class of problems because it is easy to implement and scales favorably. In LBM a probability distribution of discrete particle velocities is propagated on a computational lattice. Interactions between particles, boundaries, etc. are modeled by suitable collision operators. Several multiphase models were suggested in the literature for the simulation of immiscible fluids [21]. The most prominent ones are the Shan-Chen model [22, 23], the free-energy model [24, 25], and the color gradient or Rothman-Keller (RK) model [26, 27, 28]. A common problem of the methods is the numerical stability and accuracy in the simulation of systems with high density and viscosity ratios which also includes the air-water system. Recent publications present modifications of the models which are able to overcome this limitation [29, 30]. This extends the applicability of the method to technical systems like gas diffusion media of polymer electrolyte fuel cells and the number of publications in the field increased rapidly [31, 32, 33, 34, 35, 36]. However, to our knowledge this is the first publication of pore-scale LBM simulations of gas diffusion electrodes for metal-air batteries. In our work we use an RK type multiphase model to simulate multiphase flow in two and three dimensions. The model is based on the work of Leclaire [37, 38, 39] and Liu et al. [40]. In their publications the authors successfully demonstrated the simulation of high density ratios. In our study we focus on aqueous electrolytes, however, the presented methodology is not limited to this case. An important feature of our model is that we explicitly take into account the non-homogeneous wetting properties of the GDE which consists of hydrophilic electrode particles and hydrophobic binder fibers. This is an important step and has been barely pursued in previous pore-scale studies of electrochemical devices [34, 41, 42, 43, 44, 45].

The focus of our article is set on the development of a methodology for the characterization of gas diffusion electrodes for metal-air batteries. First, the structure of porous Ag model electrodes is reconstructed based on FIB-SEM tomography as explained in Section 2. The reconstruction serves as simulation domain for multiphase LBM simulations. Model equations and computational details are summarized in Sections 3 and 4. 2D and 3D simulations are performed to simulate the evolution of the phase distribution towards an energetic minimum (Section 5.1). The results are evaluated to obtain characteristic pressure-saturation curves (Section 5.2) and saturation dependent transport parameters (Section 5.3).

2 Electrode reconstruction

We use samples of Ag model electrodes which were characterized regarding their electrochemical performance in our previous article [46]. The focus of this work is set on the structural characterization and investigation of transport processes on the pore-scale. In this section we explain in detail the methodology which was developed for the reconstruction of gas diffusion electrodes using FIB-SEM tomography. The suggested procedure is schematically shown in Figure 1.

2.1 FIB-SEM

To obtain the micro-structure for our simulations a dual-beam ZEISS N-Vision 40 SEM-FIB instrument is used. The SEM has a field-emission electron gun with an acceleration voltage between 5 and 30 kV and a vertical electron-optic axis. The FIB is based on a Ga+ primary ion-beam with a 30 kV acceleration voltage. The ion-optic axis is at an angle of 54°  to the electron-optic axis. The standard vacuum levels for the electron-gun and the sample chamber are 10*-9* and 10*-5* mbar, respectively. The serial sectioning of the sample is achieved using FIB and the images are acquired using the SEM at specified milling intervals.

Before the SEM imaging of the sample, it must be ensured that the SEM and FIB images correspond to the same region-of-interest of the sample. This condition is obtained in the following way: Firstly, we set the sample stage to the eucentric tilt position to avoid an offset in sample position with sample tilt. Secondly, the co-incidence of the electron and ion beam has to be ensured. This is established by tilting the sample by 54°  with a constant working-distance of about 5 mm and adjusting the Z-axis until the SEM image of the sample comes into focus. A trench is cut into the sample in the vicinity of the region-of-interest with a relatively intense ion-beam current of 6.5 nA such that the cross-sectional (CS) plane becomes visible. The CS plane is then gently polished with an ion-beam current of 300 pA. Subsequently, the total region for 3D tomography is selected for ion-milling. The two samples were filled with a low viscosity epoxy resin in order to improve the contrast of the images. This is an important step to facilitate the reconstruction process. The left panel of Figure 1 shows a representative SEM image. We note that the resin has rather well impregnated the pore space of the GDEs. Horizontal and vertical dimension are in the following named $x$ and $y$, respectively. The direction perpendicular to the $x$-$y$ plane is denoted by $z$ and represents the direction of the FIB cut. The total thickness of the cut is 10 µm. In order to optimize the total duration of milling and image-acquisition, the SEM images are obtained at every fifth FIB-slice which results in 84 images representing slices of 0.12 µm thickness. Note, that ($x$, $y$) values of the voxel-size are related by the relation $y$ = $x$/sin(54°), where the value of $x$ is determined by the magnification of the image and the 54°  origin from the tilt of the sample.

2.2 Structure generation

The second panel of Figure 1 shows a representative binarized and cropped image of the electrode micro-structure (black). The images of different slices were aligned to account for the sample-drift during imaging. The perspective correction and pixel-size adjustments were done in the software package IMOD. The epoxy resin improves the contrast in the images and helps to identify solid particles. In spots where the impregnation of the electrode is incomplete the phases are assigned manually. Further details of the methodology of reconstruction is discussed in detail elsewhere [47]. Finally, the images were stacked to a virtual structure in the commercial software GeoDict [48]. The resulting geometry can be seen in the third panel of Figure 1.

The reconstructions are mirrored at the $x-y$ plane in order to increase the simulation domain in the direction of FIB sampling ($z$-direction). This step avoids systematic errors due to periodic boundary conditions which are used in the LBM simulations. The resulting geometry is subsequently coarsened in order to decrease the computational load. The dimensions of the samples are summarized in Table 1.

As outlined above the hydrophobic binder material is important for the distribution of the liquid phase and, thus, performance of the device. Unfortunately, the binder distribution is lost in the imaging and reconstruction process. In our reconstructions we distribute binder fibers by a random walk algorithm on the electrode surface. Binder fibers crossing the void space between electrode particles are neglected. The resulting binder distribution is shown in red color in the right panel of Figure 1. It is in reasonable optical agreement with distributions observed in SEM images [46].

3 Lattice Boltzmann Method

3.1 Model description

In this work we perform two (D2Q9) and three (D3Q19) dimensional simulations of two-phase flow in porous gas diffusion electrodes. The LBM multiphase models employed in this study are based on the color-gradient approach and follow the recent publications of Liu and Leclaire et al. [37, 38, 39, 40]. A detailed derivation and discussion of the model equations can be found in their publications.

In RK type models the state of each phase $k$ ($k$=gas, liquid) is described with the help of a probability distribution function $f^{k}_{i}(\mathbf{x},t)$, where $i$ is the index representing discrete directions in velocity space (see Eq. (27) and (28)). The macroscopic properties at a lattice node $\mathbf{x}$ are given by the moments of the probability distribution. For instance, the density $\rho$ is calculated by the 0th moment which is simply the sum of the probabilities in all lattice directions $nv$

[TABLE]

The fluid velocity $\mathbf{v}$ follows as the 1st moment which can be regarded as average of the discrete velocities $\mathbf{e}$ weighted by their probability

[TABLE]

Finally, the pressure can be calculated according to

[TABLE]

where $\alpha_{k}$ is a simulation parameter given in A.

In general the temporal evolution of the probability distribution $f^{k}_{i}(\mathbf{x},t)$ is described by

[TABLE]

where the collision operator $\Omega_{i}^{k}$ is in RK type models a result of three sub operators [49, 50]

[TABLE]

The operators are applied successively to the probability distribution $f^{k}_{i}(\mathbf{x},t)$ and are defined by

Single-phase collision (SRT-BGK)

[TABLE]

In our model we apply the standard SRT-BGK approximation [51] to model fluid interactions in the same phase

[TABLE]

where $\tau_{k}$ is the relaxation time of the collision process and $f_{i}^{k,\text{eq}}$ the local equilibrium distribution of fluid velocities. The relaxation time $\tau_{k}$ can be related to the kinematic viscosity $\nu_{k}$ by

[TABLE]

In our simulations we use a single relaxation time $\tau_{\text{gas}}=\tau_{\text{liq}}=\bar{\tau}$ for both phases which is calculated from a density average of the kinematic viscosity

[TABLE]

The local equilibrium distribution $f^{k,\text{eq}}$ is given by the Maxwell distributions

[TABLE]

where $\phi^{k}_{i}$ (Eq. (31) and (32)) determines the compressibility of the fluid and $w_{i}$ is a lattice specific weighting parameter (see Eq. (29) and (30)). 2. 2.

Two-phase collision (Perturbation)

[TABLE]

The effect of surface tension between the two phases is modeled by the perturbation operator. In order to distinguish the phases it is convenient to introduce a color-field $\psi(\mathbf{x},t)$

[TABLE]

Regions of gas and liquid phase are marked by values of $\psi$ close to 1 and -1, respectively. The perturbation operator takes the form [50, 40]

[TABLE]

$A_{k}$ is a parameter controlling the surface tension (see Eq. (34)), $\triangledown\psi$ is the color gradient in the two-phase region, and $B_{i}$ is a parameter which has to be chosen in order to recover the Navier-Stokes Equations [50, 40] (see Eq. (36) and (37)). The color gradient determining the surface normal in the two-phase region is approximated by higher-order isotropic discretization schemes which significantly reduce spurious velocities [52, 40, 37]. 3. 3.

Two-phase collision (Recoloring)

[TABLE]

The perturbation operator models the effect of surface tension, however, it does not enforce phase separation. This is ensured by the recoloring operator

[TABLE]

where $\rho=\rho_{\text{gas}}+\rho_{\text{liq}}$ is the total density, $f_{i}=f_{i}^{\text{gas}}+f_{i}^{\text{liq}}$ the total probability distribution, and $\phi_{i}$ the angle between the color gradient and the lattice direction vector

[TABLE]

The operator allows a moderate mixing of the two phases at the interface which additionally reduces spurious currents. The interface thickness is controlled by the parameter $\beta$. Moreover, it was shown that the proposed operator circumvents the problem of ’lattice-pinning’ [53]. 4. 4.

Streaming

[TABLE]

Finally, the modified probability distributions are streamed to their new positions on the computational grid and the corresponding boundary conditions are applied (see Section 3.2).

3.2 Boundary conditions

There are two types of boundaries in micro-structure resolved simulations of porous media: boundaries of the computational domain and boundaries at the fluid-solid interface. For the first type we employ simple periodic boundary conditions. This choice has implications on the simulation procedure and will be discussed in Section 4.1.

At the fluid-solid interface we use a simple bounce back scheme. The contact angle is adjusted by assigning an effective density to the probability distribution at solid nodes [54]. This has an effect on the local color-gradient and determines if the surface is wetting or non-wetting. In our model the densities on solid nodes were adjusted to match the contact angles which are observed experimentally (see Figure 2). Parameters can be found in Table 2.

3.3 Model parameterization

In this study we parameterize our models for the air-water system in order to represent metal-air batteries with aqueous electrolytes. It has to be noted that other electrolyte systems using organic solvents or ionic liquids can be treated in the same framework. An analysis of dimensionless numbers is helpful to evaluate the relevant forces of the physical problem. Fluid flow in porous media is mainly determined by gravitational, viscous, inertial, and surface forces. The dimensionless numbers describing the ratio between these forces are

Bond number

[TABLE]

• 2.

Capillary number

[TABLE]

• 3.

Reynolds number

[TABLE]

The characteristic length $l$ is in our studies given by the mean pore diameter $d_{50}=0.87$ µm of the electrode. The flow velocity $\mathbf{v}$ can be approximated by the flow rate of volume replacement experiments and is typically less than the $1\cdot 10^{-6}$ m s*-1* used in above calculations [55]. Moreover, in the limit of fluid mechanical equilibrium it will tend to zero. The small values of the dimensionless numbers demonstrate that gravitational and viscous forces are negligible compared to the strong influence of surface forces. Therefore, we use unity density and viscosity ratios in our simulations to improve the numerical stability. In a post processing step suitable scaling rules are applied to the simulation results [56]. The Laplace equation

[TABLE]

allows to deduce a relationship between the capillary pressure of the ’physical’ ($\Delta p_{c}$) and simulated system ($\Delta p^{\text{LBM}}_{c}$) according to

[TABLE]

where $lu$ is one lattice unit, $\Delta x$ the size of one voxel and $\sigma$ and $\sigma^{\text{LBM}}$ the surface tension of the physical system and LBM simulation, respectively. The parameters of this study are summarized in Table 2.

3.4 Model validation

In order to validate our model we present two simple test cases.

Bubble test

The pressure difference across the interface of a steady bubble can be calculated by the Laplace equation (Eq. (21)). In order to validate our model we perform 2D simulations on a lattice with 100x100 nodes and varying bubble diameter. Figure 2 a) shows simulation results for diameters ranging from 10 to 90 lattice units. The size of one lattice unit corresponds to the voxel size of the reconstruction (see Table 1). The simulated pressures reproduce favorably the Laplace equation. At small diameters a minor deviation of simulation results is observed. This can be assigned to the increasing ratio of surface to bulk voxels which is important to notice for the calculation of $p_{c}-s$ curves.

Static contact angle

The second case is the simulation of static contact angles as presented in Figure 2 b). The fluid density of solid nodes was adjusted to reproduce the contact angle of water on Ag [57] and PTFE surfaces [45]. The Figure demonstrates that our model is able to simulate the wetting characteristics of the Ag GDE sample.

4 Simulation methodology

4.1 Initial conditions

The $p_{c}-s$ curves of porous media are commonly determined with the method of standard porosimetry (MSP) [58, 59] or dynamic volume replacement experiments [60, 61, 62, 55]. After filling the porous medium to a defined saturation level the corresponding average capillary pressure follows as difference between the pressure in the gas and the liquid phase. It was found that the direction of the process (injection/imbibition or removal/drainage of the liquid phase) causes a hysteresis in the resulting $p_{c}-s$ curves. It is therefore expected that a different protocol for loading the porous medium with liquid or gas phase will lead to slightly different $p_{c}-s$ curves. Since the conditions for the standard dynamic volume experiments do not match the conditions in the GDE under dynamic operations, we chose a different simulation set up to determine $p_{c}-s$ curves. This setup corresponds to establishing a force balance in the volume of the porous medium after initializing the simulation with a maximally wetting or non-wetting condition respectively for a given volume fraction of liquid in the GDE. We argue, that the thus obtained pressure saturation curves are closer to the real operating condition of a GDE in a metal-air cell, where electrolyte is pushed out of the GDE due to the occurrence of solid reaction products in the bulk of the porous medium and not due to the application of external pressure forces [46]. The different initial configurations within our LBM model with periodic boundary conditions are shown in the first row of Figure 3. During drainage the electrolyte will remain preferentially on the hydrophilic electrode particles. We maximize the interface between liquid phase and hydrophilic Ag particles by randomly placing liquid droplets on the electrode surface or, if all surface sites are occupied, in contact with another liquid droplet (configuration I). In the second configuration the liquid phase is introduced in the pore space as one solid block (Figure 3 right column - configuration II) in order to minimize the interfacial area between gas and liquid phase. This situation is comparable to the imbibition process in the experimental setup where the electrolyte is forced into the porous structure and also occupies areas which are energetically not favorable. Our simulations are able to reproduce the hysteretic behavior which is also observed in the forced drainage-imbibition experiments. These results will be discussed in more detail in Section 5.2.

4.2 Pressure-saturation curves

For the determination of $p_{c}-s$ curves we perform independent 2D and 3D simulations at various saturation levels of the GDE sample. The number of iterations was chosen sufficiently large to achieve fluid mechanical equilibrium. At the end of the simulation the capillary pressure follows as $\Delta p_{c}^{\text{LBM}}=\bar{p}^{\text{liq}}-\bar{p}^{\text{gas}}$, where $\bar{p}^{k}$ is averaged over the simulation domain. In the case of 2D simulations we conduct at each saturation 10 independent simulations on randomly chosen two-dimensional slices of the reconstruction in order to get statistically significant results.

Pressure saturation curves are commonly described with the so called Leverett function [63, 64]

[TABLE]

where $B_{0}$ is the permeability and $\varepsilon_{0}$ the porosity of the electrode. In this study we use a modification of the standard formulation according to Hao et al. [34] to account for the heterogeneous wetting properties of the GDE

[TABLE]

4.3 Transport parameters

The final phase distribution at the end of the 3D simulations is used as input for the calculation of saturation-dependent effective transport parameters and surface areas in the software GeoDict [48]. The final fluid distribution is basically regarded as ’frozen’ in the porous structure. For the determination of effective transport parameters in the gas phase the voxels of the liquid phase are then treated as solids and vice versa. In their study Garcia et al. [65] took a similar approach to determine the diffusivity of partially saturated GDLs. However, they obtained the distribution of the liquid phase from tomographic images instead of two-phase simulations.

In the modeling of electrochemical devices the Bruggeman correlation is often used for the approximation of effective transport parameters. The diffusivity which is the ratio between effective and bulk diffusion coefficients is in this approach described by

[TABLE]

where $\varepsilon_{k}$ is the volume fraction of the transporting phase and $\beta$ the so-called Bruggeman coefficient. The volume fraction of the liquid and gas phase are calculated according to $\varepsilon_{\text{liq}}=\varepsilon_{0}s$ and $\varepsilon_{\text{gas}}=\varepsilon_{0}(1-s)$, respectively.

5 Results and discussion

5.1 Electrolyte distribution

Figure 3 shows the phase distribution during 3D simulations of configuration I (left) and configuration II (right) at a saturation of the pore space with liquid electrolyte of 50%. The first row illustrates the initial conditions which were applied to mimic the drainage and imbibition of the liquid phase. (cf. Section 4.1). In configuration I the electrolyte is initially distributed on the electrode surface. During the simulation the distribution evolves towards a local minimum of the free energy. Our simulations reach a quasi stationary state after about 1$\cdot$104 iterations. Further changes are only marginal. This is also reflected in the corresponding pressure signal which is shown as function of iterations in Figure 4 b). A visually similar steady-state is also observed for configuration II. This indicates that the initial condition which we apply in our simulations has only a minor influence on the saturation dependence of effective transport parameters. This effect will be discussed in Section 5.3.

5.2 Pressure-saturation curves

Figure 4 shows the capillary pressure signal during 2D and 3D simulations representing the drainage process. In both cases the pressure signal approaches a constant value at the end of the simulations which indicates that the system is in a local minimum of the free energy. The 3D simulations give a clear trend of increasing capillary pressure with increasing liquid phase saturation. It has to be noted that the 2D simulations shown in Figure 4 a) were not performed on the same slice of the reconstruction. Therefore, the fluctuations in pressure which occur are due to different simulation domains and a clear trend can not be deduced from the graph. In our approach we take the average of 10 simulations to capture this effect statistically. Figure 5 shows $p_{c}-s$ curves for configuration I (left) and configuration II (right). The results of the 2D simulations are shown in red color and the error bars represent the corresponding standard deviation. The error bars are larger at low and high saturation which can be explained by a reduced configurational freedom. Meaning, reconfigurations of the phases are suppressed and they remain in their initial (random) configuration. The calculated pressures of 3D simulations are included as blue triangles. There is a systematic deviation at high saturation, however, the general agreement between 2D and 3D simulations is favorable. This is an important result which shows that a series of computationally efficient 2D simulations can be used for a screening of different electrode structures. In a previous publication [46] we predicted a liquid phase saturation of $\approx$ 50 % at $\Delta p=0$ by a fit of continuum simulations to electrochemical data. This is in excellent agreement with the results presented in this work and indicates the validity of our approach.

Pressure saturation curves are commonly presented as dimensionless Leverett functions (see Eq. (24)). Hao et al. showed that an expression of the form

[TABLE]

is able to give a good representation of $p_{c}-s$ data. We fit the parameters of Eq. (26) to the $p_{c}-s$ curves of our 2D simulations. Figure 6 shows Leverett functions of configuration I (left) and configuration II (right). Similar to the experimental work on GDLs we observe a hysteresis at intermediate saturation [60, 55]. This can be explained by looking at the experimental procedure and taking into account micro-structural effects in the GDE. Some areas are connected to the pore network only through narrow pores. According to the Laplace equation they are filled with the liquid phase only at high capillary pressures. On the reverse process a higher negative pressure difference is needed to withdraw the electrolyte from this part of the pore network. A similar reasoning can be made for the simulation approach presented in this study. In the simulations starting with configuration I all parts of the pore network are accessible for the liquid phase. This results in a lower capillary pressure compared to simulations with initial configuration II where the access is restricted to surrounding pores. This result justifies our choice of initial conditions to mimic the two processes. However, it has to be noted that the hysteresis is within the standard deviation of the simulations. Experimental studies on the electrodes will be needed to validate our simulation methodology. The parameters of Eq. (26) resulting from the fit to the simulation data are summarized in Table 3.

5.3 Effective transport parameters

Figure 7 shows the diffusivity in the gas and liquid phase for configuration I (left) and configuration II (right). Results are presented individually for the three spatial coordinates $x$, $y$, and $z$. The graphs show that the transport behavior is anisotropic. Especially, in $z$-direction which is the direction of FIB sampling the transport parameters are smaller. This can be attributed to a lower spatial resolution and indicates that the anisotropy might be an artifact of the reconstruction. The overall trend of the diffusivity with saturation is the same for both initial configurations. The diffusivity of the gas phase decreases with an increasing saturation of the pore space with the liquid phase. This corresponds to the behavior predicted by the Bruggeman correlation (cf. Eq. (25)) presented in Figure 8 a).

The main transport direction in the GDE during operation is in our notation denoted by $y$ and we limit our discussion in the remainder of this section to this most important case. Figure 8 a) shows the calculated diffusivity in the gas and liquid phase for simulations with initial configuration I and II. As mentioned in Section 5.1 the final liquid phase distribution is similar for both cases and, thus, also saturation dependent transport coefficients are comparable. The dotted line represents the Bruggeman correlation with a coefficient of $\beta$=1.7. The graph demonstrates that the correlation gives a reasonable representation of the simulated values although the coefficient is slightly higher than the standard value of 1.5 which is usually assumed in the literature.

An important parameter for the performance of the GDE is the length of the triple phase boundary (TPB). Here, gas, liquid, and solid phase are in close contact and reaction kinetics are facile. The dependence of the TPB on the saturation is shown in Figure 8 b) for both, configuration I and configuration II. The length of the TPB first increases with saturation. Then, it reaches a plateau because large parts of the pore space are completely filled with electrolyte. The same effect leads to a decrease in TPB length at high saturation. This shows that at a high to moderate saturation the performance of the GDE is improved. The deviation between configuration I and configuration II origins from the more homogeneous distribution of the liquid phase in configuration I.

The correlations determined above are important input for simulations on the continuum scale [66, 46] and allow an efficient improvement of gas diffusion electrodes. An overview of this multi-scale approach can also be found in Ref. [67].

6 Conclusions

In this article we present a methodology for the characterization of gas diffusion electrodes for metal-air batteries. First, we take FIB-SEM images of an Ag model electrode. In a subsequent step the images are binarized and stacked to a virtual reconstruction of the electrode. On this geometry we performed 2D and 3D LBM simulations with a RK type multiphase model. To the best of our knowledge this study is the first publication investigating the complex multiphase transport mechanism in GDEs using LBM on the pore-scale. An important feature of our model is that we take into account the non-homogeneous wetting behavior of hydrophilic Ag particles and hydrophobic binder. The model is parametrized to represent an aqueous electrolyte system, however, other electrolytes (e.g. organic solutions or ionic liquids) can be treated in the same framework. Results of the simulations are evaluated for the determination of $p_{c}-s$ curves and saturation dependent transport parameters. Our simulations demonstrate that a series of 2D simulations is an efficient tool for the screening of $p_{c}-s$ characteristics of GDEs. Transport parameters and specific surface areas are determined from the final phase distributions of 3D simulations. The diffusivity generally follows the trends predicted by the Bruggeman correlation with a coefficient of 1.7. The results of this study provide important input parameters for simulations on the continuum scale. Moreover, the methodology presented here can be used as a tool for the optimization of GDEs for metal air batteries.

Acknowledgments Funding from the European Union’s innovation programme Horizon 2020 (2014-2020) under grant agreement n°646186 (ZAS) is gratefully acknowledged. The authors would also like to thank Norbert Wagner, Institute for Engineering Thermodynamics, DLR Stuttgart for providing the samples of Ag model electrodes and Prasanth Balasubramanian, Ulm University for help in tomogram reconstruction. Santhana Eswara would like to thank Jörg Bernhard and Prof. Ute Kaiser of Ulm University for fruitful discussions. Timo Danner would like to thank Prof. Wolfgang Bessler, University of Applied Sciences Offenburg, for discussions and comments.

Appendix A LBM parameters

A.1 Computational lattice

The D2Q9 and D3Q19 lattices [68] in cartesian coordinates are given by

[TABLE]

and

[TABLE]

respectively.

A.2 Single-phase collision

For the D2Q9 and D3Q19 lattices $w_{i}$ is given by

[TABLE]

and

[TABLE]

respectively.

The parameter $\phi^{k}_{i}$ is related to the compressibility of the fluid and, thus, to the speed of sound $c^{k}_{s}$ and hydrostatic pressure in phase $k$. For the D2Q9 and D3Q19 lattice $\phi^{k}_{i}$ is given by [39, 40]

[TABLE]

and

[TABLE]

respectively, where one of the $\alpha_{k}$ is a free parameter setting the pressure level in the system. The values of $\alpha_{k}$ are related by

[TABLE]

in order to guarantee a stable interface ($p_{\text{gas}}=p_{\text{liq}}$). In our simulations we set $\alpha_{\text{gas}}$ to 4/9.

A.3 Two-phase collision: Perturbation

The parameter $A$ controlling the surface tension is defined as [40]

[TABLE]

where $\sigma^{\text{LBM}}$ is the surface tension and $\bar{\tau}$ the density averaged relaxation time. In our simulations we use a fourth-order isotropic approximation for the calculation of the color gradient [40]

[TABLE]

In the literature approximations of higher order are reported [52], however, the computational load increases significantly.

The parameters $B_{i}$ depend on the lattice and have to be chosen in order to ensure the conservation of mass in the perturbation step. For the D2Q9 and D3Q19 lattice the parameters are given by [50, 40]

[TABLE]

and

[TABLE]

respectively.

A.4 Two-phase collision: Recoloring

The only additional parameter appearing in the recoloring operator is $\beta=0.85$ which determines the thickness of the inter-facial layer.