Solution of the Brinkman Equation: Multiple Scale Porous Media

Modeling fluid flow in porous heterogeneous materials with more than one typical pore size (e.g. concrete, microporous rocks and fractured materials) presents a challenge because it is difficult to simultaneously resolve all the microstructural features of the porous medium that are at different length scales. One possible approach is to divide the porous medium into two regions: ($1$) the larger pores and ($2$) homogeneous regions of smaller pores. In the larger pores, the Stokes' equations for incompressible flow hold:

$$\nabla p = \mu \nabla^2 {\bf v}$$ (12)

$$\nabla \cdot {\bf v} = 0$$ (13)

where $p$ is the pressure, ${\bf v}$ is the fluid velocity and $\mu$ is the fluid viscosity. Regions with the smaller pores are treated as a permeable medium and flow is described by Darcy's law. The two boundary conditions to be satisfied at the pore/permeable medium interface are continuity of the fluid velocity and the shear stress [16,17]. Darcy's law alone is not sufficient to satisfy these boundary conditions. The Brinkman equation [17] is a generalization of Darcy's law that facilitates the matching of boundary conditions at an interface between the larger pores and the permeable medium. Brinkman's equation is

$$\langle \nabla p \rangle = -\frac{\mu}{k}{\bf v}+ \mu_e \nabla^2\langle {\bf v} \rangle$$ (14)

where ${\bf v}$ is the fluid velocity, $\mu$ is the fluid viscosity, and $\mu_e$ is an effective viscosity parameter. The so-called effective viscosity should not be thought of as the viscosity of the fluid but only a parameter that allows for matching of the shear stress boundary condition across the free-fluid/porous medium interface. That is, ( $\mu d\langle v \rangle/dy$ $(y=0^+)$ $=\mu_e d\langle v \rangle/dy$ $(y=0^-)$ where $y=0$ specifies the location of the interface for this example. The $+$ and $-$ refer to regions in the free-fluid and porous medium, respectively.

Although the Brinkman equation is semi-empirical in nature, it has been validated by detailed numerical solution of the Stokes' equations in regions near the interface between dissimilar regions [18]. Numerical solution of the Brinkman equation by more traditional computational methods (e.g. finite difference and finite element) is certainly possible. However, a recent lattice Boltzmann (LB) based model by Spaid and Phelan [19], along with recent improvements [20], has proven to be a simple and computationally efficient method to numerically approximate fluid flow described by the Brinkman equation.

To produce flow consistent with the Brinkman equation, a dissipative forcing ${\bf F}=-\frac{\mu {\bf v}}{k}$ is used. Originally, this forcing was incorporated into a LB model, normally used to approximate the Navier Stokes equations, by introducing a velocity shift, $\Delta {\bf v}=\tau {\bf F}/n$, ($\tau$ is a relaxation parameter and $n$ is the density) in the Boltzmann equilibrium distribution according to the method of Shan and Chen [4]. However it is well known that this approach will produce errors of order $\tau^2 F^2$ in the pressure tensor [20]. Such errors can have a significant impact on the fluid dynamics of such systems. Hence, it can be advantageous to instead apply the force in the body force term.

Figure 7: Velocity field of a sheared system next to a porous medium. The filled triangles and circles represent data from the lattice Boltzmann simulation using ($\mu _e/\mu =1$ and $\mu _e/\mu =4$, respectively). The solid lines are analytic solutions of the Brinkman equation. The region below the dashed line $y=34.5$ (in units of lattice spacing) corresponds to a porous medium. The moving wall is at $y=44$.

Figure 8: Bulk permeability, $k_b$ vs permeability assigned to the normally solid portion of the microstructure. The curves (top to bottom) correspond to microstructres whose initial porosity was 40 % (squares), 22.5 % (Xs), 13.0 % (circles), and 7.5 % (triangles). The isolated data points on the left represent the case where $k_s = 0$.

To first validate this model, a simple Couette flow geometry was used (see Figure 7). Starting with a parallel plate geometry, a permeable medium is positioned such that there is a gap between the permeable medium and the upper plate. The upper plate is given a velocity $V_w$ to the right. Analytic solution of the Brinkman equation predicts a linear velocity profile in the gap and an exponentially decaying velocity profile in the porous medium. The rate of decay depends on the value of $\sqrt{\frac{\mu_e}{\mu}}$ [16]. In Figure 7, velocity profiles are compared for the case of $\frac{\mu_e}{\mu}$=4 and the assumption of $\frac{\mu_e}{\mu}=1$. The solid line is the analytic solution of the Brinkman solution. Clearly, there is excellent agreement between simulation and theory and there can be a considerable change in the velocity profile when $\frac{\mu_e}{\mu}\ne1$. In addition, the lattice Boltzmann method also does a reasonably good job capturing the discontinuity of the gradient of the velocity field at the free-fluid/porous medium interface for the case of $\frac{\mu_e}{\mu}$=4. Note that this is achieved without direct incorporation of the stress boundary condition in the simulation model.

For this test case, $k=1/11$ in units of lattice spacing squared. Such a choice of $k$, ignoring tortuosity effects, corresponds to a porous medium with a typical pore size of order a lattice spacing as can be seen by noting that the permeability associated with a cylindrical tube is $k=r^2/8$ where $r$ is the tube radius.

We next consider solution of the Brinkman equation using the Fontainebleau sandstone, described earlier, as a porous medium where the solid phase is now assigned a permeability $k_s$. Although this may not be the case for the original rock, the sandstone image serves as a convenient "random" pore structure to use. Four different permeable media were used with porosity ranging from about $7$ % to $40$ %. Here, the porosity refers to the original pore structure. In Figure 8, we plot the bulk permeability, $k_b$, of the overall system vs $k_s$. Clearly, $k_b$ increases with $k_s$. At lower values of $k_s$, the higher porosity system appears to be less sensitive to $k_s$ as a result of the the larger pores carrying most of the flow. Fluid flow in the lower porosity systems are much more influenced by $k_s$, as a larger fraction of the system is composed of the permeable medium. It should also be pointed out that an alternate version [21] of the lattice Boltzmann method was used to determine results for $k_s < 0.1$ as the above described model is unstable in this regime. The instability is a result of the fixed time step used in the LB method. The alternate version allows for introducing smaller time steps so that the instability is avoided. The reader is referred to [21] for more details.