Mesoscopic Modeling of Multiphysicochemical Transport Phenomena in Porous Media

Abstract

We present our recent progress on mesoscopic modeling of multiphysicochemical transport phenomena in porous media based on the lattice Boltzmann method. Simulation examples include injection of CO₂-saturated brine into a limestone rock, two-phase behavior and flooding phenomena in polymer electrolyte fuel cells, and electroosmosis in homogeneously charged porous media. It is shown that the lattice Boltzmann method can account for multiple, coupled physicochemical processes in these systems and can shed some light on the underlying physics occurring at the fundamental scale. Therefore, it can be a potential powerful numerical tool to analyze multiphysicochemical processes in various energy, earth, and environmental systems.

1. Introduction

Multiphysicochemical transport phenomena in porous media are ubiquitous, particularly in various energy, earth, and environment systems. One example is the disposal of supercritical CO₂ in geologic formations, the most promising near-term solution to the problem of reducing carbon emissions into the atmosphere [1]. Experimental analyses of the long-term fate of CO₂ after injection into the geologic formations are not possible with relatively short-term laboratory experiments. Therefore it is necessary to employ comprehensive numerical models that incorporate multiple physicochemical processes involving interactions between the injected CO₂, the brine in the pore spaces, and the minerals lining the pores. Supercritical CO₂, as a buoyant and slightly miscible fluid, once injected, displaces brine from the pore space in a complex pattern. At the interface with brine, CO₂ dissolves into the brine to form carbonic acid that can react with and dissolve minerals eventually leading to mineral precipitation further along the flow path. Clearly, there are multiple physics processes involved, including hydrodynamics, thermodynamics, chemical dynamics, and electrodynamics (because the surface of most natural media is charged). All these processes are ultimately governed by pore-scale interfacial phenomena, which occur at scales of microns. However, because of the wide disparity in scales ranging from pore to field, a continuum formulation based on spatial averages taken over length scales much larger than typical pore and mineral grain sizes is often utilized, implying the existence of a representative elemental volume (REV) [2]. As a result, spatial heterogeneities at smaller scales are unresolved and the aggregate effects of the pore-scale (mesoscopic scale) processes are taken into account through various effective constitutive parameters. One of the goals of performing pore-scale simulations is to obtain values for these constitutive parameters through upscaling the pore-scale results. Other goals are to identify key parameters and physicochemical processes that control macroscopic phenomena, and to validate continuum descriptions.

Another example is fuel cells, and in particular polymer electrolyte fuel cells (PEFCs). In PEFCs, the catalyst layer (CL) is the host to several competing transport mechanisms involving charge (proton and electron), species (oxygen, nitrogen, and water vapor), and liquid water transport. The multi-faceted functionality of a gas diffusion layer (GDL) includes reactant distribution, liquid water transport, electron transport, heat conduction, and mechanical support to the membrane-electrode-assembly. Despite tremendous recent progress in enhancing the overall cell performance, a pivotal performance limitation in PEFCs is manifested in terms of mass transport loss originating from suboptimal liquid water transport and resulting flooding in the constituent components [3]. In recent years, several macroscopic computational models for multiple-physicochemical transport processes in PEFCs [4–10] have been developed. These macroscopic models again are based on the theory of volume averaging and treat the catalyst layer and gas diffusion layer as macrohomogeneous porous layers. Due to their macroscopic nature, the current models fail to resolve the influence of the structural morphology of the CL and GDL on the underlying physics. Mesoscopic modeling is critical to understanding of the overall structure-wettability-transport interactions as well as the underlying multiphysicochemical processes in the CL and GDL, and hence is a useful tool for design and optimization of microstructures of electrodes for better performance and durability.

In this paper, we review our recent work on mesoscopic modeling of multiphysicochemical processes in porous media, based on the lattice Boltzmann method (LBM), a relatively new numerical method for simulating fluid flows and modeling physics in fluids [11]. Originating from the classical statistical physics, LBM is a mesoscopic method based on simplified kinetic equations. In the LBM, the fluid is modeled as a collection of fictitious particles propagating and colliding over a discrete lattice domain. Mesoscopic continuity and momentum equations can be obtained from this propagation-collision dynamics through a rigors mathematical analysis. The particulate nature and local dynamics provide advantages for complex boundaries and parallel computation. In addition, the kinetic nature of the LBM makes it easy to account for new physics in the LBM framework, which is especially useful for modeling multiphysicochemical phenomena. In Section 2, the partial differential equations governing fluid flow, transport of reactive species and electric potential, as well as mineral reactions in porous media will be given. In Section 3, the implementation of the LBM to solve these governing equations will be presented. Some simulation examples will be given in Section 4 and concluding remarks in Section 5.

2. Governing Equations

Consider the electrolyte fluid flowing through solid porous media. Although the pore scale in this study may be at micrometer scale, the fluid can still be treated as a continuum Newtonian fluid since the characteristic size is much larger than the molecular diameters [12–14].

2.1. Continuity and Momentum Equations

For isothermal incompressible fluid flow, the continuity and momentum equations can be written as [15]

$\begin{array}{c}\nabla ·\mathbf{u}=0,\\ \rho \frac{\partial \mathbf{u}}{\partial t}+\rho \mathbf{u}·\nabla \mathbf{u}=-\nabla p+\rho \nu {\nabla }^2\mathbf{u}+\mathbf{F},\end{array}$

(1)

where ρ represents the density of the fluid, t the time, u the velocity vector, p the pressure, ν the kinetic viscosity, and F the body force density which may include all the effective body forces.

For the hydrodynamic boundary condition, we use the nonslip model at the solid surfaces. The slip boundary conditions have been adopted in some recent studies [16, 17]; however, a careful molecular study showed that the hydrodynamic boundary condition, slip or not, depended on the molecular interactions between fluid and solid and on the channel size [18, 19]. For the flow in porous media considered in this work, the nonslip boundary condition is still valid.

2.2. Transport Equations for Aqueous Species and Electrical Potential

For the ith ion species in the solute, the mass conservation equation describing transport and reaction can be written in the general form [20]

$\frac{\partial C_i}{\partial t}+\nabla ·{\mathbf{J}}_i+{\lambda }_i{C}_i=R_i,$

(2)

where C _i denotes the ionic concentration, J _i the species flux, λ _i a radioactive decay constant, and R _i the rate at which the ith species is produced or consumed by chemical reactions. The flux J _i , consisting of contributions from advection, diffusion, and electrochemical migration terms has the form [20]

${\mathbf{J}}_i=-\frac{e{z}_i{D}_i}{kT}{C}_i\nabla \mathrm{\Psi }-D_i\left(\nabla C_i+C_i\nabla \ln {\gamma }_i\right)+C_i\mathbf{u},$

(3)

where the first term on the right refers to electrochemical migration, the second term to aqueous diffusion, and the last term to advective transport. Here z _i , D _i , and γ _i denote the ion algebraic valence, the diffusivity, and the activity coefficient of the ith species, respectively; and e, k, and T denote the absolute charge of electron, the Boltzmann constant, and the absolute temperature, respectively. The quantity Ψ represents the local electrical potential caused by the ionic distribution which is governed by the Poisson equation

$\nabla ·\left({\epsilon }_r{\epsilon }_0\nabla \mathrm{\Psi }\right)=-{\rho }_e=-{\displaystyle \sum }_{i}e{z}_i{C}_i,$

(4)

where ε _r is the local dimensionless fluid dielectric constant, ε₀ the permittivity of a vacuum, and ρ _e the net charge density. Assuming no radiation and constant activity coefficient and substituting (4) into (2), we have

$\frac{\partial C_i}{\partial t}+\mathbf{u}·\nabla C_i=D_i{\nabla }^2{C}_i+\frac{e{z}_i{D}_i}{kT}\nabla ·\left(C_i\nabla \mathrm{\Psi }\right).$

(5)

This is the Nernst-Planck equation [21], where F can be any kind of effective body force. In this contribution we only consider the static electrical force from an external electric field. The general form of electrical force in electrokinetic fluids can be expressed as

${\mathbf{F}}_E=-{\rho }_e\nabla {\mathrm{\Psi }}_{\text{ext}},$

(6)

where Ψ_ext is the external electrical potential field.

When the ionic convection is negligible and the electric potential is continuously derivable, (5) has a simple steady-state solution for dilute electrolyte solutions:

$C_i=C_{i,\infty }{e}^{-e{z}_i\mathrm{\Psi }/kT}.$

(7)

Substituting (7) into (4) yields the nonlinear PB equation [22]

${\nabla }^2\mathrm{\Psi }=-\frac{1}{{\epsilon }_r{\epsilon }_0}{\displaystyle \sum }_{i}e{z}_i{C}_{i,\infty }\exp \left(-\frac{e{z}_i}{kT}\mathrm{\Psi }\right).$

(8)

2.3. Equations for Mineral Reaction Rates

Heterogeneous reactions between aqueous species and minerals at the pore-mineral interface are given by [23]

$AD{\frac{\partial {\mathrm{\Psi }}_j}{\partial n}|}_s={\displaystyle \sum }_s^{}{\nu }_{js}{A}_s{\stackrel{\circ }{I}}_s^{}\left(s\right),$

(9)

where n denotes the unit normal perpendicular to the fluid-mineral interface pointing toward the fluid phase, D denotes the aqueous diffusion coefficient assumed to be the same for all species for simplicity, and ${\stackrel{}{I}}_s^{}(s)$ denotes the reaction flux of the sth mineral at its surface, taken as positive for precipitation and negative for dissolution. The total surface area A across which diffusion takes place equal to the sum of individual mineral surface areas A _s is given as

$A={\displaystyle \sum }_s^{}{A}_s.$

(10)

3. Lattice Boltzmann Model Implementation

3.1. Incompressible Lattice Boltzmann Model for Single-Phase Flow

In order to eliminate compressible effects in the conventional LBM, here we use an incompressible LB model constructed by Guo et al. [24]. The pore-scale flow of a single aqueous fluid phase is simulated by the following evolution equation:

$f_{\alpha }\left(\mathbf{x}+{\mathbf{e}}_{\alpha }{\delta }_t,t+{\delta }_t\right)=f_{\alpha }\left(\mathbf{x},t\right)-\frac{f_{\alpha }\left(\mathbf{x},t\right)-{f_{\alpha }}^{\text{eq}}\left(\mathbf{x},t\right)}{\tau }.$

(11)

In the above equation, δ _t is the time increment, f _α the distribution function along the α direction in velocity space, ${f_{\alpha }}^{\text{eq}}$ the corresponding equilibrium distribution function, and τ the dimensionless relaxation time. For the commonly used two-dimensional, nine-speed LB model (D2Q9), the discrete velocities e _α have the following form:

${\mathbf{e}}_{\alpha }=\{\begin{array}{cc}\left(\mathrm{0,0}\right),\hfill & \alpha =0,\hfill \\ \left(\cos {\theta }_{\alpha },\sin {\theta }_{\alpha }\right)c,\hfill & {\theta }_{\alpha }=\frac{\left(\alpha -1\right)\pi }{2},\alpha =1-4,\hfill \\ \sqrt{2}\left(\cos {\theta }_{\alpha },\sin {\theta }_{\alpha }\right)c,\hfill & {\theta }_{\alpha }=\frac{\left(\alpha -5\right)\pi }{2}+\frac{\pi }{4},\alpha =5-8.\hfill \end{array}$

(12)

For the incompressible LB model, the equilibrium distribution is defined by Guo et al. [24]:

${f_{\alpha }}^{\text{eq}}=\{\begin{array}{cc}-4\sigma \frac{p}{\rho c^2}+s_{\alpha }\left(\mathbf{u}\right),\hfill & \alpha =0,\hfill \\ \lambda \frac{p}{\rho c^2}+s_{\alpha }\left(\mathbf{u}\right),\hfill & \alpha =1-4,\hfill \\ \gamma \frac{p}{\rho c^2}+s_{\alpha }\left(\mathbf{u}\right),\hfill & \alpha =5-8,\hfill \end{array}$

(13)

where σ,λ, and γ are the parameters satisfying

$\begin{array}{c}\lambda +\gamma =\sigma ,\\ \lambda +2\gamma =\frac{1}{2},\\ s_{\alpha }\left(\mathbf{u}\right)={\omega }_{\alpha }\left[3\frac{{\mathbf{e}}_{\alpha }·\mathbf{u}}{c^2}+9\frac{{({\mathbf{e}}_{\alpha }·\mathbf{u})}^2}{2c^4}-\frac{3{\mathbf{u}}^2}{2c^2}\right].\end{array}$

(14)

In these equations, c = δx/δt, where δx is the space increment, and p and u are the pressure and velocity of the fluid, respectively. The corresponding weight coefficients are

${\omega }_{\alpha }=\{\begin{array}{cc}\frac{4}{9},\hfill & \alpha =0,\hfill \\ \frac{1}{9},\hfill & \alpha =1-4,\hfill \\ \frac{1}{36},\hfill & \alpha =5-8.\hfill \end{array}$

(15)

Equation (11) has been shown to recover (1) [24], with the velocity and pressure given by

$\begin{array}{c}\mathbf{u}={\displaystyle \sum }_{\alpha =1}^8{\mathbf{e}}_{\alpha }{f}_{\alpha }\end{array}$

(16)

$\begin{array}{c}\frac{p}{\rho }=\frac{c_{}^2}{4\sigma }\left[{\displaystyle \sum }_{\alpha =1}^8{f}_{\alpha }+s_0\left(\mathbf{u}\right)\right],\end{array}$

(17)

respectively.

3.2. Lattice Boltzmann Model for Two-Phase Flow

The interaction-potential model, originally proposed by Shan and Chen [25], and henceforth referred to as the S-C model, introduces k distribution functions for a fluid mixture comprising of k components. Each distribution function represents a fluid component and satisfies the evolution equation. The nonlocal interaction between particles at neighboring lattice sites is included in the kinetics through a set of potentials. The LB equation for the kth component can be written as

$f_i^k\left(\mathbf{x}+{\mathbf{e}}_i{\delta }_t,t+{\delta }_t\right)-f_i^k\left(\mathbf{x},t\right)=-\frac{f_i^k\left(\mathbf{x},t\right)-f_i^{k\left(\text{eq}\right)}\left(\mathbf{x},t\right)}{{\tau }_k},$

(18)

where $f_i^k(\mathbf{x},t)$ is the number-density distribution function for the kth component in the ith velocity direction at position x and time t, and δ _t is the time increment. In the term on the right-hand side, τ _k is the relaxation time of the kth component in lattice unit, and $f_i^{k(\text{eq})}(\mathbf{x},t)$ is the corresponding equilibrium distribution function.

The phase separation between different fluid phases, the wettability of a particular fluid phase to the solid, and the body force are taken into account by modifying the velocity used to calculate the equilibrium distribution function. An extra component-specific velocity due to interparticle interaction is added on top of a common velocity for each component. Interparticle interaction is realized through the total force, F _k , acting on the kth component, including fluid/fluid interaction, fluid/solid interaction, and external force. More details can be found in [26, 27].

The continuity and momentum equations can be obtained for the fluid mixture as a single fluid using Chapman-Enskog expansion procedure in the nearly incompressible limit:

$\begin{array}{c}\frac{\partial \rho }{\partial t}+\nabla ·\left(\rho \mathbf{u}\right)=0,\\ \rho \left[\frac{\partial \mathbf{u}}{\partial t}+(\mathbf{u}·\nabla )\mathbf{u}\right]=-\nabla p+\nabla ·\left[\rho \nu \left(\nabla \mathbf{u}+\mathbf{u}\nabla \right)\right]+\rho \mathbf{g},\end{array}$

(19)

where the total density and velocity of the fluid mixture are given, respectively, by

$\begin{array}{c}\rho ={{\displaystyle \sum }}_k{\rho }_k,\\ \rho \mathbf{u}={{\displaystyle \sum }}_k{\rho }_k{\mathbf{u}}_k+\frac{1}{2}{{\displaystyle \sum }}_k{\mathbf{F}}_k\end{array}$

(20)

with a nonideal gas equation of state given by [28].

3.3. Lattice Boltzmann Model for Transport of Reactive Solutes

In a previous article, Kang et al. [29] have derived the following LB equation for the total primary species concentrations for chemical systems with reactions written in canonical form:

$\begin{array}{cc}\hfill G_{\alpha j}^{}\left(\mathbf{x}+{\mathbf{e}}_{\alpha }{\delta }_t,t+{\delta }_t\right)=G_{\alpha j}^{}\left(\mathbf{x},t\right)-\frac{G_{\alpha j}^{}\left(\mathbf{x},t\right)-G_{\alpha j}^{\text{eq}}\left({\mathrm{\Psi }}_j,\mathbf{u}\right)}{{\tau }_{\text{aq}}},& \hfill \\ \hfill \left(j=1,\dots ,N_C\right),& \hfill \end{array}$

(21)

where N _C is the number of primary species, Ψ _j is the total concentration of the jth primary species, G _α,j is its distribution function along the α direction, $G_{\alpha j}^{\text{eq}}$ is the corresponding equilibrium distribution function, e _α are velocity vectors, δ _t is the time increment, and τ_aq is the dimensionless relaxation time.

It has been shown that the above equation can recover the following pore-scale advection-diffusion equation for Ψ _j [30]:

$\frac{\partial {\mathrm{\Psi }}_j}{\partial t}+\left(\mathbf{u}·\nabla \right){\mathrm{\Psi }}_j=\nabla ·\left(D\nabla {\mathrm{\Psi }}_j\right).$

(22)

This equation is the same as (5) except that here the electrochemical migration is neglected. Assuming that the homogeneous reactions are in instantaneous equilibrium, we have the following mass action equation [31, 32]:

$C_i={\left({\gamma }_i\right)}^{-1}{K}_i{\displaystyle \prod }_{j=1}^{N_c}{\left({\gamma }_j{C}_j\right)}^{{\nu }_{ji}},$

(23)

where ν _ji are the stoichiometric coefficients, K _i is the equilibrium constant of the ith homogeneous reaction, γ _i is the activity coefficient of the ith secondary species, and C _j and C _i are solute concentrations for primary and secondary species, respectively. They are related by

${\mathrm{\Psi }}_j=C_j+{\displaystyle \sum }_{i=1}^{N_R}{\nu }_{ji}{C}_i,$

(24)

where N _R is the number of independent homogeneous reactions, or, equivalently, secondary species.

More details on the heterogeneous reactions between aqueous species and minerals at the pore-mineral interface described by (9), and on the update of solid phase, can be found in [29, 33].

3.4. Lattice Poisson Method

To solve the Poisson equation with strong nonlinearity, (8), we adopt the lattice Poisson method (LPM) developed previously [34, 35], which tracks the electrical potential distribution transporting on the discrete lattices. By expanding (8) into the time-dependent form

$\frac{\partial \psi }{\partial t}={\nabla }^2\psi +g_{\text{rhs}}(\mathbf{r},\psi ,t)$

(25)

with $g_{\text{rhs}}=(1/\epsilon {\epsilon }_0){\sum }_i^{}{z}_{i}e{n}_{i,\infty }\exp (-(z_{i}e/k_{b}T)\psi )$ representing the negative right-hand side (RHS) term of the original (8), we get the discrete evolution equation for the electrical potential distribution

$\begin{array}{cc}\hfill & g_{\alpha }(\mathbf{r}+\mathrm{\Delta }\mathbf{r},t+{\delta }_{t,g})-g_{\alpha }(\mathbf{r},t)\hfill \\ \hfill & \hspace{1em}\hspace{1em}=-\frac{1}{{\tau }_g}\left[g_{\alpha }(\mathbf{r},t)-g_{\alpha }^{\text{eq}}(\mathbf{r},t)\right]+(1-\frac{0.5}{{\tau }_g}){\delta }_{t,g}{\omega }_{\alpha }{g}_{\text{rhs},}\hfill \end{array}$

(26)

where $g_{\alpha }^{\text{eq}}$ is the equilibrium distribution of the electric potential evolution variable. The time step for the electrical potential evolution is

${\delta }_{t,g}=\frac{{\delta }_x}{c^{\prime }},$

(27)

where c′ is a pseudo-sound speed in the potential field. After evolving on the discrete lattices, the mesoscopic electrical potential can be calculated using

$\psi ={\displaystyle \sum }_{\alpha }(g_{\alpha }+0.5{\delta }_{t,g}{g}_{\text{rhs}}{\omega }_{\alpha }).$

(28)

Although the electrical potential evolution equations are in an unsteady form, only the steady state result is realistic, because the electromagnetic susceptibility has not been considered. Although the lattice evolution method for nonlinear Poisson equation is not as efficient as the multi-grid solutions due to its long wavelength limit, it has the advantages of suitability for geometrical complexity and parallel computing.

4. Simulation Examples

4.1. Injection of CO₂ into a Limestone Rock

We first present some modeling results on the injection of a fluid saturated with 170 bars CO₂(g) into a limestone rock at the pore scale. The pore structure was derived from a digitized image of a limestone rock thin section with 640 × 480 pixels (Figure 1). We reduced the original resolution to save computational time (Figure 2). The chemical system of Na⁺-Ca²⁺-Mg²⁺-H⁺-SO₄²-Cl⁻-CO₂ with the reaction of calcite to form dolomite and gypsum is considered. Secondary species included in the simulation are OH⁻, ${\text{HS}{\text{O}}_4}^{-}$ , H₂SO₄(aq), ${\text{C}{\text{O}}_3}^{2-}$ , ${\text{HC}{\text{O}}_3}^{-}$ , CaCO₃(aq), ${{\text{CaHCO}}_3}^{+}$ , CaOH⁺, CaSO₄(aq), MgCO₃(aq), ${{\text{MgHCO}}_3}^{+}$ , MgSO₄(aq), MgOH⁺, NaCl(aq), NaHCO₃(aq), and NaOH(aq). Initial fluid composition is pH 7.75 and 2.69 m NaCl brine, equilibrium with minerals calcite, dolomite and gypsum at 25°C. Initial rock composition is calcite. Secondary minerals include dolomite and gypsum. For boundary conditions, a constant pressure gradient is imposed across the domain for flow. When flow reaches steady state, a fluid with a pH of 3.87 and in equilibrium with 179 bars CO₂(g) and minerals dolomite and gypsum is introduced at the inlet. Zero gradient boundary conditions are imposed at the outlet. Two different cases are considered with different mineral reaction rates to show their effects on solution concentration, mineral deposition, and change in geometry.

Figure 1

Photographic image of a limestone rock thin section (640 × 480 pixels).

Figure 2

Digitized image with reduced resolution.

Resulting geometries at time = 15625 seconds for two different mineral reaction rate constants are plotted in Figure 3. Damkohler is 7.375 for calcite and gypsum, 0.7375 for dolomite for the faster mineral reactions, $7.375\times {10}^{-2}$ for calcite and gypsum, and $7.375\times {10}^{-3}$ for dolomite for slower reactions. Concentration distribution of total Ca²⁺, Mg²⁺, and ${\text{S}{\text{O}}_4}^{2-}$ , pH, and reaction rates for calcite, dolomite, and gypsum for the slower reactions are plotted in Figure 4. As can been seen from the figures, as the reaction rate constants decrease, the deposition of dolomite becomes more uniform surrounding the dissolving calcite grains. Only a small amount of gypsum forms on top of dolomite. At some point in the simulation, the major pores for flow become blocked halting further fluid flow through the medium. The pH is uniform over the entire pore space. All reaction rates have finite values at the mineral surface in the whole domain, outlining the solid geometry. The reaction rate is negative for calcite and positive for dolomite and gypsum, confirming that calcite is dissolving while dolomite and gypsum are precipitating.

Resulting geometries at time = 15625 seconds for two different mineral reaction rate constants: (a) large reaction rate constants; (b) small reaction rate constants.

Distribution of solute concentrations, pH value, and reaction rates at time = 15625 seconds for small reaction rate constants.

4.2. Two-Phase Behavior and Flooding Phenomena in Polymer Electrolyte Fuel Cells

In this Section, we present some results for two-phase flow through the porous CL and the fibrous GDL in a PEFC. Details can be found in [36]. Figure 5 displays the steady state invading liquid water fronts corresponding to increasing capillary pressures from the primary drainage simulation in the reconstructed CL microstructure characterized by slightly hydrophobic wetting characteristics with a static contact angle of 100°. At lower capillary pressures, the liquid water saturation front exhibits finger like pattern, similar to the displacement pattern observed typically in the capillary fingering regime. The displacing liquid water phase penetrates into the body of the resident wetting phase (i.e., air) in the shape of fingers owing to the surface tension driven capillary force. However, at high saturation levels, the invading nonwetting phase tends to exhibit a somewhat flat advancing front. This observation, as highlighted in Figure 5(b), indicates that with increasing capillary pressure, even at very low capillary number (Ca), several penetrating saturation fronts tend to merge and form a stable front. The invasion pattern transitions from the capillary fingering regime to the stable displacement regime and potentially lies in the transition zone in between. In an operating fuel cell, the resulting liquid water displacement pattern pertaining to the underlying pore-morphology and wetting characteristics would play a vital role in the transport of the liquid water and hence the overall flooding behavior.

Advancing liquid water front with increasing capillary pressure through the initially air-saturated reconstructed CL microstructure from the primary drainage simulation.

Figure 6 shows the liquid water distribution as well as the invasion pattern from the primary drainage simulation with increasing capillary pressure in the initially air-saturated reconstructed carbon paper GDL characterized by hydrophobic wetting characteristics with a static contact angle of 140°. The reconstructed GDL structure used in the two-phase simulation consists of 100 × 100 × 100 lattice points in order to manage the computational overhead to a reasonable level. Physically, this scenario corresponds to the transport of liquid water generated from the electrochemical reaction in a hydrophobic CL into the otherwise air-occupied GDL in an operating fuel cell. At the initially very low capillary pressure, the invading front overcomes the barrier pressure only at some preferential locations depending upon the pore size along with the emergence of droplets owing to strong hydrophobicity. As the capillary pressure increases, several liquid water fronts start to penetrate into the air occupied domain. Further increase in capillary pressure exhibits growth of droplets at two invasion fronts, followed by the coalescence of the drops and collapsing into a single front. This newly formed front then invades into the less tortuous in-plane direction. Additionally, emergence of tiny droplets and subsequent growth can be observed in the constricted pores in the vicinity of the inlet region primarily due to strong wall adhesion forces from interactions with highly hydrophobic fibers with the increasing capillary pressure. One of the several invading fronts finally reaches the air reservoir, physically the GDL/channel interface, at a preferential location corresponding to the capillary pressure and is also referred to as the bubble point. It is important to note that the mesoscopic LB simulations provide fundamental insight into the pore-scale liquid water transport through different GDL structures and would likely enable novel GDL design with better water removal and flooding mitigation.

Figure 6

Advancing liquid water front with increasing capillary pressure through the initially air-saturated reconstructed GDL microstructure from the primary drainage simulation.

4.3. Electroosmosis in Homogeneously Charged Micro- and Nanoscale Porous Media

In this section, we briefly present some simulation results on electroosmotic flows (EOFs) in charged micro porous media using the lattice Poisson-Boltzmann method (LPBM), with geometry effects, solution and surface charge effects considered. More details can be found in [35]. We focus on a cubic system of which each side is 1 micron long. A 60 × 60 × 60 uniform grid is used. We change microstructure geometries of porous media by varying the porosity ε from 0.1 to 0.9. The average characteristic length of particles varies from 20 to 150 nm. Figure 7 shows Schematics of the generated porous structures for porosity 0.3 and 0.6. The bulk ionic concentration n _∞ varies from 10⁻⁶ to 10⁻³ M and the surface zeta potential from 0 to −100 mV. The other properties and parameters used in this work are the fluid density ρ = 999.9 kg/m³, the dielectric constant ${\epsilon }_r{\epsilon }_0=6.95\times {10}^{-10}$ C²/Jm, the dynamic viscosity μ = 0.889 mPas, the temperature T = 273 K and the external electrical field strength E = 1 × 10⁴ V/m.

Schematics of the generated porous structures on 60 × 60 × 60 grid systems. The white is solid particles and the dark is fluid: (a) porosity = 0.6, (b) porosity = 0.3.

First, the geometry effects on the electroosmotic permeability in micro porous media are investigated by changing volume fraction and particle size (or number density) of the solid phase. We define the electroosmotic permeability, κ _e , as

${\kappa }_e=\frac{\overline{u}}{E},$

(29)

where $\overline{u}$ is the averaged velocity of EOF along the direction of the driving electrical field E. The coefficients of electroosmotic permeability (κ _e ) for different porosities (ε) of porous media are shown in Figure 8. The bulk molar concentration c_∞ = 10⁻⁴ M, and ζ = −50 mV. The electroosmotic permeability increases with the porosity monotonically. The increasing rate rises with the porosity as well which is very low when the porosity is smaller than 0.5 and becomes sharply high when the porosity is larger than 0.7. The predicted electroosmotic permeability is in the order of 10⁻⁹ m²/sV, which is consistent with the existing experimental measurements.

Figure 8

Predicted electroosmotic permeabilities for various porosities of porous media at c_∞ = 10⁻⁴ M, ζ = −50 mV, E = 1 × 10⁴ V/m.

Figure 9 shows the predicted electroosmotic permeability versus the bulk ionic concentration of the electrolyte solution. We used ε = 0.38. The electroosmotic permeability κ _e increase monotonically with the bulk ionic concentration c _∞ as c _∞ varies from 10⁻⁶ to 10⁻³ M. This result can be explained by the undeveloped electrical potential distributions in narrow channels, whose similar results can be found in [37, Figure 2] and [38, Figures 1 and 2]. When c _∞ is lower than 10⁻⁴ M, the electroosmotic permeability is nearly proportional to the bulk ionic concentration. When c _∞ is higher, the increasing rate becomes a little smaller.

Figure 9

The electroosmotic permeability changing with the bulk ionic concentration for ε = 0.38, ζ = −50 mV, and E = 1 × 10⁴ V/m.

Zeta potential on solid surfaces of porous media affects EOF permeability directly. Simple proportional relationships have been obtained between the electroosmotic permeability and the zeta potential for electrical transports in soils and in polymer composites recently based on the boundary-layer theory. Here we analyze such effects using our numerical methods.

Figure 10 shows the calculated electroosmotic permeability versus the zeta potential on solid surfaces of porous media. All surfaces are homogeneously charged with a same value of ζ. The other parameters used are: c_∞ = 10⁻⁴ M, c _d = 0.1 and ε = 0.38. The zeta potential ζ changes from 0 to 90 mV. It shows that the proportionally linear relationship between electroosmotic permeability and zeta potential is accurate only when ζ is very small (<30 mV). The permeability increases much sharper when the zeta potential ζ is larger than 40 mV.

Figure 10

The electroosmotic permeability versus the zeta potential for ε = 0.38, c_∞ = 10⁻⁴ M, and E = 1 × 10⁴ V/m.

5. Conclusions

We have presented our recent work on mesoscopic modeling of multiphysicochemical processes in porous media, based on the LBM. For the problem of injecting CO₂ saturated brine into a limestone rock, it is shown that the LBM is able to provide detailed information on fluid velocity, solute concentration, mineral composition, and reaction rates, as well as the evolution of the porous media geometry, and therefore can shed some light on the fundamental physics occurring at the pore scale for reactive transport involved in geologic CO₂ sequestration. For two-phase behavior and flooding phenomena in PEFCs, the LBM is a powerful tool to study the influence of the pore structure and surface wettability on liquid water transport and interfacial dynamics in the PEFC catalyst layer and gas diffusion layer. Particularly, the two-phase regime transition phenomenon in the capillary dominated transport in the CL and the influence of the mixed wetting characteristics on the flooding dynamics in the GDL are demonstrated. For electroosmotic flows in charged porous media, the strongly nonlinear governing equations of electroosmosis in three-dimensional porous media are solved using the LPBM. The effects of pore structure, bulk ionic concentration, and the surface charge on electroosmotic permeability are carefully investigated. It is concluded that the LBM is a powerful numerical tool to simulate multiphysicochemical processes in porous media at the pore-scale.