## Abstract

The accuracy of solid wall treatment in the lattice Boltzmann method (LBM) simulation of porous structures affects different hydraulic parameters including integral properties, such as permeability, or local phenomena, such as apparent slip. Based on an analysis of the advantages and disadvantages of the current methods, a new technique is introduced for exact boundary extraction from binary representation. Using this technique, the LBM model can simultaneously benefit from the advantages of existing approaches, i.e. the real micro-/nanostructure obtained with X-ray computed tomography, and a reduction in the resolution requirement. To evaluate the technique, permeability and slip length on the solid walls are investigated for a porous gas diffusion layer. The results show acceptable accuracy improvement balanced with computational costs.

## 1. Introduction

Porous media applications can be found in many different fields of research. Consequently, understanding and predicting fluid flow through structurally complex porous media is of common interest in related areas. Among the numerical methods proposed for porous media fluid-flow simulation during recent years, the lattice Boltzmann method (LBM) has gained the most popularity because of its inherent features. Some of the attributes of this method are its simple numerical algorithm, intrinsic parallelism, its capability for handling multiple fluid species/phases and also non-continuum flows, and finally, and most importantly, the straightforward implementation of complex solid boundaries.

The key characteristic of a porous medium that makes it distinguishable from others is its geometrical microstructure. It has been reported that porous media samples with different microstructures, characterized by the same porosity and tortuosity, show different hydraulic characteristics [1].

Since permeability depends on the microstructure of porous media, the geometry has to be well accounted for at the pore scale, otherwise calculated values may differ drastically from those that are measured, especially if the anisotropic behaviour of permeability is considered. Therefore, one of the most challenging parameters that can be used to evaluate the ability of a numerical method simulating a porous medium to handle complex geometries is the permeability calculation. In addition, some defects in the LBM solid boundary, models, such as apparent slip or viscosity-dependent boundaries, have been reported by several authors [2–5]. Such a dependence could be a good benchmark for studying the capabilities of the proposed technique for implementing the solid boundaries, especially for cases such as porous media in which the wall effects are dominant.

Two different approaches can be found in the literature to enhance the accuracy of solid boundary models in porous media by the LBM. The first is based on an X-ray computed tomography (XCT) imaging technique, hereinafter called the image-based approach. In this approach, capturing the shape of the boundaries is inherently limited to a staircase approximation because of the cubic shape of the voxels that are the smallest accessible units of data in XCT images. A comprehensive literature review addressing advances in combining XCT with computer simulation, including the LBM, for the analysis of porous materials has been presented by Moreno-Atanasio *et al.* [6]; this review also includes a procedure for using binary representation of the real geometry of porous structures as the input for LBM simulation. The limiting factor in this approach is not the spatial resolution of the XCT device but the computational cost of the LBM simulation if the same LBM grid resolution is to be used.

In the second approach, hereinafter called the model-based approach, the geometry is reconstructed by stochastic numerical techniques. Since the real geometry is modelled via a stochastic layout of some known geometric entities such as spheres, ellipsoids and cylinders, the exact profile of a solid boundary is known everywhere. A detailed description and comparison of such methods can be found in Junk & Yang [7]. More specific comparisons of these methods in the field of porous media simulations are presented by Pan *et al.* [5] and Maier & Bernard [8]. By using higher order solid boundary implementation, resolution requirements and convergence times could be considerably reduced without degrading solution accuracy, as reported by Pan *et al.* [5]. This point becomes very important when the computational cost is the main barrier to simulation. However, the requirement for knowing the exact boundary shapes is the main limitation for this approach. All of the literature introducing or using such methods, to the authors’ knowledge, uses the model-based approach for geometry reconstruction, which itself might be considered a source of error owing to deviation from the real geometry. Therefore, in order to benefit from the full potential of the improvement in accuracy of this approach, one should use higher order solid boundary methods to implement the real geometry obtained from XCT images.

In this paper, a new technique is introduced for the exact boundary extraction (EBE) from binary representation by an efficient algorithm and without the need to implement extra data structure in the developed computer code. To evaluate this new technique, the permeability calculation and apparent slip are investigated and the apparent slip length for complex porous media is computed.

## 2. Geometry binary representation

In the image-based approach, the geometry of the flow field is provided for the solver using a three-dimensional binary array in which one (1) corresponds to pore regions and zero (0) to solid ones. Figure 1*a* shows such a data structure schematically.

This three-dimensional binary array is the product of the binarization of a three-dimensional greyscale image obtained from XCT. As the original XCT data are 8-bit greyscale, a threshold value should be selected in order to convert these data to binary format. This threshold value needs to be determined based on an independent estimate of the sample porosity. Ostadi *et al.* [9] showed that the threshold value has a considerable influence on the permeability calculation. In this investigation, the XCT image of a 0.512^{3} mm carbon paper sample scanned at a maximum resolution of 500 nm is used. The X-ray nano-tomography process was carried out using a Phoenix Nanotom S with a 180 kV/15 W ultra-high-performance nano-focus X-ray tube. The captured image is depicted in figure 1*b*. The average pore diameter of this kind of carbon paper with a porosity of 87 per cent is approximately 5 *μ*m, so there are about 10 voxels for each pore diameter, which is more than twice the acceptable resolution (4 voxels/pore diameter), as described by Videla *et al.* [10]. The binary array of size 1024^{3} provided in this step is used as an input file containing geometrical data for LBM simulations. In addition, three other computed tomography (CT) image datasets with resolutions of 512^{3}, 256^{3} and 128^{3} were also prepared. These CT images were used to study the impact of different tomographic resolutions for a given LBM grid.

## 3. Simulation model description

As reported by He *et al.* [2], the lattice Bhatnagar–Gross–Krook model, in which the collision operator is approximated by a single relaxation time, has some defects such as the dependency of boundary locations on viscosity. The viscosity-dependent boundary conditions pose a significant problem for simulating flow through porous media because the permeability becomes viscosity dependent, whereas it should be only a characteristic of the physical properties of the porous medium. This problem can be partially solved by using a multi-relaxation time (MRT)–LBM and will be discussed more in §5*b*. In the current work, a D3Q19 MRT–LBM is used with the set of relaxation parameters proposed by Ginzburg & d’Humieres [3]. The formulation and details of the method are well described by Pan *et al.* [5] and will not be repeated here. The set of relaxation parameters used in this study is as follows:
3.1
The volume-averaged Reynolds number is kept below 0.1 for all simulations and *τ*=1/*S*_{ν}=1.0 in all cases unless otherwise mentioned.

The bounce-back of the non-equilibrium distribution rule introduced by Zou & He [11] is used for the velocity and pressure boundary conditions at the inlet and outlet. The details of implementation for cases similar to the present work can be found elsewhere [12]. The other surfaces of the fluid domain are considered as periodic boundaries.

The fluid–solid boundaries are implemented by three methods, including the standard bounce-back (SBB) scheme, the linearly interpolated bounce-back (LIBB) scheme and the multi-reflection (MR) scheme. The description of these methods and some other higher order methods can be found in several papers [3,5,8,13].

## 4. Exact boundary extraction from binary representation

The key idea of the EBE technique is the extraction of a solid boundary in fine resolution and running the simulation on a coarse grid. This is schematically depicted in figure 2 for a two-dimensional lattice and extension to a three-dimensional lattice is straightforward. All of the higher order solid boundary implementations are based on the points of intersection between the lattice links and the real boundary. Therefore, the only data required to formulate such boundary methods are *q*_{i}s in each lattice node. As can be seen in figure 2, one can calculate *q*_{i}s for a coarse grid based on a fine grid but one can perform the simulation in the coarse grid with a much lower computational cost. The calculation of *q*_{i}s is as simple as counting the 0s and 1s in the direction of the desired lattice link. For example, in figure 2, the base grid is eight times finer than the coarse one (the computational grid for the LBM) and *q*_{i}s can be calculated as *q*_{2}=5/8, *q*_{6}=2/8 and *q*_{3}=4/8.

## 5. Results and discussion

The MRT–LBM model powered by the EBE technique is used to simulate the flow through the complex geometry explained in §2. Two hydraulic parameters, i.e. permeability and slip length, are investigated to evaluate the model capabilities in terms of accuracy and computational costs.

### (a) Permeability calculation

Considering the usage of GDL, permeability is calculated in only the normal direction to the carbon paper sample. This is the same direction as the plane’s normal vector on which the inlet and outlet pressure boundaries are implemented. Therefore, the permeability corresponding to the normal direction of the pressure gradient can be calculated using Darcy’s Law as *k*=−*νρu*_{d}/∇_{n}*P*, where *ν* is the kinematic viscosity, *ρ* is the density, ∇_{n}*P* is the average pressure gradient in the direction of *n* and *u*_{d} is Darcy’s velocity obtained from the volume-averaged fluid velocity over the system [5]. As neither experimental data nor an exact solution exist, the permeability calculated with the highest accessible resolution is considered as the exact reference value of the XCT sample permeability, *k*_{ref}. Since after thresholding no sub-voxel data can be retrieved, the highest possible resolution used by the LBM solver would be the original three-dimensional binary array of size 1024^{3} obtained in the binarization. If *N* is considered as the lattice Boltzmann (LB) grid cell in each dimension and *k*_{N} as the permeability calculated by the LBM using a grid size of *N*, then *k*_{ref}=*k*_{1024}. Since there are no further geometrical data on the solid boundaries’ shape except the staircase approximation, the only method that can be used to calculate *k*_{1024} is the SBB scheme. Of course in the coarser grid, all of the three methods mentioned in §3 can be used via the EBE technique introduced in §4.

As shown in figure 3*a*,*b*, the permeability is calculated in LB grid sizes of 2^{5}, 2^{6}, 2^{7}, 2^{8} and 2^{9}. Normalized permeability, *k*_{N}/*k*_{ref}, for all cases is presented in figure 3*a*. The EBE technique shows considerable improvement in accurately calculating the permeability in both the LIBB and MR schemes when compared with the SBB scheme for the same grid resolution. As shown in figure 3, the best result is obtained by the MR boundary implementation, which has also been reported by other authors [5].

In figure 3*b*, the relative error, *δ*_{k}=|*k*_{N}/*k*_{ref}−1|, of the permeability calculation is shown. In this figure, the capability of the EBE technique to reduce the resolution requirements is clearly shown. As can be seen, if *δ*_{k}<1% is considered as the criterion of the resolution-independent solution, the EBE technique using the MR boundary satisfies this requirement at a resolution of 2^{7}, while the SBB scheme does not reach a resolution-independent solution even at 2^{9}. This means that at least 4^{3} times less resolution is required using the EBE technique. It can also be seen that, at a resolution of 2^{7}, the computational cells (corresponding to the image voxels) are eight times coarser than the resolution of the reference grid (2^{10}), i.e. about 1.25 cells per pore diameter. Therefore, by using the EBE technique, there is no further requirement to set the LB resolution as high as the accepted tomographic resolution (4 voxels/pore diameter). Of course, this reduction in the resolution requirement is gained at the expense of some computational overhead, but it should be noted that most of extra effort for this method is carried out in the pre-processing step and only once in each simulation. One other point to be considered is that the usage of the LIBB and MR schemes also reduces the convergence time, as reported by Pan *et al.* [5].

Table 1 presents the relative error of permeability at different CT image tomographic resolutions, *R*, with the given LBM grid size, *N*, calculated using the MR boundary in all cases. The results imply that the errors are more sensitive to image resolution than the LB grid size. This can be explained by robust boundary reconstruction using the EBE technique, so that reducing the LB grid resolution does not affect the accuracy drastically. On the other hand, using lower tomographic image resolution results in a larger voxel size than the accepted limit, i.e. 4 voxels/pore diameter. Therefore, for an optimized simulation set-up, the lowest acceptable tomographic resolution should be examined with different LB grid sizes to satisfy the criteria for relative error, e.g. *δ*_{k}<1%. In this case, the optimized permeability calculation is obtained at *k*_{128,512}.

### (b) Slip length investigation

It is well known that the solid boundary location in the BGK–LBM depends on the viscosity [2,3]. The MRT–LBM proposed by Ginzburg & d’Humieres [3] can overcome the problem by a carefully constructed collision operator. One can choose the relaxation rates {*S*_{i}}, such that the position of zero velocity is fixed at exactly one-half the lattice spacing beyond the last fluid node for Poiseuille flow when the solid walls are parallel to the underlying lattice lines [5]. This will result in equation (3.1). Obviously, this setting only completely solves the problem (apparent slip) for the solid boundaries parallel to the lattice grid; for other wall configurations new {*S*_{i}} settings are required [14]. Such settings can only be obtained analytically for simple geometries. Therefore, in complex geometries, apparent slip should be expected. In pore-scale simulations in which the flow length scale is of the order of 1 *μ*m and lower, the apparent slip phenomenon becomes more important because of the comparable order of magnitude with physical parameters. Moreover, an accurate understanding of the apparent and physical slip becomes more important in the present study, as the flow is at the onset of the slip regime, i.e. *K*_{n}∼0.01. These points show that the apparent slip is a suitable parameter to evaluate the solid wall implementation technique for this kind of problem.

In this work, the apparent slip is calculated via the slip length, *b*, which is depicted schematically in figure 4*a*. The method of slip length calculation is explained in figure 4*b* and is applied for all boundary nodes with sufficient neighbour fluid nodes for second-order extrapolation and then the slip length is averaged over the whole domain. Figure 5*a* shows the sensitivity of the slip length calculation versus the LB grid size for different solid boundary models. The effects of exact boundary reconstruction on the reduction in apparent slip length is demonstrated well by the results. The effects of the relaxation factor on the slip length are presented in figure 5*b*. As can be seen, by using the MRT–LBM, viscosity dependency remains relatively low even in the SBB boundary model, as reported by Pan *et al.* [5].

## 6. Conclusion

A new computationally efficient approach is developed in the present study for solid boundary implementation in porous complex geometry, represented by binary three-dimensional images obtained with XCT. Using this technique, the LBM model can benefit simultaneously from the real micro-/nanostructure obtained from XCT images and from a reduction in the resolution requirement. The results of the present study show an approximately 4^{3} reduction in the resolution requirement by using the EBE technique for the permeability calculation. The apparent slip resulting from the LBM solid–boundary modelling in complex porous structures is also studied. The slip length investigation shows the capability of the EBE technique to handle this numerical artefact better than the current solid boundary implementation scheme in complex geometries.

## Acknowledgements

The authors acknowledge the kind support from the Vehicle, Fuel and Environment Research Institute (VFERI) of the University of Tehran.

## Footnotes

One contribution of 26 to a Theme Issue ‘Discrete simulation of fluid dynamics: methods’.

- This journal is © 2011 The Royal Society