0
Research Papers: Fundamental Issues and Canonical Flows

# An Efficient Method of Generating and Characterizing Filter Substrates for Lattice Boltzmann AnalysisOPEN ACCESS

[+] Author and Article Information
John Ryan Murdock

Department of Mechanical
Engineering–Engineering Mechanics,
Michigan Technological University,
Houghton, MI 49931
e-mail: jrmurdoc@mtu.edu

Aamir Ibrahim

Department of Mechanical
Engineering–Engineering Mechanics,
Michigan Technological University,
Houghton, MI 49931
e-mail: aamiri@mtu.edu

Song-Lin Yang

Professor
Fellow ASME
Department of Mechanical
Engineering–Engineering Mechanics,
Michigan Technological University,
Houghton, MI 49931
e-mail: slyang@mtu.edu

Contributed by the Fluids Engineering Division of ASME for publication in the JOURNAL OF FLUIDS ENGINEERING. Manuscript received April 9, 2017; final manuscript received September 28, 2017; published online November 23, 2017. Assoc. Editor: Praveen Ramaprabhu.

J. Fluids Eng 140(4), 041203 (Nov 23, 2017) (11 pages) Paper No: FE-17-1216; doi: 10.1115/1.4038167 History: Received April 09, 2017; Revised September 28, 2017

## Abstract

To provide porous media substrates that are quick to generate and characterize for lattice Boltzmann analysis, we propose a straightforward algorithm. The method leverages the benefits of the lattice Boltzmann method (LBM), and is extensible to multiphysics flows. Several parameters allow for simple customization. The generation algorithm and LBM are reviewed, and suggested implementation covered. Additionally, results are discussed and interpreted to evaluate the approach. Several verification tools are employed such as Darcy's law, the Ergun equation, the Koponen correlation, a newly proposed correlation, and experimental data. Agreement and repeatability are found to be excellent, suggesting this relatively simple method is a good option for engineering studies.

<>

## Introduction

Disordered media flow analysis can provide significant design direction and research findings across a variety of fields such as oil and gas, chemical, automotive, and environmental engineering and sciences. However, the multiscale nature of porous flow presents challenges not present in many areas of computational fluid dynamics. Three length scales broadly describe most media: domain scale, representative elementary volume, and pore scale. Domain scale provides the least detail and effectively acts as a black box. Behavior may be governed by simple pressure drop versus flow rate assignments. The representative elementary volume scale provides more details governed by semi-empirical correlations such as Darcy's law, the Kozeny–Carman relation, and the Klinkenberg, Forchheimer, and Brinkman considerations [1,2].

It is at the pore scale where more fundamental governing equations can describe and display the flow for deeper understanding. Unlike the previous two descriptions, the pore structures must be characterized, and put in a form for a discretized solution method. There are two main approaches: Acquisition of the solid matrix by an imaging device, such as X-ray computed tomography (CT scan), as used in Ref. [3] and construction of the solid matrix by an algorithm that follows certain properties [48]. The focus here is on a simpler, but still effective algorithm construction to provide greater flexibility, speed, and variety than a scanning process can.

Not all mathematical techniques for resolving flow fields are ideal for pore-scale flows. The choice of technique depends on important considerations like handling of complex boundary conditions and the ease with which the scheme can be parallelized [9]. The lattice Boltzmann method (LBM) utilizes the relatively simple bounce-back scheme on arbitrarily complex solid surfaces and momentum inlets [10,11], is relatively easy to implement [9], and due to its linear nature, scales well when parallelized [12]. The lattice Boltzmann equation (LBE), however, introduces a compressibility error for incompressible flows [13], such as most filter flows [14]. In order to reduce this error, an incompressible method (marked with an “i”) is utilized here [15].

In this study, a simpler algorithm for filter substrate generation for lattice Boltzmann analysis is proposed. A filter substrate refers to the underlying material, which provides the pore structure, commonly cordierite for diesel particulate filtration. The generated substrates are compared against semi-empirical models with an incompressible LBE to verify the algorithm. Flows and conditions are representative of common filter conditions, such as $Re ≲100$. Capability of the generated substrate to handle thermal and multiphase flows and sweeps of various porosities is discussed. Modifications for various types of porous media and three-dimensional flows are also considered for future work.

## Substrate Generation and Measurement Algorithm

Defining the simulation domain is fundamentally a task of marking nodes of a lattice as solid or fluid on an Nx × Ny state matrix. This state matrix must provide the ability to control porosity, and ideally the ability to control a characteristic size. It must also be possible to make the necessary measurements such as porosity, wetted perimeter, and flow area to obtain characteristics like hydraulic radius and Reynolds number (Re).

To comply with the variety of languages for LBM code implementation, a pseudocode with broad applicability is available at the website link.1 The code provides functions not only for the substrate generation but also for measurement, and to flag boundary nodes. By flagging boundary nodes during preprocessing, LBE streaming and bounce back occurs only on the necessary nodes, introducing an efficiency improvement.

Generation is based on the random placement of ellipses that can overlap to form more complex shapes. Each ellipse has a major and minor axis of random radii within a designated range of values. Nodes that fall within the ellipse are flagged as “solid” in the state matrix. Only certain portions of the lattice are allowed to contain solids to provide inlet and outlet areas, both of which can be dictated. To ensure these areas are obeyed, the portions of ellipses that were assigned to the inlet and outlet areas are reset to “fluid.” The size of the domain, number of ellipses, and radii ranges are important choices, and will be discussed further in Secs. 4 and 5.

This compact method produces substrates with porosity within 1% of that assigned, once calibrated. Shape of the solid substrate is more complex than if circles or squares formed the fundamental substrate. Since an ellipse can easily be described by an equation with a third dimension, the generation method can be readily extended. An extension to three-dimensional would provide a percolation threshold below 0.4, as opposed to this initial study on two dimension, which limits the porosity to values greater than 0.4 to avoid the zero permeability limit.

Boundaries are marked with a search algorithm, and the link(s) of the lattice between fluid and solid nodes are recorded. As a result, operations in the LBM are performed only where necessary. Since the perimeter computation is done in the initialization stage, it is not repeated during simulation.

Finally, measurement is conducted based on inputs from the previous functions. The discrete nature of the lattice makes characterization fundamentally a counting task. By providing porosity and wetted perimeter, the nondimensional values that verify and characterize the filter can be computed. These equations and implementations are discussed in the sections on methodology and results.

Figures 1 and 2 demonstrate samples of the generated substrates.

## The Lattice Boltzmann Method

By solving the Boltzmann equation for streaming and colliding particles, fluid flow can be resolved in a more fundamental manner than the Navier–Stokes equations. Discretization results in the lattice Boltzmann equation, and allows for a computational approach to fluid dynamics.

###### The Incompressible Lattice Boltzmann Equation.

The LBE is stated as [16] Display Formula

(1)$fix+ciδt,t+δt=fix,t-1τfix,t-fieqx,t$

where δt is the time-step, f is the velocity distribution function, c = δx/δt is the lattice speed, δx is the spatial discretization, i is the index of the lattice velocity, τ is the collision relaxation time related to viscosity, and f(eq) is the equilibrium distribution function based on the Maxwell–Boltzmann distribution.

The form of the equilibrium distribution function largely dictates the characteristics of the flow solution. The standard form is [16] Display Formula

(2)$fi(eq)=wiρ1+(ci.u)cs2+(ci.u)22cs4-(u.u)2cs2$

where w is the velocity weight, ρ is density, and cs is the speed of sound. Macroscopic density and momentum are computed by Display Formula

(3)$ρ=∑ifi,ρu=∑icifi$

At small Mach numbers and density variations, the error in simulating incompressible flows is limited. For porous flows, this can pose a problem due to large pressure and density gradients [14].

Past filtration studies have utilized an LBE, which limits the compressibility error [17] with a pseudo-compressibility approach. However, improvement in achieving an incompressible behavior for steady and transient flows can still be made [15]. The incompressible singular equilibrium form is Display Formula

(4)$fieq=-5P3c2+Siu,i=0P3c2+Siu,i=1,2,3,4P12c2+Siu,i=5,6,7,8$
with Display Formula
(5)$Siu=wi3(ci.u)c+92(ci.u)2c2-32(u.u)c2$

where P is pressure. Density is no longer a variable, and a speed of sound is no longer present. The macroscopic variables are computed as Display Formula

(6)$u=∑icifi,P=−c25[2c2(u.u)+3f0(eq)]$

Through the Chapman–Enskog expansion, the incompressible Navier–Stokes equations are recovered as Display Formula

(7)$∇.u=0∂u∂t+∇.uu=-∇P+ν∇2u$

where viscosity and the collision relaxation time are related by Display Formula

(8)$ν=δx23δtτ-12$

When the LBM has been applied to pore-scale simulations, convergence and accuracy were found to be compromised due to artifacts in the boundary conditions. Pan et al. recognized this issue and found that a multiple relaxation time (MRT) form of the LBE avoided the problem [18]. More specifically, viscosity in the single relaxation time (SRT) form has an effect on where the boundaries are located [10]. While the permeability should scale linearly with viscosity (all other quantities being equal), Pan et al. find unwanted deviations using SRT, in particular at low viscosities. Implementing the SRT form in our own models results in the same trend.

The MRT–LBE replaces the single value τ with a collision matrix, which allows each mode to relax at a different rate. Due to these benefits, the flow solver for this study is the MRT-iLBE. d'Humieres introduces the MRT–LBE, and Murdock et al. provide derivations, forms, and results of using the MRT-iLBE, which is second-order accurate in time and space [1921].

###### The Lattice.

As Secs. 1 and 2 have discussed, the substrate matrix is formed on, and the solution takes place on a lattice of nodes. For this study, that lattice is D2Q9, where D labels the number of spatial dimensions, and Q labels the number of discrete velocities. Figure 3 shows the lattice structure used here.

Values discussed in the section on the LBE, such as w, c, and the subscript i, all made reference to properties of this lattice. For D2Q9, wi = {4/9, 1/9, 1/9, 1/9, 1/9, 1/36, 1/36, 1/36, 1/36} for i ={0–8}. Since this is a square lattice, c = 1 and ci = (0,0),(1,0), (0,1),(−1,0),(0,−1),(1,1),(−1,1),(−1,−1),(1,−1) for i = {0–8}. In the standard LBE case, cs2 = 1/3.

###### Boundary Conditions.

For porous media flows, the major benefit of the LBM is boundary condition implementation. The bounce-back method is commonly implemented in LBM for complex geometries. The algorithm is simple. As shown in Fig. 3, if a value for f8 streams into a wall, it is simply returned along the lattice vector 6. Even without the perimeter and link search of the algorithm presented in Sec. 2, the boundary conditions are easy to implement and efficient just by knowing which lattice nodes are solid and which are fluid.

Boundaries conform to the Cartesian grid, and form a stair-step approximation to the pore structures. While this is an approximation to potentially more complex surfaces, it allows for implementation ease and computational efficiency, and is frequently applied to LBM filter simulations with success [14,17]. Implementing this approach for common filter flow conditions, we find good agreement with verification techniques in the later sections of this work. Accuracy of these complex boundaries is improved with increasing global grid refinement. This can be expensive; however, the need to refine the entire domain to capture the random distribution of small pore passages partially negates the benefits of a local refinement approach. This simple method is applied on all substrate surfaces in the domain.

The same boundary method can be modified for symmetry and inlet boundary conditions. To provide the free-slip behavior of symmetry bounds, the value of f8 is specularly reflected along the lattice vector 5. For the top and bottom boundaries of the substrate domain, this method is implemented to simulate the length of an entire filter wall. At the inlet, the value of fi bounced back into the domain is modified with a momentum source term to generate the desired velocity [11]. The boundary (wall, symmetry, or inlet) is actually placed at the point halfway between solid and fluid node, and this is accounted for in our algorithm. With this half-way interpolation, the boundaries are second-order accurate [10,21].

Outlets are the only boundary, which cannot utilize a bounce-back scheme. Instead, an extrapolation scheme [22], similar to the finite difference method, is implemented, which is second-order accurate.

###### Multiphysics Capability.

Since transport phenomena are described by equations fundamentally similar to the Navier–Stokes equations, the lattice Boltzmann equation described in Secs. 3.13.3 can be readily extended to additional physics [23]. Further, since the substrate is generated to be compliant with the LBM, these additional transport LBEs can be employed. Nonisothermal incompressible flow is straightforward to implement, and since velocity can simply be turned to zero at solids, conjugate heat transfer on the substrate is not difficult to implement [24]. While viscosity is replaced by thermal diffusivity for energy transport, the same equation can be employed with mass diffusivity to transport the species to be filtered [25].

## Verification Methodology

The following empirical and derived equations and correlations provide the basis for utilizing the algorithm presented in Sec. 2 for filter studies.

###### Darcy's Law.

Darcy's law is a widely employed model relating the porosity (φ) and the ease with which flow passes through disordered media, measured by permeability (κ). For a fluid of viscosity ν, density ρ, and velocity u, the pressure drop should result in a permeability given by Ref. [2]. Pressure drop in Eq. (9) is calculated using the average value of pressure at the inlet and outlet, respectively Display Formula

(9)$κ=uνρ∇P$

Since the pore-scale flow is being resolved, permeability is a direct result of the simulation algorithm, as opposed to an input or tuned result. Further, with the use of the MRT form, permeability behaved linearly with viscosity changes in our models. The relationship between porosity and permeability is linear on a semilog plot. Simple calculation and expectations make this a good initial methodology for verification of pore-scale simulations.

###### Ergun Equation.

The Ergun equation relates pressure losses in the form of the friction factor (f) to the Reynolds number of a packed column and fluidized beds. The equation is a correlation without the foundation of the Navier–Stokes equations like Darcy's law. However, it is a useful tool in determining the validity of methods for analyzing fluid flows in porous media. The equation is stated here as [26] Display Formula

(10)$f=150Re+1.75$

where Re is defined as Display Formula

(11)$Re=uDpν(1-φ)$

For three dimensions, the hydraulic diameter is defined as Display Formula

(12)$Dp=6Rh(1-φ)φ$

where the hydraulic radius Rh is the ratio of flow volume to wetted surface area.

Since two dimensions are of interest in this initial study, the definition of Dp must be modified. Since the radius is determined by the ratio of volume to surface area in three dimensions, it is possible to “remove” one dimension from each measure in two dimensions, yielding the ratio of flow area to wetted perimeter. This ratio reduces to Dsphere/6 in three dimensions and Dcircle/4 in two dimensions. As a result, the two-dimensional hydraulic diameter is defined as Display Formula

(13)$Dp=4Rh(1-φ)φ$

Friction factor from the simulation is computed as Display Formula

(14)$f=∇PDpφ3u(1-φ)$

and compared with Eq. (14). On a log–log plot, the Ergun equation forms a more complex nonlinear shape, providing a more challenging verification test.

###### Koponen Model.

A simple exponential model is fitted to the porosity–permeability data as proposed by Koponen for creeping flow through large three-dimensional random fiber webs [27]. However, there is no theoretical basis for this model in literature as also stated in Ref. [27]. This model is valid for values of porosity between 0.42 and 0.85. The model has an essential shortcoming. The zero permeability limit is not reached at a finite porosity value, called the percolation threshold, φc. Despite these shortcomings, the model provides a good verification test for our random substrate generator, the model being valid over a wide range of porosities and also having shown good agreement with experimental data [28]. The relation for the model expressed in nondimensional form is Display Formula

(15)$κDp2=αexpβ1-φ-1$

The following values are found for constants α and β in Eq. (15):

$α=2.44,β=13.94$

###### A New Semi-Empirical Model.

A number of correlations have been proposed since the late 1940s trying to relate porosity and permeability in porous media. Broadly, these correlations can be classified into three types:

1. (1)Relations based on the flow through conduits model.
2. (2)Relations based on the flow past a submerged body model.
3. (3)Relations based on the cell model theory.

At low and intermediate porosities, the porous flow is modeled better by the flow through conduits approach while at higher values of porosity, the flow past a submerged body model or the cell model theory gives better results [29]. The model usually used for flow through filters, called the Kubawara model [30], is based on the cell model theory. Hence, it slightly overestimates the permeability values at low porosities, as it does not take into account the effect due to surrounding fibers in the porous structure. A model based on the conduit-flow approach would give a better approximation at low and intermediate porosity values. Though at higher values of porosity, the Kubawara model gives very good results [31]. Apart from this, the Kubawara model does not take into account the percolation threshold (the minimum porosity of the medium above which it becomes permeable). The idea to improve the permeability model is based on using a weighted combination of two models, a conduit-based model for low porosity values and a submerged body model or a cell model for higher porosities.

###### The Conduit-Model: Gebart Model.

Gebart's formulation is based on modeling the porous structure as a parallel array of cylinders perpendicular to the flow direction and modeling the flow through the voids as flow through conduits of very small cross section or diameter [32]. If Gebart's derivation is followed for cylinders with elliptical cross section (the shape of the solid substrate used in the substrate generator being an ellipse) instead of circles in hexagonal arrangement, the following result is reached: Display Formula

(16)$κ=49π61-φc1-φ-152a32b3(a+b)52$

where a and b are the semi-major and the semiminor axis of the ellipses, respectively. The second term on the right-hand side in the above equation in brackets () is simply a length scale squared. This can be replaced by a characteristic length scale, Dp2 (chosen to be the hydraulic diameter as defined in Eq. (13)), and the above equation can be written in nondimensional form as follows: Display Formula

(17)$κDp2=49π61-φc1-φ-152$

###### Flow Past a Submerged Body Model and the Cell Model.

The model chosen for modeling permeability at higher porosities is the Brinkmann model [33]. This model is based on flow past a swarm of spherical particles. The resulting equation for permeability based on this model in nondimensional form is Display Formula

(18)$κDp2=1723+41-φ-381-φ-3$

The Kubawara model, based on the cell model theory [31], can also be used in place of the Brinkmann model for intermediate and higher porosities and gives exactly the same results. The Kubawara model also improves the range of validity, as the Brinkmann model is not valid for porosities lower than 0.33 [29]. The relation for the Kubawara model in nondimensional form is given by Display Formula

(19)$κDp2=292-951-φ13-φ-(1-φ)21-φ$

###### Tortuosity Factor.

The models mentioned above in Eqs. (17)(19), like most permeability models that suffer from the drawback of not including the effects of the twisted flow path through the porous structure, are quantified by the tortuosity of the porous structure, T. In the present model, an empirical relation for tortuosity is used to incorporate its effect. The following correlation is used to model tortuosity as used by Carman [29]: Display Formula

(20)$T=LeL2$

where $L$ is the length of the porous structure in the macroscopic flow direction, and the value of $Le$ depends on the type of porous media under consideration. The value of $L$ for the present simulation is 200 units.

###### Blending Functions/Basis Functions.

The two models are finally blended together using the idea of blending functions or basis functions. The particular choice of the basis set allows the Gebart form to be the major contributor to the relation at low porosities and the Brinkmann model or the Kubawara model to take over at high porosities >0.75. The basis set chosen is as follows: Display Formula

(21)$w1φ=a1+a2ea3(1-φ)$
Display Formula
(22)$w2φ=b1+b2eb3φ$

The blending function constants a1, a2, a3, b1, b2, and b3 are found using the values of the blending functions, W1(φ) and W2(φ) at the three nodal values that include the end nodes φ = φc (the percolation threshold) and φ = 0.9 (the last value of φ in our numerical simulation results) and a midnode chosen to be φ = 0.75. The values of the blending functions at the nodes are: W1(φc) = 1, W1(0.9) = 0, W2(φc) = 0, W2(0.9) = 1 and W1(0.75) = 0.4 = W2(0.75) = 0.4. The blending functions, W1(φ) and W2(φ), with $φc=0.135$ are plotted against porosity in Fig. 4.

###### The Model Equations.

The new model is based on the idea of modeling the permeability at low porosities using the conduit flow approach and the permeability at high porosities using the flow past a submerged body or the cell-model approach. The combined model is based on a more flexible form of the Gebart relation given by Nabovati et al. [34], and the Brinkmann relation given by Eq. (18) or the Kubawara relation given by Eq. (19). The following nondimensional model equations are arrived at based on above discussions for the Brinkmann and Kubawara models, respectively. Both equations are divided by the tortuosity factor, T as mentioned in Eq. (20). The relation for the combined Gebart–Brinkmann (GB) model is Display Formula

(23)$κDp2=1α/L2w1φ49π61-φc1-φ-1γ+1β/L2w2φ1723+41-φ-381-φ-3$

and for the combined Gebart–Kubawara (GK) model is Display Formula

(24)$κDp2=1α/L2w1φ49π61-φc1-φ-1γ+1β/L2w2φ292-951-φ13-φ-(1-φ)21-φ$

where W1(φ) and W2(φ) are given by Eqs. (21) and (22), respectively. Gebart derives the value of the percolation threshold in the above equations for a hexagonal arrangement of circles [32]. A more exact value of φc = 0.135 is chosen from a numerical simulation of random arrangement of ellipses as shown in Ref. [35] and used in Eq. (24). The values of the blending function constants for the GK model with φc = 0.135 and specific nodal conditions are a1 = 1.12, a2 = −1.49, a3 = −2.91, b1 = −0.01, b2 = 0.005, and b3 = 5.93, respectively. As for the Gebart–Brinkmann model in Eq. (23), the value of φc = 0.33 is used, this being the lower limit in the range of validity of the Brinkmann model. The values of the blending function constants for the GB model with φc = 0.33 are found to be a1 = 1.38, a2 = −1.72, a3 = −2.25, b1 = −0.04, b2 = 0.0064, and b3 = 5.65, respectively. A curve of the form Eqs. (23) or (24) is then fit to the data generated from the iLBM simulation using the nonlinear least-squares method. The following values are found for constants α, β, and γ in Eqs. (23) and (24), respectively:

$α=140.90,β=106.18,γ=2.17$
$α=133.71,β=171.51,γ=3.84$

The values found for γ are around the same range as found by Nabovati et al. [34] and Gebart [32].

## Results

For each study point, a new filter is generated to demonstrate repeatability. Unless otherwise marked, units are nondimensional lattice Boltzmann units. Inlet velocities are varied between 0.001 and 0.1 depending on porosity, for low Mach numbers, and outlet pressures are 1.0. All values are initialized to 0 and the solution proceeds until convergence is achieved at an root-mean-square velocity change of $10-9$. In all cases, τ = 1.0, so Re is computed and reported.

A number of grid sensitivity tests indicate that porosities greater than about 0.8 can utilize a grid of only 300 × 401 nodes, while lower porosities require a grid of 300 × 601 to ensure effective flow passages. 50 nodes upstream and downstream of the filter section bring the boundary conditions sufficiently far from the high gradients of the complex filter flow, leaving a 200 node-wide effective filter wall. The domain can be visualized as shown in Fig. 5. While both the number and radii of ellipses can be used to control porosity, the number of ellipses is modified and radii ranges are assigned as 5–10 nodes here. This forms substrates more in line with cordierite, but other disordered media may be better matched by adjusting the allowed radii.

Velocity fields of the sample filter substrates in Figs. 1 and 2 are visualized in Figs. 6 and 7, respectively. Velocities are normalized by the inlet value. The apparent rivulets are the flow paths through the substrate, with darker colors representing higher velocities, and lighter colors signifying low velocity (the zero velocity of the substrate solids is the lightest). Visually, the flow paths are sensible, and the greatest mass flow occurs where resistance is lowest.

###### Verification by Darcy's Law.

In Fig. 8, permeability is plotted as a function of porosity for five filters at different porosities, at three different velocities. The goal is to evaluate if permeability of the generated filter is intrinsic, and if the relationship is linear on a semilog plot.

In all cases, it is shown that permeability is independent of the specific filter and velocity. Permeability is only a function of porosity, within numerical error.

###### Verification by Ergun Equation.

For five filters of porosities 0.5, 0.6, 0.7, 0.8, and 0.9 a variety of inlet velocities were applied to vary Re. Friction factor at each Re was computed by Eq. (14), and plotted on a log–log scale in Fig. 9. The Ergun equation, Eq. (10) is plotted with the data as a solid line.

The results provide a good match to the correlation. Since visual comparison to the empirical correlation is important in this case, it is worth considering the results of other studies. When set beside Fig. 6.4-2 of Bird's Transport Phenomena [26], the filter substrate generated by the algorithm of this study, and analyzed by the MRT-iLBE, is at least as good a match as experimental results from previous works.

###### Verification by Porosity–Permeability Correlations.

Permeability is evaluated at different values of porosity generated using the random substrate generator. The values of permeability are nondimensionalized using the hydraulic diameter, Dp defined in Eq. (13). The nondimensional permeability, κ/Dp2 is then plotted against porosity, and compared against existing Koponen correlation [27] and the new porosity–permeability models developed in Sec. 4.

###### Verification by Koponen Model.

An equation of the form Eq. (15) is fitted to the iLBM simulation data and the values of nondimensional permeability, κ/Dp2 from the numerical simulation produce an excellent fit with the Koponen correlation within its valid range (0.42<φ < 0.85) [27], as seen in Figs. 10 and 11, cementing the effectiveness of the presented random substrate generator.

###### Verification by Gebart–Brinkmann Model.

As explained in Sec. 4.4, the GB model better approximates the underlying physics of flow through porous media capturing the essential concept of Percolation threshold at low porosities and at the same time providing sound results at high porosities in agreement with previous models. The model strongly verifies our substrate generator showing that it is in compliance with the physics of porous flows at low as well as high porosities as can be seen in Figs. 12 and 13.

###### Verification by Gebart–Kubawara Model.

The GK model fit to the data in Figs. 14 and 15 essentially shows the same trend as the results in Sec. 5.3.2 for the GB model in Figs. 12 and 13. This result is expected as explained in Sec. 4.4. The flow at high porosities in porous media can be modeled using the Cell Model theory as in the GK model or the flow past a submerged body theory as in the GB model. Both models are known to show similar results as stated in Ref. [31]. An additional benefit of the GK model is its wider range of validity at low porosities in comparison to the GB model.

###### Experimental Verification.

Finally, the data for nondimensional permeability found using the iLBM simulation and the random substrate generator are compared against data from experiments in Refs. [27] and [28]. The experiments in Refs. [27] and [28] were performed with fiber mats and fibrous filters, and the permeability was nondimensionalized using a characteristic diameter, Dc = 5. Therefore, the experimental results were scaled by a factor of Dc2/Dp,avg2 and then compared with our simulation results. The value of Dp,avg is taken to be the average of all values of Dp in the plotted range of φ. For the present simulation, the value of Dp,avg is 12.43. The experimental results show very good agreement with our simulation data as shown in Fig. 16. This confirms the idea that our substrate generator can simulate different types of porous media effectively.

## Conclusions

Filtration and flows through disordered porous media play a substantial role in many physical and engineering processes. However, methods for obtaining or generating substrate geometry for the well-suited lattice Boltzmann method can be cumbersome or physically prohibitive. This work provides a pseudocode and description of a method fundamentally built on randomly sized and placed ellipses, which is flexible, quick, and easy to implement for engineering studies. The subsequent results section utilized well-established semi-empirical tools to test the methods validity, such as the Ergun equation and the Koponen correlation. The data collected from such runs produced intricate well-visualized images for pore-scale understanding, something not possible with model equations or representative volume analyses. More importantly, the results showed excellent agreement ranging from the fundamental Darcy's law to the most complex correlations and multisubstrate experimental data. Due to the method's ability to produce very complex shapes and patterns from what is a fundamentally simple shape, the ellipse, the performance matches physical scan and more complex generation methods without over complicating a key step in many studies. Thus, it is recommended as a porous media study tool.

Not only is the method shown to repeatedly produce valid substrates, but also the study reveals important computational fluid dynamics methodology in pore-scale analysis. A singular form incompressible MRT lattice Boltzmann equation was important in limiting compressibility effects and small geometry feature artifacts. The bounce-back scheme, and its variants, handles the inherently random geometry in a tractable manner. Reasonable lattice Boltzmann unit boundary condition values for primitive variables such as velocity and pressure are also a result of the successful analyses. This study adds merit to the growing interest in the lattice Boltzmann method in complex geometry flows.

Of additional interest is that the study of this method resulted in a secondary development in the form of a blending function-based correlation describing permeability variations in different types of porous media over a wider porosity range. Each of the contributing model equations describes the phenomena in different porosity regimes, but together increases the breadth of validity. The blending function correlations are valid at porosities as high as 0.9 as can be seen from Figs. 12 and 14. The model also successfully incorporates the phenomenon of zero permeability at a low finite value of porosity (percolation threshold, $φc$), as seen in Figs. 13 and 15. This essentially helps to differentiate between porosity and effective porosity, which is of utmost significance in all real porous media.

Since porous media is commonly a filtration mechanism, multiphysics capabilities are an interest in future studies. Conjugate heat transfer in nonisothermal filter flows should be straightforward based on the recent lattice Boltzmann literature. A multitude of methods are available for multiphase transport, deposition, and erosion, including layered advection–diffusion lattice Boltzmann equations and lattice gas cellular automata, both compatible with the substrates generated by the studied algorithm. Also of future interest is the extension to three dimensions, which is possible due to the fundamental use of basic mathematical equations valid in two and three dimensions to turn simple shapes into complex patterns.

## Acknowledgements

The authors gratefully acknowledge the Michigan Technological University Computing & Visualization resource, Superior, which was the vehicle for performing these simulations.

## Nomenclature

• c =

lattice speed

• cs =

speed of sound

• Dp =

hydraulic diameter

• f =

friction factor

• fi =

velocity distribution function

• i =

incompressible form

• Re =

Reynolds number

• T =

tortuosity

• w =

lattice weight

• δ =

discrete unit

• κ =

permeability

• ν =

viscosity

• ρ =

density

• τ =

collision relaxation time

• φ =

porosity

## References

Bear, J. , 1972, Dynamics of Fluids in Porous Media, Dover Publications, Mineola, NY.
Whitaker, S. , 1986, “ Flow in Porous Media—I: A Theoretical Derivation of Darcy's Law,” Transp. Porous Media, 1(1), pp. 3–25.
Elkatatny, S. , Mahmoud, M. , and Nasr-El-Din, H. A. , 2013, “ Filter Cake Properties of Water-Based Drilling Fluids Under Static and Dynamic Conditions Using Computed Tomography Scan,” ASME J. Energy Resour. Technol., 135(4), p. 042201.
Pilotti, M. , 1998, “ Generation of Realistic Porous Media by Grains Sedimentation,” Transp. Porous Media, 33(3), pp. 257–278.
Maier, R. S. , Kroll, D. M. , Benard, R. S. , Howington, S. E. , Peters, J. F. , and Davis, H. T. , 2000, “ Pore-Scale Simulation of Dispersion,” Phys. Fluids, 12(8), pp. 2065–2079.
Madadi, M. , and Sahimi, M. , 2003, “ Lattice Boltzmann Simulation of Fluid Flow in Fracture Networks With Rough, Self-Affine Surfaces,” Phys. Rev. E, 67, p. 026309.
Zhang, H. , Ge, X. , and Ye, H. , 2006, “ Randomly Mixed Model for Predicting the Effective Thermal Conductivity of Moist Porous Media,” J. Phys. D, 39(1), pp. 220–226.
Wang, M. , Wang, J. , Pan, N. , and Chen, S. , 2007, “ Mesoscopic Predictions of the Effective Thermal Conductivity for Microscale Random Porous Media,” Phys Rev. E, 75(3), p. 036702.
Guo, Z. , and Shu, C. , 2013, Lattice Boltzmann Method and Its Applications in Engineering, World Scientific Publishing, Hackensack, NJ.
He, X. , Zou, Q. , Luo, L. S. , and Dembo, M. , 1997, “ Analytic Solutions of Simple Flow and Analysis of Non-Slip Boundary Conditions for the Lattice Boltzmann BGK Model,” J. Stat. Phys., 87(1), pp. 115–123.
Ladd, A. , 1994, “ Numerical Simulations of Particulate Suspensions Via a Discretized Boltzmann Equation—Part I: Theoretical Foundation,” J. Fluid Mech., 271, pp. 285–310.
Wellein, G. , Lammers, P. , Hager, G. , Donath, S. , and Zeiser, T. , 2006, “ Towards Optimal Performance for Lattice Boltzmann Applications on Terascale Computers,” International Conference Parallel Computational Fluid Dynamics (ParCFD), Indianapolis, IN, May 14–17, pp. 31–40.
He, X. , Doolen, G. , and Clark, T. , 2002, “ Comparison of the Lattice Boltzmann Method and the Artificial Compressibility Method for Navier–Stokes Equations,” J. Comp. Phys., 179(2), pp. 439–451.
Kang, Q. , Lichtner, P. , and Janecky, D. , 2010, “ Lattice Boltzmann Method for Reacting Flows in Porous Media,” Adv. Appl. Math. Mech., 2(5), pp. 545–563.
Murdock, J. , and Yang, S. , 2016, “ Alternate and Explicit Derivation of the Lattice Boltzmann Equation for the Unsteady Incompressible Navier–Stokes Equation,” Int. J. Comput. Eng. Res., 6(12), pp. 47–59.
Hudong, C. , Chen, S. , and Matthaeus, W. , 1992, “ Recovery of the Navier–Stokes Equations Using a Lattice Boltzmann Method,” Phys Rev. A, 45(8), p. R5339. [PubMed]
Yamamoto, K. , and Ohori, S. , 2012, “ Simulations on Flow and Soot Deposition in Diesel Particulate Filters,” Int. J. Engine Res., 14(4), pp. 333–340.
Pan, C. , Luo, L. S. , and Miller, C. , 2006, “ An Evaluation of Lattice Boltzmann Schemes for Porous Medium Flow Simulations,” Comput. Fluids, 35(8–9), pp. 898–909.
d'Humieres, D. , 1992, “ Generalized Lattice-Boltzmann Equations,” Prog. Astro. Aeron., 159, pp. 450–458.
Murdock, J. , Ickes, J. , and Yang, S. , 2017, “ Transition Flow With an Incompressible Lattice Boltzmann Method,” Adv. Appl. Math. Mech., 9(5), pp. 1271–1288.
Kruger, T. , Kusumaatmaja, H. , Kuzmin, A. , Shardt, O. , Silva, G. , and Viggen, E. M. , 2017, The Lattice Boltzmann Method: Principles and Practice, Springer International Publishing, Cham, Switzerland.
Chen, S. , Martinez, D. , and Mei, R. , 1996, “ On Boundary Conditions in Lattice Boltzmann Methods,” Phys. Fluids, 8(9), pp. 2527–2531.
Yuan, P. , and Schaefer, L. , 2005, “ A Thermal Lattice Boltzmann Two-Phase Flow Model and Its Application to Heat Transfer Problems—Part 2: Integration and Validation,” ASME J. Fluids Eng., 128(1), pp. 151–156.
Karani, H. , and Huber, C. , 2015, “ Lattice Boltzmann Formulation for Conjugate Heat Transfer in Heterogeneous Media,” Phys. Rev. E., 91(2), p. 023304.
Inamuro, T. , Yoshino, M. , Inoue, H. , Mizuno, R. , and Ogino, F. , 2002, “ A Lattice Boltzmann Method for a Binary Miscible Fluid Mixture and Its Application to a Heat-Transfer Problem,” J. Comp. Phys., 179(2), pp. 201–215.
Bird, R. , Stewart, W. , and Lightfoot, E. , 2006, Transport Phenomena, Wiley, Hoboken, NJ.
Koponen, A. , 1998, “ Permeability of Three-Dimensional Random Fiber Webs,” Phys. Rev. Lett., 80(4), pp. 716–719.
Jackson, G. W. , and James, D. F. , 1986, “ The Permeability of Fibrous Porous Media,” Canad. J. Chem. Eng., 64(3), pp. 364–374.
Dulien, F. , 1979, Porous Media Fluid Transport and Pore Structure, Academic Press, New York.
Kubawara, S. , 1959, “ The Forces Experienced by Randomly Distributed Parallel Circular Cylinders of Spheres in a Viscous Flow at Small Reynolds Numbers,” J. Phys. Soc. Jpn., 14(4), pp. 527–532.
Hutten, I. M. , 2007, Handbook of Nonwoven Filter Media, Butterworth-Heinemann, Burlington, MA.
Gebart, B. , 1992, “ Permeability of Unidirectional Reinforcements for RTM,” J. Comp. Mater., 26(8), pp. 1100–1133.
Brinkmann, H. C. , 1949, “ A Calculation of the Viscous Force Exerted by a Flowing Fluid on a Dense Swarm of Particles,” Flow Turb. Combust., 1(1), pp. 27–34.
Nabovati, A. , Edward, W. , and Llewellin, A. C. , 2009, “ A General Model for the Permeability of Fibrous Porous Media Based on Fluid Flow Simulations Using the Lattice Boltzmann Method,” Compos. Part A: Appl. Sci. Manuf., 40(6–7), pp. 860–869.
Delaney, D. , Weaire, S. H. , and Murphy, S. , 2005, “ Random Packing of Elliptical Disks,” Phil. Mag. Lett., 85(2), pp. 89–96.
View article in PDF format.

## References

Bear, J. , 1972, Dynamics of Fluids in Porous Media, Dover Publications, Mineola, NY.
Whitaker, S. , 1986, “ Flow in Porous Media—I: A Theoretical Derivation of Darcy's Law,” Transp. Porous Media, 1(1), pp. 3–25.
Elkatatny, S. , Mahmoud, M. , and Nasr-El-Din, H. A. , 2013, “ Filter Cake Properties of Water-Based Drilling Fluids Under Static and Dynamic Conditions Using Computed Tomography Scan,” ASME J. Energy Resour. Technol., 135(4), p. 042201.
Pilotti, M. , 1998, “ Generation of Realistic Porous Media by Grains Sedimentation,” Transp. Porous Media, 33(3), pp. 257–278.
Maier, R. S. , Kroll, D. M. , Benard, R. S. , Howington, S. E. , Peters, J. F. , and Davis, H. T. , 2000, “ Pore-Scale Simulation of Dispersion,” Phys. Fluids, 12(8), pp. 2065–2079.
Madadi, M. , and Sahimi, M. , 2003, “ Lattice Boltzmann Simulation of Fluid Flow in Fracture Networks With Rough, Self-Affine Surfaces,” Phys. Rev. E, 67, p. 026309.
Zhang, H. , Ge, X. , and Ye, H. , 2006, “ Randomly Mixed Model for Predicting the Effective Thermal Conductivity of Moist Porous Media,” J. Phys. D, 39(1), pp. 220–226.
Wang, M. , Wang, J. , Pan, N. , and Chen, S. , 2007, “ Mesoscopic Predictions of the Effective Thermal Conductivity for Microscale Random Porous Media,” Phys Rev. E, 75(3), p. 036702.
Guo, Z. , and Shu, C. , 2013, Lattice Boltzmann Method and Its Applications in Engineering, World Scientific Publishing, Hackensack, NJ.
He, X. , Zou, Q. , Luo, L. S. , and Dembo, M. , 1997, “ Analytic Solutions of Simple Flow and Analysis of Non-Slip Boundary Conditions for the Lattice Boltzmann BGK Model,” J. Stat. Phys., 87(1), pp. 115–123.
Ladd, A. , 1994, “ Numerical Simulations of Particulate Suspensions Via a Discretized Boltzmann Equation—Part I: Theoretical Foundation,” J. Fluid Mech., 271, pp. 285–310.
Wellein, G. , Lammers, P. , Hager, G. , Donath, S. , and Zeiser, T. , 2006, “ Towards Optimal Performance for Lattice Boltzmann Applications on Terascale Computers,” International Conference Parallel Computational Fluid Dynamics (ParCFD), Indianapolis, IN, May 14–17, pp. 31–40.
He, X. , Doolen, G. , and Clark, T. , 2002, “ Comparison of the Lattice Boltzmann Method and the Artificial Compressibility Method for Navier–Stokes Equations,” J. Comp. Phys., 179(2), pp. 439–451.
Kang, Q. , Lichtner, P. , and Janecky, D. , 2010, “ Lattice Boltzmann Method for Reacting Flows in Porous Media,” Adv. Appl. Math. Mech., 2(5), pp. 545–563.
Murdock, J. , and Yang, S. , 2016, “ Alternate and Explicit Derivation of the Lattice Boltzmann Equation for the Unsteady Incompressible Navier–Stokes Equation,” Int. J. Comput. Eng. Res., 6(12), pp. 47–59.
Hudong, C. , Chen, S. , and Matthaeus, W. , 1992, “ Recovery of the Navier–Stokes Equations Using a Lattice Boltzmann Method,” Phys Rev. A, 45(8), p. R5339. [PubMed]
Yamamoto, K. , and Ohori, S. , 2012, “ Simulations on Flow and Soot Deposition in Diesel Particulate Filters,” Int. J. Engine Res., 14(4), pp. 333–340.
Pan, C. , Luo, L. S. , and Miller, C. , 2006, “ An Evaluation of Lattice Boltzmann Schemes for Porous Medium Flow Simulations,” Comput. Fluids, 35(8–9), pp. 898–909.
d'Humieres, D. , 1992, “ Generalized Lattice-Boltzmann Equations,” Prog. Astro. Aeron., 159, pp. 450–458.
Murdock, J. , Ickes, J. , and Yang, S. , 2017, “ Transition Flow With an Incompressible Lattice Boltzmann Method,” Adv. Appl. Math. Mech., 9(5), pp. 1271–1288.
Kruger, T. , Kusumaatmaja, H. , Kuzmin, A. , Shardt, O. , Silva, G. , and Viggen, E. M. , 2017, The Lattice Boltzmann Method: Principles and Practice, Springer International Publishing, Cham, Switzerland.
Chen, S. , Martinez, D. , and Mei, R. , 1996, “ On Boundary Conditions in Lattice Boltzmann Methods,” Phys. Fluids, 8(9), pp. 2527–2531.
Yuan, P. , and Schaefer, L. , 2005, “ A Thermal Lattice Boltzmann Two-Phase Flow Model and Its Application to Heat Transfer Problems—Part 2: Integration and Validation,” ASME J. Fluids Eng., 128(1), pp. 151–156.
Karani, H. , and Huber, C. , 2015, “ Lattice Boltzmann Formulation for Conjugate Heat Transfer in Heterogeneous Media,” Phys. Rev. E., 91(2), p. 023304.
Inamuro, T. , Yoshino, M. , Inoue, H. , Mizuno, R. , and Ogino, F. , 2002, “ A Lattice Boltzmann Method for a Binary Miscible Fluid Mixture and Its Application to a Heat-Transfer Problem,” J. Comp. Phys., 179(2), pp. 201–215.
Bird, R. , Stewart, W. , and Lightfoot, E. , 2006, Transport Phenomena, Wiley, Hoboken, NJ.
Koponen, A. , 1998, “ Permeability of Three-Dimensional Random Fiber Webs,” Phys. Rev. Lett., 80(4), pp. 716–719.
Jackson, G. W. , and James, D. F. , 1986, “ The Permeability of Fibrous Porous Media,” Canad. J. Chem. Eng., 64(3), pp. 364–374.
Dulien, F. , 1979, Porous Media Fluid Transport and Pore Structure, Academic Press, New York.
Kubawara, S. , 1959, “ The Forces Experienced by Randomly Distributed Parallel Circular Cylinders of Spheres in a Viscous Flow at Small Reynolds Numbers,” J. Phys. Soc. Jpn., 14(4), pp. 527–532.
Hutten, I. M. , 2007, Handbook of Nonwoven Filter Media, Butterworth-Heinemann, Burlington, MA.
Gebart, B. , 1992, “ Permeability of Unidirectional Reinforcements for RTM,” J. Comp. Mater., 26(8), pp. 1100–1133.
Brinkmann, H. C. , 1949, “ A Calculation of the Viscous Force Exerted by a Flowing Fluid on a Dense Swarm of Particles,” Flow Turb. Combust., 1(1), pp. 27–34.
Nabovati, A. , Edward, W. , and Llewellin, A. C. , 2009, “ A General Model for the Permeability of Fibrous Porous Media Based on Fluid Flow Simulations Using the Lattice Boltzmann Method,” Compos. Part A: Appl. Sci. Manuf., 40(6–7), pp. 860–869.
Delaney, D. , Weaire, S. H. , and Murphy, S. , 2005, “ Random Packing of Elliptical Disks,” Phil. Mag. Lett., 85(2), pp. 89–96.

## Figures

Fig. 1

Sample generated substrate geometry: φ = 50%

Fig. 2

Sample generated substrate geometry: φ = 80%

Fig. 3

D2Q9 lattice structure

Fig. 4

Blending functions, W1 (ϕ) and W2 (ϕ) versus porosity (ϕ)

Fig. 5

Computational filter domain and boundaries

Fig. 6

Velocity field at 50% porosity (Re 2.6)

Fig. 7

Velocity field at 80% porosity (Re 49)

Fig. 8

Porosity versus permeability at multiple velocities

Fig. 9

Re versus friction factor at multiple velocities and porosities, compared to Ergun equation

Fig. 10

Nondimensional permeability, K/Dp2 versus ϕ, fit to the Koponen model

Fig. 11

Semilog plot: Nondimensional permeability, K/Dp2 versus ϕ, fit to the Koponen model

Fig. 12

Nondimensional permeability, K/Dp2 versus ϕ, fit to the GB model

Fig. 13

Semilog plot: Nondimensional permeability, K/Dp2 versus ϕ, fit to the GB model

Fig. 14

Nondimensional permeability, K/Dp2 versus ϕ, fit to the GK model

Fig. 15

Semilog plot: Nondimensional permeability, K/Dp2 versus ϕ, fit to the GK model

Fig. 16

Semilog plot: Nondimensional permeability, K/Dp2 versus porosity, ϕ compared with experimental data

## Discussions

Some tools below are only available to our subscribers or users with an online account.

### Related Content

Customize your page view by dragging and repositioning the boxes below.

Related Journal Articles
Related Proceedings Articles
Related eBook Content
Topic Collections