Simulation of gas flow in micro-porous media with the regularized lattice Boltzmann method

Fuel 205 (2017) 232–246

Contents lists available at ScienceDirect

Fuel journal homepage: www.elsevier.com/locate/fuel

Full Length Article

Simulation of gas flow in micro-porous media with the regularized lattice Boltzmann method Junjian Wang a, Qinjun Kang b, Yuzhu Wang a, Rajesh Pawar b, Sheik S. Rahman a,⇑ a b

School of Petroleum Engineering, University of New South Wales, Sydney, NSW 2033, Australia Earth and Environmental Sciences Division, Los Alamos National Laboratory, Los Alamos, NM 87545, USA

a r t i c l e

i n f o

Article history: Received 3 April 2017 Received in revised form 22 May 2017 Accepted 25 May 2017 Available online 1 June 2017 Keywords: Lattice Boltzmann method Micro-porous media Rarefied gas flow

a b s t r a c t One primary challenge for prediction of gas flow in the unconventional gas reservoir at the pore-scale such as shale and tight gas reservoirs is the geometric complexity of the micro-porous media. In this paper, a regularized multiple-relaxation-time (MRT) lattice Boltzmann (LB) model is applied to analyze gas flow in 2-dimensional micro-porous medium reconstructed by quartet structure generation set (QSGS) on pore-scale. The velocity distribution inside the porous structure is presented and analyzed, and the effects of the porosity and specific surface area on the rarefied gas flow and apparent permeability are examined and investigated. The simulation results indicate that the gas exhibits different flow behaviours at various pressure conditions and the gas permeability is strongly related to the pressure. The increased porosity or the decreased specific surface area leads to the increase of the gas apparent permeability, and the gas flow is more sensitive to the pore morphological properties at low-pressure conditions. Ó 2017 Elsevier Ltd. All rights reserved.

1. Introduction Successful hydrocarbon production of unconventional gas such as shale gas and tight gas in the United States in the last decade has initiated interest around the world. Europe, Australia, and China are now starting to evaluate and explore their unconventional reservoirs [1,2]. To reduce uncertainties and increase the efficiency of production related to these systems, accurate production prediction over time is necessary [3]. This is typically achieved by running large-scale fluid flow simulations on the order of meters to kilometers using constitutive equations. They, therefore, require the input of effective material properties such as porosity and permeability, which characterize the formation where the hydrocarbon production takes place [4]. The unconventional gas reservoirs are composed of various porous components and micro-fractures. Properties of the pore structure and morphology, such as porosity, surface area, connectivity have a significant effect on the formation’s performance. The internal micro-structure and the macroscopic behaviour are indeed strongly and directly related to each other. Therefore, understanding why and how the macroscopic features of the unconventional gas reservoir vary over space and time requires the detail examination of their micro-structure. ⇑ Corresponding author. E-mail address: [email protected] (S.S. Rahman). http://dx.doi.org/10.1016/j.fuel.2017.05.080 0016-2361/Ó 2017 Elsevier Ltd. All rights reserved.

Given the complexity and the heterogeneity of these natural micro-porous media, such a deep understanding is crucial. However, experimental measurements of transport properties in this context are often difficult, expensive and time-consuming since these reservoirs are extremely tight with the average pore space in the matrix on the same order of magnitude as the mean free path of the gas molecules. Therefore, some analytical or semianalytical models that derived based on straight tube models [5–8] and the image-based simulation techniques at the scale of the material’s micro-structure become an attractive tool to supplement direct measurements, and to enhance the understanding of a material’s behaviour. Among the image-based simulation tools, the lattice Boltzmann (LB) method is attracting increasing attention in the recent decade. Generally speaking, both representative elementary volume (REV)-scale and pore-scale LB approaches have been proposed to study the gas transport in micro-porous media. The REV-scale simulation is carried out with the application of generalized lattice Boltzmann model (GLBM) [9], which ignores the detailed structure of the porous media by including an additional term in the LB equation to account for the presence of a porous medium. Although this approach cannot obtain the detailed flow information, it can be used for systems of large size and produces reasonable results with appropriate models for the porous medium. Chen et al. [10,11] combined the Beskok and Karniadakis-Civan’s

J. Wang et al. / Fuel 205 (2017) 232–246

correlation for local apparent permeability with GLBM, and performed series of simulations to discuss the slip effect on gas transport in different reconstructed shale samples. Wang et al. [12] later improved the GLBM of Chen et al. by applying dusty gas modelgeneralized Maxwell Stefan model to study the different flow behaviour of organic matter and inorganic mater of shale gas reservoir. Different from the REV-scale LB simulation, the pore-scale LB direct simulation is the most straightforward way to apply LBM to porous flows based on the detailed geometry. The main advantage of the pore-scale LB simulation is that the local flow information can be obtained which is beneficial to reveal the fluid flow physics and to study the macroscopic relations. Pore-scale LB simulations have be regarded as currently the standard method of choice to calculate the single phase flow and transport in continuum flow regime [13], and many attempts in regarding developing slip boundary conditions and effective viscosity/mean free path for LB equation have extended its application in rarefied gas flows in slip flow and transition flow regimes [14,15]. With the implementation of slip boundary conditions and appropriate expressions of relaxation parameters with effective viscosity or effective mean

233

free path, it has been confirmed that the slip-based LB models can obtain reasonable results for rarefied gas flow in microchannels compared with direct simulation Monte Carlo (DSMC) results or the analytical solution of linear Boltzmann equation. However, in spite of the efforts having been conducted to reveal the rarefied gas flow behaviour with LB model, the coupled effects of rarefaction and gas–solid interaction on the gas flow in microporous media at the pore scale are still less understood. To date, little attention has been paid to the influence of stochastic geometric morphologies of micro-porous media on gas flow on the pore scale. Especially, how the detailed nature of pore structure for porous media affects the gas flow behaviour through a porous medium is yet to be explored. The challenges of applying slip-based LB models to simulate gas flow in micro-porous media include but are not limited to: 1) The irregular boundary of solid obstacles increases the complicity of the application of slip boundary condition. 2) Since Kn is a geometry-dependent parameter, the tortuous connected pore spaces and the confinement of the random distributed solid obstacles raise the difficulties for its calculation. To find the local Kn for gas molecules flow in the micro-porous media, Zhao et al. [16,17] applied a medial axis method and suggested to

Fig. 1. Illustration of the 4-CCL algorithm. Background pixels in black are the solid phase. Object pixels in white are the void phase.

234

J. Wang et al. / Fuel 205 (2017) 232–246

(a) Reconstructed structure based on QSGS

(b) Structure after removing the occluded pores

Fig. 2. Comparison of structures before and after removing the occluded pores.

use the distance to the molecule’s closest solids or the average distance to the surrounding solids as the characteristic length. Instead, Landry et al. [18] defined an averaged exponential wall function to calculate the local mean free path. In this study, in order to examine and investigate the effect of rarefaction and pore configuration of the porous structures on gas flow, and to solve the mentioned two problems, we present a novel adaptation of the regularized multiple-relaxation-time (MRT)-LBM for gas flow in 2-dimensional (2D) reconstructed micro-porous media. The rest of the paper is as follows. The reconstruction of the micro-porous media is introduced in Section 2, and the mathematical model for predicting gas flow is introduced in Section 3. In Section 4, LB model is validated with previous results for gas flow in the micro-channel with/without obstacles, and the gas flow in reconstructed micro-porous media is analyzed in Section 5. Some conclusions are drawn in Section 6.

2. Shale structure reconstruction with improved quartet structure generation set(QSGS) Very recently the developing characterization techniques such as FIB-SEM and Micro-CT provide new insights into the microstructures of organic shale. It is found that the pore structure of organic shale is very complex and the morphological properties of shale present multi-scale and heterogeneous features [19]. The porosity of organic shale varies in a wide range at micro-scale. Loucks et al. [20] reported the porosity of Barnett shale is 0–30%, while Sondergeld et al. [21] reported a porosity of 55% in kerogen based on 2-D SEM images. Also, the micro-structure of shale matrix is of high heterogeneity and diversity. The characterization of organic shale based on FIB-SEM shows that the interconnected micro-pores with irregular shapes in different components provide the main pathways for gas flow in shale matrix [22]. In order to generate porous morphological features closely resembling the porous shale matrix, the quartet structure generation set (QSGS) method proposed by Wang et al. [23] is adopted in the present study. Since the porous shale matrix consists of solid and pore phases, correspondingly, a two-phase QSGS process is utilized, and the solid and the pores are set as the growing phase (marked with ‘‘0”) and the non-growing phase (marked with ‘‘1”), respectively. The QSGS process is described as follows: (a) Randomly distribute seeds in the domain based on a distribution probability, cd , and the value of the cd is far less than its desired volume fraction. (b) Based on a given directional growth probability Di , the solid phase grows from each seed along its eight directions. To the each growing element, new random numbers are assigned to its neighbouring cells. (c) The process (b) is repeated until the volume fraction of the growing phase reaches the desired value. The spaces not occupied at the end represent the porosity /.

Although QSGS has been widely used to generate the structure of the porous media, this method generates occluded pores which are not connected to the main void space and do not contribute to flow. Considering flow in a porous medium, only the interconnected pore spaces from the inlet to the outlet are of interest, a 4-connected component labelling (4-CCL) algorithm [24] is applied to identify inter-connected pores and occluded pores. The algorithm is described in the following steps: (a) Classify the pixels of the porous structure into two types: object pixel (void phase) and background pixel (solid phase) (Fig. 1(a)). (b) First pass labelling: Scan the image and label each object pixel with a temporary number by checking its east, west, south and north neighbourhoods. If none of the neighbours is labelled, assign a new label to this pixel and set the label as the pixel’s class. If one neighbour is labelled, assign the label and the class of this neighbour to the current pixel. If two or more neighbours are already labelled, assign the current pixel the minimum of its neighbours’ label values, and set the class of the minimal neighbour as the class of this pixel and other labelled neighbours (Fig. 1(b)). (c) Second pass labelling: Replace each temporary label by the smallest label of its equivalence class (Fig. 1(c)). (d) After a unique value for each inter-connected void phase is defined, set the occluded pores as solids and remain the void space that connects the inlets and outlets, then update the volume fraction or the porosity / (Fig. 1(d)). A comparison of the structure generated by QSGS before and after the application of component labelling algorithm is shown in Fig. 2. Although the reconstructed model may be considered somewhat artificial, it nevertheless captures much of the geometric complexity of the natural micro-porous media such as the shale micro-porous structure. 3. Lattice Boltzmann model 3.1. LB equation The multiple-relaxation-time (MRT)-LB equation derived from Boltzmann equation can be written as [25]:

f i ðx þ ci dt ; t þ dt Þ  f i ðx; tÞ ¼ Xi ðf Þ;

ð1Þ

where f i ðx; tÞ is the discrete distribution function for particles with velocities ci at position x and time t; dt is the time step, and Xi ðf Þ is the MRT discrete collision operator. In D2Q9 model, the discrete velocity ci are defined as c1 ¼ ð0; 0Þ; c2 ¼ c4 ¼ ð1; 0Þc; c3 ¼ c5 ¼ ð0; 1Þc; c6 ¼ c8 ¼ ð1; 1Þc; c7 ¼ c9 ¼ ð1; 1Þc. where c ¼ dx =dt is the lattice speed with dx being the lattice spacing. The MRT collision operator is given by:

235

J. Wang et al. / Fuel 205 (2017) 232–246

Xi ðf Þ ¼ ðM1 SMÞij ðf j  f eq j Þ; where

M

is

a

ð2Þ

orthogonal

transformation

matrix

and

1

S ¼ diagðsq ; se ; se ; sj ; sq ; sj ; sq ; ss ; ss Þ is a diagonal matrix. The distribution function f i and its equilibrium distribution eq function f i can be projected onto the moment space, and the following results can be obtained: T

m ¼ Mf ¼ ðq; e; e; jx ; qx ; jy ; qy ; pxx ; pxy Þ ;

ð3Þ

3.2. Boundary conditions 3.2.1. Wall boundaries In order to capture the gas slippage near the solid boundary, the kinetic boundary condition proposed by Tang et al. [29] is introduced in this study, and the distribution functions of the boundary nodes f i ðxÞ are a combination of the specular reflection distribution sr functions f i with the completely diffusive distribution functions dr

f i which can be expressed as: dr

m ¼ Mf eq

eq

¼ ðq; e ; e eq

eq

¼ qð1; 2 þ 3juj2 ; 1  3juj2 ; u; v ; u; v ; u2  v 2 ; uv Þ ; T

T

eq

sr

f i ðxÞ ¼ rf i ðxÞ þ ð1  rÞf i ðxÞ; ðci  uw Þ  n > 0;

eq eq eq T ; jx ; qeq x ; jy ; qy ; pxx ; pxy Þ

ð4Þ

where the specularly reflected distribution function sr

f i ðxÞ ¼ f i0 ðxÞ;

eq eq eq T ðf 1 ; f 2 ; . . . ; f 9 Þ ;

ð9Þ sr fi

satisfies:

ð10Þ

where f ¼ ðf 1 ; f 2 ; . . . ; f 9 Þ and f ¼ q is the density, e is the energy mode, e is related to energy square, ðjx ; jy Þ are the momentum components, ðqx ; qy Þ correspond to energy flux, ðpxx ; pxy Þ are related to the diagonal and off-diagonal components of the stress tensors; u and v are the components of the velocity u, respectively. Therefore, the MRT-LB equation can be written as:

with ci0 ¼ ci  2½ðci  uw Þ  nn being the specular velocity of ci , and n is the unit normal vector. The diffusive distribution part is derived by Ansumali and Karlin [30] and can be calculated as:

j fðx þ ci dt; t þ dtÞi j fðx; tÞi ¼ M1 S½j mðx; tÞi  meq ðx; tÞi:

To match the second order slip boundary condition with a halfway reflection mode, the portion of the diffusive reflection part in the combination (r) in Eq. (9) and the relaxation time sq must be chosen as [31]:

ð5Þ

For isothermal gas flow, the macroscopic density, momentum, and the shear viscosity are recovered by:

q¼

X

f i;

qu ¼

i

X ci f i ;



l ¼ qc2s ss 

i

 1 dt ; 2

ð6Þ

pffiffiffi pffiffiffiffiffiffi where cs ¼ RT ¼ c= 3 is the sound speed, in which R is the gas constant and T is the temperature. For gas flow in micro-porous tight reservoir, one important characteristic parameter of the flow is the Knudsen number. Thus, the viscosity or relaxation parameter(s) in the LB models and the Knudsen number/mean free path should be appropriately related for the simulation of micro-scale gas flows. Based on gas kinetic theory, the relaxation time ss can be given as [26]:

ss

1 ¼ þ 2

rffiffiffiffi 6 k ; p dx

ð7Þ

1 2

rffiffiffiffi rffiffiffiffi 6 kuðdi ; kÞ 1 6 kloc ¼ þ ; 2 p dx p dx

ð8Þ

where uðdi ; kÞ is the correction function which considers the effect of wall surface on the mean free path k; di is the distance to the wall and kloc is the local mean free path.

ðcj uw Þn<0 jðcj

ðci uw Þn>0 jðci

 uw Þ  njf j ðxÞ

1 2

sq ¼ þ

eq

 uw Þ  njf i ðxÞ

rffiffiffiffi 1  p r ¼2 1þ ; A1 6

pA2 ð2ss  1Þ2 þ 3 ; 8ð2ss  1Þ

eq

f i ðxÞ:

ð11Þ

ð12Þ

ð13Þ

where A1 and A2 are functions of tangential momentum accommodation function, r [26]:

A1 ¼ A2 ¼

where k is the mean free path. In micro-scale wall bounded geometries, the collisions between gas molecules are not sufficient due to the wall confinement effects, especially in the near wall region. Thus, there exists a so-called Knudsen layer near the wall in which the local mean free path of the total molecules is smaller than that in the unbounded systems due to the wall effects. To account for the effect of gas molecule/wall interactions, generally two kinds of corrections on mean free path and/or viscosity are proposed in the literature. Guo et al. [25] and Li et al. [27] averaged the effect of the Knudsen layer to the whole flowing area and proposed the simplified effective Knudsen number expressions. Besides, Stops [28] implemented a ‘‘geometry dependent” correction function to account for the wall confinement. More discussions on Knudsen layer effect can be found from our previous review paper [14]. In this study, the correction function approach is used, and the local relaxation time ss satisfies:

ss ¼ þ

P dr

f i ðxÞ ¼ P

2r

r 2

ð1  0:1817rÞ;

1 þ A21 :

p 2

ð14Þ ð15Þ

From Eqs. (12) and (13) it can be seen that different from the SRT-LB model, in which the control parameter r in the slip boundary condition depends on the local Knudsen number, the interaction parameters A1 and A2 , and the lattice size (see Eq. (26) in Ref. [32]), with the application of MRT collision operator, the control parameter r in the LB boundary condition depends on only the gas–solid interaction parameters A1 and A2 , which eliminates the difficulties of calculating the local characteristic length and the local Knudsen number of gas molecules flow in complex structures. The application of slip boundary condition of LBM (Eq. (9)) requires the identification of the normal direction to the interfacial node. The approach for boundary identification and classification in porous media is described as follows: (a) Read into the data information of the phase functions of the porous media generated by the QSGS method. (b) Enlarge the geometry by the multiple of 4. This treatment is conducive to excluding the fewest number of unexpected boundary nodes in following steps. Meanwhile, it is contributive to improving the grid resolution and the computational accuracy. (a) Identify the fluid boundary nodes. If a fluid boundary node has at least one solid phase node as its adjacent point, then this fluid node is marked as a boundary node.

236

J. Wang et al. / Fuel 205 (2017) 232–246

f i ðNx þ 1; yÞ ¼ 2f i ðNx ; yÞ  f i ðNx  1; yÞ;

i ¼ f4; 7; 8g:

ð17Þ

After this, in order to realize the pressure boundary condition, the extrapolated density needs to be rescaled to the prescribed density:

Ny qin

qnew ð1; yÞ ¼ PNy

y¼1

qð1; yÞ

qð1; yÞ;

N qout

qnew ðNx ; yÞ ¼ PNy y y¼1

qð1; yÞ

ð18Þ

qðNx ; yÞ:

ð19Þ

3.3. Regularization procedure

Fig. 3. Illustration of the 12 types of boundary nodes, namely, four flat wall nodes (F1)–(F4), four convex corner nodes (C1)–(C4), and four concave corner nodes (C5)– (C8). The blocks in grey, downward diagonal and black are solid, fluid boundary nodes and fluid, respectively.

Generally speaking, the LBM is carried out in alternate streaming and collision steps. When this dynamic procedure is viewed as a projection of the LB equation onto the Hermite space HN , it introduces an error analogous to the aliasing effect because the distribution function f i does not automatically lie entirely within HN . Such an error is usually very small and negligible. However, it can be no longer neglected when the system is away from equilibrium. Therefore, a ‘‘regularization procedure” is introduced to resolve this error and to improve the numerical stability [34–36]. The procedure is implemented as follows: First, the post-streaming distribution function f i is divided into eq the equilibrium part f i and the non-equilibrium part f i 0: eq

f i ¼ f i þ f 0i : As

eq fi

ð20Þ

already lies entirely in the subspace H , then f i 0 needs to N

be converted to a new distribution ~f 0i which lies within the subspace spanned by the higher-order Hermite polynomials. Using the second-order Hermite polynomials, for the D2Q 9 models, ~f 0 is i

expressed as:

" ~f 0 ¼ w i i

Fig. 4. A molecule confined between two planar walls with spacing H [37]. The gas molecule has an equal probability to travel in any zenith angle hþ or h and r is the travelling distance.

#  Q1 1 ð2Þ ci X 0 ~ f cj cj : H 2c2s cs j¼0 j

ð21Þ

where Hð2Þ is the second order Hermite polynomial, and is given ð1Þ

ð2Þ

explicitly as Hð0Þ ðcÞ ¼ 1; Hi ðcÞ ¼ ci and Hij ðcÞ ¼ ci cj  dij , with dij 0 is the Kronecker delta function. By replacing f with ~f 0 , the i

i

MRT-LB model becomes: (d) Classify all the fluid boundary nodes into 12 types, i.e., four flat wall nodes, four convex corner nodes and four concave corner nodes. Fig. 3 presents an illustration of the 12 types of fluid boundary nodes. Flat wall nodes ðF1—F4Þ are within a straight-line boundary. Convex nodes ðC1—C4Þ and concave nodes ðC5—C8Þ are observed at the points of intersection between two straight-line boundaries.

eq j fðx þ ci dt; t þ dtÞi ¼j f ðx; tÞiþ j ~f 0 ðx; tÞi  M1 SM j ~f 0 ðx; tÞi   S M j Fidt; ð22Þ þ M1 I  2

and this regularization procedure filters out the higher-order non-equilibrium moments that contain strong discrete artefacts. 3.4. Local mean free path

Specifically, this approach is not only applicable to reconstructed images, but also applicable to the scanned images whose solid and void phases are marked differently. 3.2.2. Inlet/outlet In order to specify the pressure on flow boundaries, the extrapolated boundary conditions [33] are employed at inlet/outlet ends of the channel. The inlet and outlet are located at x ¼ ð1; yÞ and x ¼ ðN x ; yÞ with y 2 ½1; N y , the ghost nodes are defined at the positions x ¼ ð0; yÞ and x ¼ ðN x þ 1; yÞ. The following equations are used to determine the unknown distribution based on a linear extrapolation scheme:

f i ð0; yÞ ¼ 2f i ð1; yÞ  f i ð2; yÞ;

i ¼ f2; 6; 9g;

ð16Þ

In micro-scale wall bounded geometries, through the wall function approach, the local mean free path is a function of the mean free path of gas molecules in the unbounded system and the distance of these molecules to the solid walls. Some studies on gas flow in simple micro-scale geometries such as channels (2D) and tubes(3D) considering the local mean free path of gas molecules have been proposed in the literature. However, to our best knowledge, the wall function approach to porous media remains a challenge. Dongari et al. [37] proposed a power-law (PL) based correction function and carefully validated their data against molecular dynamics (MD) data for a two planar wall case (Fig. 4), and they claim that their model can capture the Knudsen layer in relative high Knudsen numbers compared with Arlemark et al.’s

237

J. Wang et al. / Fuel 205 (2017) 232–246

kloc;planar ¼

kþ ðyÞ þ k ðyÞ kðuðdþ ; kÞ þ uðd ; kÞÞ ¼ ; 2 2

ð23Þ

with dþ and d are the distances to the upper and lower walls and wall function, u satisfies:

" 1n 1n 8  X 1 dm dm uðdm ; kÞ ¼ 1  þ4 1þ 1þ 48 a acos½ð2i  1Þp=32 i¼1 #   1n 7 X dm 1þ þ2 ; ð24Þ acos½i p=16 i¼1

Fig. 5. A molecule confined by solid obstacles in porous media. The gas molecule has an equal probability to travel in any zenith angle h and r is the travelling distance.

model [38] and Stop’s model [28]. For gas molecules confined in the space shown in Fig. 4, the local mean free path given by Dongari et al. [37] is the average of the mean free path confined by each wall:

where m indicates the wall confinement from + or  direction, a ¼ kðn  2Þ and n ¼ 3 which has been tested by Ref. [39] and is applied in this study. For gas flow in reconstructed porous media, considering that the gas molecules are confined by solid walls in different directions and can move towards either direction with the same probability (Fig. 5), in this study, the wall confinement in xþ; x; yþ and y directions are considered and the local mean free path for all the molecules in the flow domain can be determined by averaging these four parts [37]:

kloc ¼

4 4 1X 1 X kloc;m ¼ k uðdm ; kÞ; 4 m¼1 4 m¼1

(a) Knout =0.0194

(b) Knout =0.194

(c) Knout =0.388 Fig. 6. Normalized streamwise velocity at different Knudsen number.

ð25Þ

238

J. Wang et al. / Fuel 205 (2017) 232–246

the relaxation time ss , the MRT-LB model combined with powerlaw function to calculate the mean free path applied in this study allows us to consider the wall confinement from each direction in details with the ease of eliminating the calculation of characteristic length.

4. Validation

Fig. 7. Schematic graph of the geometry of the flow around a square cylinder in a microchannel [43].

where m indicates xþ; x; yþ or y direction. From Eq. (24) and (25), it is easy to see that, the PL model applied in this study satisfies the physically intuitive requirements for Kn ! 0, i.e., ujwall  1=2 and ujbulk  1, and it is consistent with the phenomenological viscosity model derived by Fichman and Hetsroni [40]. To calculate the distance of each lattice to the solid boundaries dm , a window growth method is proposed and applied. The implementation is as follows: (a) Read into the data information of the phase functions of the porous media generated by the QSGS method. The inlet/outlet boundaries are set to‘‘3”, which is different from solid phase and void phase. (b) Identify the fluid nodes, set the fluid node as the center of a square window. Starts growing the window from the centre to find solid nodes in xþ; x; yþ and y directions. If a solid boundary node is found, its corresponding distance to the centre is recorded, and the window stops growing in this direction. Otherwise, moves to the next step. (c) If the radius meets the inlet/outlet, an extremely large distance is given to this direction, to indicate the gas molecules are not bounded in this direction. (d) Repeat process (b) and (c) until the distances of every fluid nodes to solids or inlet/outlet boundaries in four directions have been found. The relaxation time ss of each lattice node can be obtained by substituting Eq. (24) and (25) to Eq. (8). Different from the most recent methodologies such as the use of the closest distance to the solids [16] or the average distance to the solids [17] to calculate

(a) Knout =0.0194

The first validation is carried out by gas flow in a single channel, and the results of the streamwise velocity profile are compared with that of DSMC and IP-DSMC results as reported by Shen et al. [41], which are suitable for modelling micro-gaseous flow and have been extensively used for such validation purposes. The LB simulation results without regularization procedure are also presented for comparison. In this case, the ratio of the length to the height L=H is 100, and a uniform lattice N x  N y ¼ 2100  21is adopted. The outlet gas density qout is 1, and the pressure ratio of the inlet and outlet pin =pout is taken as 1.4 for Knout ¼ 0:0194 and 2 for Knout = 0.194 and Knout = 0.388, respectively. Extrapolated boundary conditions are applied to the inlet and outlet of the channel, and the proposed boundary condition is applied to the solid wall. The relaxation times are specified as follows: ss is determined by Eq. (8) with k ¼ Knout H to simulate the gas flow at the disired Knout ; sq is given by Eq. (13), se ¼ 1:1; se ¼ 1:2, and sq ¼ sj ¼ 1 [32]. The tangential momentum accommodation coefficient r is set to be 1 [27,42]. The convergence criterion towards the steady state is set as follows:

sP ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 i juðxi ; t n Þ  uðxi ; t n  1000dtÞj < 1012 P 2 juðx ; t Þj i n i

ð26Þ

The normalized streamwise velocity U ¼ u=umax predicted by this study for Knout ¼ 0:0194; Knout ¼ 0:194 and Knout ¼ 0:388 are shown in Fig. 6. It can be seen that the velocity profile of nonregularized and regularized LBMs are both in good agreement with those of the DSMC/IP-DSMC at Kn ¼ 0:0194 and Kn ¼ 0:194, which indicates that the aliasing error of the distribution function f is negligible for slip flow in parallel wall confined channel at low Kn. However, with the increase of Kn, the system is away from the equilibrium and the aliasing error increases. When Kn reaches 0.388, the velocity profile of the regularized LB model shows a better match with the DSMC/IP-DSMC results compared with that of the non-regularized LB model. In order to verify further the applicability of the current model to simulate gas flow in complex geometries, gas flow around a square cylinder at Kn ¼ 0:084 in a micro-channel [34] is considered. The geometry of the simulation domain in shown in Fig. 7. In the LB simulation, N x  N y ¼ 100  100 is adopted to eliminate

(b) Knout =0.194

(c) Knout =0.388

Fig. 8. Distribution of the streamwise velocity magnitude at different Kn with non-regularized LBM.

239

J. Wang et al. / Fuel 205 (2017) 232–246

(a) Knout =0.0194

(b) Knout =0.194

(c) Knout =0.388

Fig. 9. Distribution of the streamwise velocity magnitude at different locations with regularized LBM.

Fig. 10. Validation of normalized streamwise velocity profiles at Kn ¼ 0:084. Dash-line is the regularized LB simulation result, solid-line is the non-regularized LB simulation result. Dots are MD simulation results.

Fig. 11. The 2-D porous structure generated with QSGS.

the grid dependency, and the pressure boundary condition is applied to the inlet and outlet boundaries to set the nondimensional acceleration value as 0.001. The contour maps of the normalized streamwise velocity with and without regularization at different Kn are presented. The streamwise velocity profiles are also compared with molecular dynamics (MD) simulations. The distributions of the normalized streamwise velocity magnitude at different Kn are shown in Fig. 8 and Fig. 9. It clearly shows that, even though the LB model without the regularization works reasonably well at relative low Kn in the single channel case. When

obstacles are introduced in the system, it presents a non-physical pattern of the flow field around the square cylinder at the same Kn. Also, this pattern becomes more obvious with the increase of Kn (From Fig. 8(a) to (c)). As can be seen in Fig. 9, the regularized LB model filters out this momentum oscillation, and the flow field becomes continuous at a wide range of Kn. As discussed in Section 3, without the regularization procedure, the collision step in the LB model introduces an error because the distribution functions do not automatically lie entirely within the Hermite space. This error is amplified by the complex structure of the porous

240

J. Wang et al. / Fuel 205 (2017) 232–246

Fig. 12. Results of the velocity vector profile for water and methane flow in the porous structure. P out ¼ 20 MPa; rx P ¼ 0:1 MPa=m; Width ¼ 100 nm and Length ¼ 400nm.

(a) x/ L = 0.2

(b) x/ L = 0.4

(c) x/ L = 0.6

(d) x/ L = 0.8

Fig. 13. Velocity profile of methane and water at different locations.

media compared to simple geometry. Therefore, the performance of the majority of the LB models proposed for channel flow in simulating gas flow in complex micro-geometry could be questionable, and the regularisation procedure is required to improve the accuracy of the model. Fig. 10 presents the non-dimensional streamwise velocity normalized by the bulk mean velocity at x=H ¼0, 0.25, 0.5 and 0.75. Compare with the predictions of molecular dynamic (MD) simulations, the solution of velocity without the regularization becomes oscillated and has non-physical sharp kinks. However, the general velocity profiles are well predicted by the regularized D2Q9 LB model at different locations in the simulation domain, which

confirms that the regularized LB model can be applied to study the gas flow in complex geometries. 5. Results and discussion In this section, the gas flow in reconstructed shale samples at different pressure are discussed in details. The extrapolated boundary conditions are applied to the inlet and outlet boundaries and the slip boundary condition is applied to the solid wall. The relaxation times and the TMAC are the same as that applied eq in Section 4. Before starting the simulation, f with u = 0 and q satisfying linear distribution from qin and qout are adopted as

241

J. Wang et al. / Fuel 205 (2017) 232–246

Fig. 14. Results of the velocity vector profile for methane at different pressure condition in the porous structure. rx P ¼ 0:1MPa=m; Width ¼ 100nm and Length ¼ 400nm.

(a) x/ L = 0.2

(c) x/ L = 0.6

(b) x/ L = 0.4

(d) x/ L = 0.8

Fig. 15. Distribution of the streamwise velocity magnitude at different Kn with regularized LBM.

the initial value of distribution f i at each lattice. The conversion from lattice unit to physical unit is given in the following way. The unit of length is given by the lattice size dx, and the unit of time dt is derived from the kinematic viscosity. Considering

the kinematic viscosity has the dimensions of length squared over time, the velocity in physical unit satisfies uphy ¼ udx2 =dt. The scaling for other parameters such as velocity and permeability can be obtained in the same way.

242

J. Wang et al. / Fuel 205 (2017) 232–246

Fig. 16. Variation of the apparent permeability K app with outlet pressure.

Table 1 Structure parameters for the 2-D reconstructed micro-porous media. Structure

G1-1

G1-2

G1-3

G1-4

G2-1

G2-2

G2-3

G2-4

/ S

0.64 0.38

0.64 0.26

0.64 0.16

0.64 0.11

0.37 0.10

0.41 0.10

0.43 0.10

0.46 0.10

5.1. Velocity distribution The velocity distribution of the porous media provides the understanding of the fluid flow physics in the confined space. In order to gain a deeper insight into the gas flow behaviour in unconventional gas reservoir, the methane gas flow in a 2-D structure is studied, and the water flow in the same structure is also presented for comparison, since the slip effect is not important for liquid flow in micro-scale porous system as the mean free path of fluid molecules is still very small. The structure of organic shale is generated with QSGS which is introduced in Section 2. The irregular solid obstacles are randomly distributed in the reconstructed structure to reflect the geometric complexity of the shale samples. In order to maintain the connectivity of the 2-D reconstructed structure, a

relatively high porosity is assumed. The 4001  1001 lattice is adopted for LB simulation to represent a porous structure with Length  width = 400 nm  100 nm (See Fig. 11). The porosity of the structure is 0.66, the temperature of the system is set as 338.7 K, the outlet pressure is 20 MPa, and the pressure gradient rP is 0:1MPa=m. The mean free path of methane at the outlet is 0.307 nm. The density, viscosity and the molecular weight of methane and water are calculated based on the Peace software [44]. The velocity distribution of water and methane in the 2-D reconstructed structure is shown in Fig. 12, and the velocity profile of water and methane at x=L = 0.2, 0.4, 0.6, 0.8 is present in Fig. 13, respectively. As expected, the existence of the solid obstacles obviously disturbs the velocity distribution inside the porous media. The streamlines are distorted, especially for the regions adjacent to the solid obstacles. This phenomenon indicates that the presence of the obstacles suppresses the flow and contributes to extra energy loss to the system. Fig. 12 also shows that methane can easily flow through some narrow throats compared to water as a result of the rarefaction effect, and these narrow throats connect with each other to provide extra pathways for methane flow (see locations marked with red circles and dashline in Fig. 12(b)). These further lead to different velocity distributions for methane and water flows. Take the velocity profile at x=L ¼ 0:4 in Fig. 13 as an example, because methane prefers to flow through the void space marked by red dashline in Fig. 12(b), the velocity of methane is larger than that of water at y=H = 0.25 to 0.35, and is smaller than water at y=H ¼ 0:5 to 0. The gas slippage also results in higher overall methane velocity along the streamwise direction compared to water. In addition, as the pressure decreases from the inlet to the outlet, the general velocity profile of methane increases markedly toward the exit by the constraint of mass conservation. The similar trend of the general velocity profile has also been observed in gas micro-channel flows [25].

5.2. Influence of rarefaction Based on the results of Section 5.1, it can be concluded that unlike the classical continuum flow, the rarefaction effect typically must be taken into account for gas flow in micro-porous media since the molecular free path is comparable to the system’s charac-

Fig. 17. Results of the velocity vector profile for methane flow in micro-porous structure with different specific surface area. / = 0.64, P out = 10 MPa.

J. Wang et al. / Fuel 205 (2017) 232–246

243

Fig. 18. Results of the velocity vector profile for methane flow in the micro-porous structure with the different specific surface area. / = 0.64, P out = 30 MPa.

Fig. 19. Results of the velocity vector profile for methane flow in the micro-porous structure with the different specific surface area. / = 0.64, P out = 50 MPa.

teristic dimensions. In this context, to investigate the rarefaction effect on the extraction of gas from the unconventional gas reservoir, the methane flow in the pore structure shown in Fig. 11 at a wide range of outlet pressure from 10 MPa to 50 MPa and a constant pressure gradient rP ¼ 0:1  106 Pa=m is studied. The temperature of the system is set as 338.7 K, the density, viscosity and the molecular weight of methane at different pressure are calculated based on the Peace software [44]. Fig. 14 presents the velocity vector distribution of methane at different pressure conditions, and the velocity profile at x=L = 0.2, 0.4, 0.6 and 0.8 is presents in Fig. 15, respectively. As shown in this figure, because the gas molecules collide with the solid surface more frequently than they collide with each other, the rarefaction effect is considered to be important for gas flow in micro-porous media. As a result, the non-equilibrium distribution of gas molecular velocities is induced and gives rise to the boundary velocity slippage phenomenon. This becomes evident at low pressure or

Fig. 20. Results of the apparent permeability for methane flow in the micro-porous structure with the different specific surface area. / = 0.64.

244

J. Wang et al. / Fuel 205 (2017) 232–246

Fig. 21. Results of the velocity vector profile for methane flow in the micro-porous structure with the different specific surface area. S = 0.10, P out = 10 MPa.

Fig. 22. Results of the velocity vector profile for methane flow in the micro-porous structure with the different specific surface area. S = 0.10, P out = 30 MPa.

high Kn, and the bulk velocity increases as the pressure decrease, especially near the outlet. In order to quantitatively evaluate the gas flow in micro-porous structure as affected by the rarefaction effect, the apparent permeability K app which equals 2uP out lL=ðP2in  P 2out Þ is investigated. According to the definition, the higher the apparent permeability, the more significant effect the rarefaction induces. Fig. 16 shows the variation of K app as a function of the outlet pressure. As it can be seen, because the decreases in outlet pressure with a constant pressure gradient lead to a large Kn, there is a trend towards the increase of the K app as the outlet pressure decreases, and the increase is faster at the low-pressure range. Therefore, the gas permeability must be considered as a dynamic reservoir parameter and updated accordingly as the reservoir is being depleted.

5.3. Influence of pore morphology In the confined micro-space, dispersed solids result in large surface-to-volume ratios, causing wall effects to be of significant importance in gas flows. For gas flow in micro-porous media, the surface-gas interaction by its very nature is highly dependent on the porous structure and its morphology. In this section, two important morphological properties: the porosity / and the specific surface area S which is the ratio between the total interstitial surface area and the bulk volume (the total interfacial length and bulk area for the 2D case), are discussed. Two groups of microporous structures are reconstructed: one group has a different surface area with the same porosity, whereas the other group has different porosity with the identical specific surface area. The

J. Wang et al. / Fuel 205 (2017) 232–246

245

Fig. 23. Results of the velocity vector profile for methane flow in the micro-porous structure with the different specific surface area. S = 0.10, P out = 50 MPa.

Fig. 24. Results of the apparent permeability for methane flow in the micro-porous structure with different porosity. S = 0.10.

decreases, the gas-surface interaction becomes increasingly significant for gas flow in micro-structures. At low pressure, the rarefaction and the gas-surface interaction starts to dominate the gas flow. Figs. 21–23 show the velocity vector profile of structure G2-1 to G2-4 which have different porosity but the same specific surface area. With an increased porosity, the velocity vector length increases and the velocity field becomes more intense under different pressure conditions. The high porosity attributes to the low resistance for mass transfer, which is associated with a faster flow rate. The corresponding apparent permeability is shown in Fig. 24. Clearly, a large porosity gives rise to a large apparent permeability. Similar to the variation of permeability with the specific surface area, the K app  / relationship varies with pressure, and K app is more sensitive to porosity at low pressure condition. This phenomenon implies that the impact of porosity on K app becomes more significant at low-pressure conditions with enhanced rarefaction.

6. Conclusion morphological features of the reconstructed micro-porous structures are listed in Table 1. The length  width is 400 nm  100 nm, and the 4001  1001 lattice is adopted for the LB simulation. The methane flow in micro-structure under P out ¼ 10 MPa; 30 MPa and 50 MPa with rP ¼ 0:1 MPa=m is considered and the properties of methane are calculated based on the Peace software [44]. Figs. 17–19 present the velocity vector profile of structure G1-1 to G1-4 which have a different specific surface area with the same porosity of 0.64 at different pressure conditions. At a certain porosity, the increase of the specific surface area increases the tortuosity of the micro-structure. As a result, the pore structure with smaller particle size and high specific surface area exerts more resistance to the gas flow and more likely to block or eliminate the flow passages. Fig. 20 quantitatively presents the apparent permeability as a function of the specific surface area at different pressures. As shown, similar to the intrinsic permeability, the apparent permeability also decreases monotonously with increasing specific surface area. In addition, the apparent permeability is more sensitive to the specific surface area at low-pressure conditions. This behaviour can be explained by the fact that, as the pressure

In this study, a regularized MRT lattice Boltzmann model is adopted to simulate the gas flow in micro-porous media. The quartet structure generation set (QSGS) method is applied for structure reconstruction, and the method is improved by introducing a 4connected component labelling algorithm to eliminate the occluded pores. The power-law (PL) wall function [37,39] for channel flow is modified and extended to porous media, to calculate the relaxation parameters for MRT-LB model. With more relaxation parameters, the current LB model eliminates the difficulties in calculating local characteristic length, and the regularization procedure improves the velocity oscillation around the solid obstacles. Based on the simulation results, it can be concluded that the gas shows a large velocity slippage at the micro-porous media, which leads to a faster mass transfer. The velocity slippage at the solid boundary becomes more obvious with the pressure decrease, which suggests that the apparent permeability must be considered as a dynamic parameter. The sensitivity studies of gas flow in different reconstructed micro-porous structure indicate that the pore morphology plays

246

J. Wang et al. / Fuel 205 (2017) 232–246

a significant role in controlling the gas flow behaviours. The morphological effect becomes significant in micro-scale gas flow at the low-pressure condition with enhanced rarefaction. In particular, the porosity and specific surface area are of considerable importance to the gas flow in the micro-porous media owing to the dominance of the surface-gas interaction. Acknowledgements The authors acknowledge the support from SCOPE, UNSW. J.W. would like to acknowledge the financial support from the China Scholarship Council (CSC) and the Postdoctoral Writing Fellowship, UNSW. Q.K. and R. P. would like to acknowledge the support from a DOE oil & gas project. J.W. also thank Prof. Kazuhiko Suga from Osaka Prefecture University, Japan, for his helpful discussions. References [1] Yuan J, Luo D, Feng L. A review of the technical and economic evaluation techniques for shale gas development. Appl Energy 2015;148:49–65. [2] Weijermars R. Us shale gas production outlook based on well roll-out rate scenarios. Appl Energy 2014;124:283–97. [3] Kim TH, Cho J, Lee KS. Evaluation of co 2 injection in shale gas reservoirs with multi-component transport and geomechanical effects. Appl Energy 2017;190:1195–206. [4] Middleton RS, Carey JW, Currier RP, Hyman JD, Kang Q, Karra S, JiménezMartínez J, Porter ML, Viswanathan HS. Shale gas and non-aqueous fracturing fluids: opportunities and challenges for supercritical co 2. Appl Energy 2015;147:500–9. [5] Javadpour F, Fisher D, Unsworth M. Nanoscale gas flow in shale gas sediments. J. Can. Pet. Technol. 2007;46(10):55–61. [6] Xu J, Wu K, Yang S, Cao J, Chen Z, Pan Y, et al. Real gas transport in tapered noncircular nanopores of shale rocks. AIChE J 2017. [7] Wu K, Chen Z, Li X, Dong X. Methane storage in nanoporous material at supercritical temperature over a wide range of pressures. Sci Rep 2016;6. [8] Wu K, Li X, Wang C, Yu W, Chen Z. Model for surface diffusion of adsorbed gas in nanopores of shale gas reservoirs. Ind Eng Chem Res 2015;54(12):3225–36. [9] Guo Z, Zhao T. Lattice boltzmann model for incompressible flows through porous media. Phys Rev E 2002;66(3):036304. [10] Chen L, Fang W, Kang Q, Hyman JD, Viswanathan HS, Tao W-Q. Generalized lattice boltzmann model for flow through tight porous media with klinkenberg’s effect. Phys Rev E 2015;91(3):033004. [11] Chen L, Kang Q, Dai Z, Viswanathan HS, Tao W. Permeability prediction of shale matrix reconstructed using the elementary building block model. Fuel 2015;160:346–56. [12] Wang J, Chen L, Kang Q, Rahman SS. Apparent permeability prediction of organic shale with generalized lattice boltzmann model considering surface diffusion effect. Fuel 2016;181:478–90. [13] Blunt MJ, Bijeljic B, Dong H, Gharbi O, Iglauer S, Mostaghimi P, Paluszny A, Pentland C. Pore-scale imaging and modelling. Adv Water Resour 2013;51:197–216. [14] Wang J, Chen L, Kang Q, Rahman SS. The lattice boltzmann method for isothermal micro-gaseous flow and its application in shale gas flow: a review. Int J Heat Mass Transf 2016;95:94–108. [15] Wang J, Kang Q, Chen L, Rahman SS. Pore-scale lattice boltzmann simulation of micro-gaseous flow considering surface diffusion effect. Int J Coal Geol 2017;169:62–73. [16] Zhao J, Yao J, Li A, Zhang M, Zhang L, Yang Y, Sun H. Simulation of microscale gas flow in heterogeneous porous media based on the lattice boltzmann method. J Appl Phys 2016;120(8):084306. [17] Zhao J, Yao J, Zhang M, Zhang L, Yang Y, Sun H, An S, Li A. Study of gas flow characteristics in tight porous media with a microscale lattice boltzmann model. Sci Rep 2016;6.

[18] Landry CJ, Prodanovic´ M, Eichhubl P. Direct simulation of supercritical gas flow in complex nanoporous media and prediction of apparent permeability. Int J Coal Geol 2016;159:120–34. [19] Tahmasebi P, Javadpour F, Sahimi M. Multiscale and multiresolution modeling of shales and their flow and morphological properties. Sci Rep 2015;5. [20] Loucks RG, Reed RM, Ruppel SC, Jarvie DM. Morphology, genesis, and distribution of nanometer-scale pores in siliceous mudstones of the mississippian barnett shale. J Sediment Res 2009;79(12):848–61. [21] Sondergeld CH, Ambrose RJ, Rai CS, Moncrieff J, et al. Micro-structural studies of gas shales. In: SPE Unconventional Gas Conference, Society of Petroleum Engineers; 2010. [22] Schieber J, et al. Common themes in the formation and preservation of intrinsic porosity in shales and mudstones-illustrated with examples across the phanerozoic. In: SPE Unconventional Gas Conference, Society of Petroleum Engineers; 2010. [23] Wang M, Wang J, Pan N, Chen S. Mesoscopic predictions of the effective thermal conductivity for microscale random porous media. Phys Rev E 2007;75(3):036702. [24] Haralock RM, Shapiro LG. Computer and robot vision. Addison-Wesley Longman Publishing Co., Inc; 1991. [25] Guo Z, Zhao T, Shi Y. Physical symmetry, spatial accuracy, and relaxation time of the lattice boltzmann equation for microgas flows. J Appl Phys 2006;99 (7):074903. [26] Guo Z, Zheng C, Shi B. Lattice boltzmann equation with multiple effective relaxation times for gaseous microscale flow. Phys Rev E 2008;77(3):036707. [27] Li Q, He Y, Tang G, Tao W. Lattice boltzmann modeling of microchannel flows in the transition flow regime. Microfluidics Nanofluidics 2011;10(3):607–18. [28] Stops D. The mean free path of gas molecules in the transition regime. J Phys D: Appl Phys 1970;3(5):685–96. [29] Tang G, Tao W, He Y. Lattice boltzmann method for gaseous microflows using kinetic theory boundary conditions. Phys Fluids (1994-present) 2005;17 (5):058101. [30] Ansumali S, Karlin IV. Kinetic boundary conditions in the lattice boltzmann method. Phys Rev E 2002;66(2):026311. [31] Guo Z, Shi B, Zhao T-S, Zheng C. Discrete effects on boundary conditions for the lattice boltzmann equation in simulating microscale gas flows. Phys Rev E 2007;76(5):056704. [32] Guo Z, Zheng C. Analysis of lattice boltzmann equation for microscale gas flows: relaxation times, boundary conditions and the knudsen layer. Int J Comput Fluid Dyn 2008;22(7):465–73. [33] Verhaeghe F, Luo L-S, Blanpain B. Lattice boltzmann modeling of microchannel flow in slip flow regime. J Comput Phys 2009;228(1):147–57. [34] Suga K, Takenaka S, Ito T, Kaneda M, Kinjo T, Hyodo S. Evaluation of a lattice boltzmann method in a complex nanoflow. Phys Rev E 2010;82(1):016701. [35] Zhang R, Shan X, Chen H. Efficient kinetic method for fluid simulation beyond the Navier–Stokes equation. Phys Rev E 2006;74(4):046703. [36] Latt J, Chopard B. Lattice boltzmann method with regularized pre-collision distribution functions. Math Comput Simul 2006;72(2):165–8. [37] Dongari N, Zhang Y, Reese JM. Modeling of knudsen layer effects in micro/nanoscale gas flows. J Fluids Eng 2011;133(7):071101. [38] Arlemark EJ, Dadzie SK, Reese JM. An extension to the navier–stokes equations to incorporate gas molecular collisions with boundaries. J Heat Transfer 2010;132(4):041006. [39] Dongari N, Agrawal A. Modeling of Navier–Stokes equations for high knudsen number gas flows. Int J Heat Mass Transf 2012;55(15):4352–8. [40] Fichman M, Hetsroni G. Viscosity and slip velocity in gas flow in microchannels. Phys Fluids (1994-present) 2005;17(12):123102. [41] Shen C, Tian D-B, Xie C, Fan J. Examination of the lbm in simulation of microchannel flow in transitional regime. Microscale Thermophys Eng 2004;8 (4):423–32. [42] Zhuo C, Zhong C. Filter-matrix lattice boltzmann model for microchannel gas flows. Phys Rev E 2013;88(5):053311. [43] Suga K. Lattice boltzmann methods for complex micro-flows: applicability and limitations for practical applications. Fluid Dyn Res 2013;45(3):034501. [44] B. Wischnewski, Online calculation of properties of water and steam, http:// www.peacesoftware.de/einigewerte/einigewerte.html.