Permeability correlation with porosity and Knudsen number for rarefied gas flow in Sierpinski carpets

Accepted Manuscript Permeability correlation with porosity and Knudsen number for rarefied gas flow in Sierpinski carpets H. Rostamzadeh, M.R. Salimi, M. Taeibi-Rahni PII:

S1875-5100(18)30292-0

DOI:

10.1016/j.jngse.2018.06.037

Reference:

JNGSE 2636

To appear in:

Journal of Natural Gas Science and Engineering

Received Date: 4 February 2018 Revised Date:

12 May 2018

Accepted Date: 22 June 2018

Please cite this article as: Rostamzadeh, H., Salimi, M.R., Taeibi-Rahni, M., Permeability correlation with porosity and Knudsen number for rarefied gas flow in Sierpinski carpets, Journal of Natural Gas Science & Engineering (2018), doi: 10.1016/j.jngse.2018.06.037. This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

ACCEPTED MANUSCRIPT

Permeability correlation with porosity and Knudsen number for

RI PT

rarefied gas flow in Sierpinski carpets

H. Rostamzadeh1, M.R. Salimi2*, M. Taeibi-Rahni3

Master of Science, Department of Aerospace Engineering, Sharif University of Technology, Azadi Ave, Tehran, Iran,

SC

1

2

M AN U

E-mail: [email protected]

*Corresponding Author, Assistant Prof., Aerospace Research Institute, 15th St., Mahestan St., Shahrak Ghods, Tehran, Iran, P.O. Box: 14665-834, E-mail: [email protected]

3

Professor, Department of Aerospace Engineering, Sharif University of Technology, Azadi Ave, Tehran, Iran,

EP

Abstract

TE D

E-mail: [email protected]

In recent years, application of porous media is highlighted among researchers due to their wide

AC C

range of usability in micro-scale problems, such as gas reservoirs, micro-filtering, heat exchangers, etc. With this respect, the accurate description of flow behavior using governing equations based on the continuum assumption is not valid since the mean free path is comparable to the characteristics length of the problem. For this purpose, a simple methodology for diffusion reflection boundary condition is developed and validated for two valuable benchmarks, namely micro-channel flow and fractal porous media, where the results were in good agreement with literature. Then, pore-scale simulation of the rarefied gas flow inside the square- and circular1

ACCEPTED MANUSCRIPT

based Sierpinski carpets in the slip and transition flow regimes (based on the developed boundary condition) for different porosities are carried out, using lattice Boltzmann method (LBM). In each case, the effects of different key parameters such as pressure gradient in the x-

RI PT

direction (body force), porosity, and Knudsen number on the velocity contours, velocity profiles, flow structures, and permeability are investigated. It is found that the velocity magnitude is reduced as Knudsen number increases, where the reduction rate of velocity in the slip flow

SC

regime is greater than that of the transition one. Also, a correlation between the permeability with porosity and Knudsen number in the Darcy flow regime ( Re << 1 ) is presented. For both

exponentially with the porosity.

M AN U

geometries, it is observed that the permeability varies logarithmically with Knudsen number and

Keywords: Sierpinski carpet; Rarefied gas flow; Lattice Boltzmann method (LBM); Pore-

TE D

scale simulation; Diffusion reflection; Nomenclature

Symbols

Re

Reynolds number

u

Velocity

Total pores area

At

Total area

wα

Weight coefficient

ci

Desceret velocity

xb

Boundary node location

xf

Fluid node location

Df fα

AC C

cs

EP

Ap

Sonic velocity Fractal dimension

Distribution function Equilibrium distribution function Post-collision distribution function

Greek Symbols

F

Body force

α

Fractal coefficient factor

Fx

Body force in x-direction

Λ

Relaxation matrix

H

Height

λ

Mean free path

2

K

Permeability

λmax

Upper limit self-similar region

Kn

Knudsen number

λmin

Lower limit self-similar region

L

Length

τ

Relaxation time

Lch

Charactristics length

ε

Porosity

LBM

Lattice Boltzmann method

ρ

Density

M

Transformation matrix

ρ0

m

Momentum vector

ω

meq

Equilibrium momentum vector

µ

MRT

Multiple relaxation time

υ

RI PT

ACCEPTED MANUSCRIPT

N

Cumulative number of capillary tubes

δt

Time increment

Ny

Number of nodes at y-direction

δx

x increment

P

Pressure

Reference density

Relaxation parameter Dynamic viscosity

M AN U

SC

Kinematic viscosity

1. Introduction

TE D

Nowadays, porous media provide a wide range of applications in engineering sciences to improve the performance of systems. Many industrial processes in chemical and metallurgy engineering involve fluid flow through the packed bed, which ultimately leads to an increase in

EP

the fluid-solid contact or mixing processes. Many applicable examples can be pointed out as follows: systems including catalytic and chromatographic reactions, packed absorption and

AC C

distillation towers, ion exchange columns, packed filters, pebble-type heat exchangers, oil reservoirs, geothermal processes, etc. The design process of these systems occurs with the mechanisms of pressure drop, fluid flow, and mass and energy transfer equations in the packed beds.

Meanwhile, real porous media have far more complex and irregular structures than other porous media. This has made it difficult to simulate the exact flow characteristic for these complex

3

ACCEPTED MANUSCRIPT

geometries[1]. Fortunately, in recent years, various studies have been carried out on the structures of various micro-porous media such as soil, sandstone, bones, industrial foams, wood, etc., which imply that all these porous media have fractal properties. Due to these complexities,

RI PT

an exact solution of the internal flows without simplification and prediction of permeability (using certain control parameters) will not be possible to achieve accurate results[2]. The utilized

geometries due to the use of empirical coefficients.

SC

mathematical models are not usually able to predict the permeability of these complex

M AN U

In 1967, the well-known French mathematician Mandelbrot presented a concept of fractal which is then resulted in the theory of fractal geometry[3]. Fractals are objects that look similar when viewed at different scales. Such objects have details at small scales that make their expression very complex in Euclidean spaces[4]. Many natural phenomena such as topographic surfaces, jagged shorelines, fickle changing clouds, and seismic activity follow this space. With this

TE D

respect, Cao et al.[5] used the fractal measure to interpret lunar fractures using qualitative (similar ratio) and quantitative (fractal dimension) approaches based on the data generated from the Lunar Reconnaissance Orbiter. They concluded that most of the lunar surface exhibits fractal

EP

behavior. Many studies in this regard indicate that pore spaces in porous media in nature have

AC C

fractal characteristics [6-8]. Also, various fractal methods have been considered by researchers for describing fluid flow behavior and predicting porous media permeability in recent decades. For example, Katz and Thompson[9] were the first people who studied in this field and showed that the pore spaces of sandstones are fractal and self-similar, where their length scales range from 100 µm to 10Å. They also showed that the volume of fractal pores between the gravels is the same. Yu et al.[10] employed the Monte Carlo method to calculate the pores diameter and permeability of the fractal porous media. Their results showed that the conducted simulation has

4

ACCEPTED MANUSCRIPT

an acceptable accuracy with respect to the analytical solutions for many problems. Also, the simulation of this group is able to predict other transport properties (such as thermal conductivity, diffusion conduction coefficient, and electrical conductivity coefficient) for fractal

RI PT

porous media in both saturated and unsaturated states. Fractal analysis is also widely used in fracture studies and morphological-textural characteristics of porous structures in the preparation of porous media[11, 12], industrial activities[13, 14], and other research works[15, 16]. For

SC

example, Othman and Helwani[17] developed a fractal permeability model for bi-dispersed porous media based on the pores fractal characteristics of the membranes. They introduced the

M AN U

utilized fractal permeability model as a suitable method for describing the observed data for quantitative analysis of the bi-dispersed porous media model.

Assumptions of the Navier-Stokes equations are only valid when three following conditions of the Newtonian condition (fluid with Newtonian viscosity), continuity approximation (the mean

TE D

free path is substantially less than the characteristic length), and thermodynamic equilibrium are met[18]. The gaseous flow in micro- and nano-scales tubes that violates these conditions are modeled based on the rarefied gas dynamics. In this case, the effects of rarefied gas flow are

EP

taken into account by Knudsen number which is defined as the ratio of the mean free path ( ) to ):

AC C

the characteristic length of the tube (

=

(1)

Rarefaction effect is pronounced at high Knudsen numbers, which can be due to an increase in the mean free path (e.g., at low pressures) or to a decrease of the characteristic length (e.g., at the micro- and nano-scales tubes) or both of them[19]. In general, various flow regimes existed in literature can be classified as follows[19, 20]: 5

ACCEPTED MANUSCRIPT < 10 ): In this regime, continuum and thermodynamic

1) Full continuum flow regime (

equilibrium assumptions are valid, and hence the flow behavior can be described by NavierStokes equations under implementation of the no-slip boundary conditions.

< 10 ): In this regime, continuum assumptions are valid

RI PT

<

2) Slip flow regime (10

within the fluid bulk, while near walls the non-equilibrium conditions dominate. 3) Transition flow region (10

<

< 10): In this regime, the fluid flow characteristics are

4) Free molecular flow regime (

SC

between slip and free molecular flows.

> 10): In this regime, the non-equilibrium effects dominate

M AN U

nearly through the all bulk fluid flow. In this case, the inter-molecular collisions are dominated by the collisions between the gas molecules and wall surfaces. > 10 , the flow behavior cannot be predicted

According to the above classification, when

by Navier-Stokes equations with no-slip boundary conditions, and hence other methods such as

TE D

molecular dynamic (MD), direct simulation of Monte Carlo (DSMC), or lattice Boltzmann method (LBM) must be used in order to describe flow characteristics more accurately. To design an optimum system for micro-scale problems, computational efficiency and simulation accuracy

EP

should be taken into account, simultaneously. Due to the widespread accessibility of the Navier-

AC C

Stokes equations, many attempts have been made to apply slip models in Navier-Stokes equations for Knudsen numbers exceeding 0.01. One possible continuous solution is solving Burnett-type equations which include second order terms[21]. The main problem in such cases is the implementation of the complex structural rules for high-order correlations. For example, the results of Bailey et al.[22] showed that the results of Navier-Stokes equations with slip boundary condition in Knudsen numbers less than 0.05 differ considerably with those of the experimental. However, the kinetic models provide a more accurate prediction. In engineering applications, for 6

ACCEPTED MANUSCRIPT

example, the macroscopic flow quantities such as shear stresses, wall slip velocities, and mass flow rates are important parameters to be determined. Therefore, methods such as MD and

RI PT

DSMC are not practical in terms of computational cost when no microscopic details are needed. On the other hand, LBM can be an interesting idea due to its inherent kinetic advantage for micro-systems, in which both macroscopic and microscopic characteristics are of particular importance. LBM has a big advantage over the Navier-Stokes solvers at high Knudsen numbers

SC

which can provide accurate results in a fluid flow beyond slip flow regime[21]. Also, the results

M AN U

have shown that LBM covers the Navier-Stokes and Burnett equations and is very suitable for the gaseous flow in the micro-systems which rarefaction effects are important at the surface of the walls[22]. The most important point is to apply an appropriate slip boundary condition in this method. With respect to this issue, the fundamental study of Nie et al. can be pointed out, in which the half-way bounce-back boundary condition is used to directly calculate the slip velocity

TE D

on the solid-fluid surface without using the slip boundary condition[23]. Later, Succi in 2002 [24] introduced specular reflection boundary condition to accurately simulate the fluid flow inside micro-channels, which was consistent with results of the molecular dynamics. Meanwhile,

EP

fluid flow simulation in the microelectromechanical systems using LBM is widely applicable,

AC C

where the results of these simulations indicate that LBM is consistent with analytical and experimental solutions[25-28]. Nevertheless, the results of studies shown that the spatial accuracy of the LBM for simulating gaseous flow in the transient flow regime is not appropriate and the solid constraints (walls) exert nonlinear effects on the relaxation time that should be considered during the simulation[29, 30]. These constraints are due to the growth of the Knudsen layer near the walls, which leads to a non-linear relation between the shear stress and the strain rate within this layer. To overcome these nonlinear effects, Zhang et al.[31] presented an

7

ACCEPTED MANUSCRIPT

effective mean free path to cover the effects of the Knudsen layer more accurately. Their results showed that in this case, the accuracy of the LBM near the walls is improved, and hence this numerical method shows highly economical results for different ranges of Knudsen numbers.

RI PT

Also, Isfahani and Soleimani[32] modified the dynamic viscosity and thermal conductivity coefficient and provided a new relationship for effective relaxation time which was able to simulate the fluid flow characteristics for different Knudsen numbers.

SC

To the best of the authors’ knowledge and by surveying the above mentioned and other similar

M AN U

literature reviews, no general relationship between permeability with porosity and Knudsen number in Sierpinski carpet for the rarefied gaseous flow in the Darcy regime is presented in a more simple form. However, this point is very important to be taken into account by simulating square- and circular-based Sierpinski carpets at different porosities and Knudsen numbers. This motivates the present study to cover shortcomings of the existing literature. In the present study,

•

TE D

the main objectives are multi-fold which can be summarized as follows: To develop a novel methodology in implementation of diffusion reflection boundary

•

EP

condition for rarefied gaseous flow in the porous media. To employ the proposed diffusion reflection boundary condition in the micro-channel

•

AC C

flow at different Knudsen numbers. To investigate the effects of rarefaction effects for two fractal geometries of the squareand circular-based Sierpinski carpets at different porosities and Knudsen numbers based on the developed boundary condition.

•

To present a general correlation between permeability with porosity and Knudsen number for both square- and circular-based Sierpinski carpets at Darcy regime.

2. Fractal models 8

ACCEPTED MANUSCRIPT

In the traditional model of the capillary tube bundle, the porous medium is usually considered as an ideal simplified model that consists of several parallel capillary tubes with a constant diameter, and hence pores size of the porous medium is homogeneous. However, for porous

RI PT

media such as gas or oil reservoirs, pore structure is very complex and the pores size is not the same. As a result, use of the traditional model of capillary tube bundle will result in some restrictions. Therefore, it is necessary to define the fractal characteristic for this type of porous

SC

media with microscopic pore structures. With this regard, the process of generating various stages of Sierpinski carpets with fractal geometries based on the square and circular shapes is

M AN U

shown in Figs. 1 and 2, respectively. In these figures, the white and black areas represent the pores sub-area and impermeable substance, respectively.

The cumulative number of capillary tubes for the unit cross-sectional area whose radii are greater than the pore radius is calculated from the following equation[33, 34]: (> ) =

TE D

where, α and

, (2)

are the factor of the fractal coefficient and fractal dimension of fractal

geometries, respectively. For two-dimensional fractal geometries 2 <

AC C

EP

dimensional fractal geometries 1 <

< 2.

< 3, while for three-

Assuming unit length for each impermeable substance, the total pore area for the first stage of Sierpinski carpet (Figs. 1(a) and 2(a)) can be expressed as follows[34]: !" = #

9

$%& $'(

)

, (3)

ACCEPTED MANUSCRIPT

where,

$%& (=

3) and

$'( (=

1) are the upper and lower limits of self-similar regions,

respectively. The upper and lower limits of self-similar regions for the second stage are 1 and 9,

The total area for these Sierpinski carpets can be expressed as[34]: $%& $'(

) . (4)

SC

!* = #

) and porosity (-) relations for 2D Sierpinski

Using Eqs. (3) and (4), the fractal dimension (

!" =# !*

$'(

$%&

)

M AN U

carpet can be expressed as: -=

RI PT

while for the third stage are 1 and 27, respectively.

,

.

=2−

-

$'(

$%&

. (5)

Eq. (5) shows that the fractal dimension for Sierpinski carpet is function of porosity and upper

TE D

and lower limits of self-similar regions. The value of

for the square- and circular-based

Sierpinski carpets is 1.8927 and 1.9168, respectively.

EP

Porosity of the nth stage of square- and circular-based Sierpinski carpets can be expressed

AC C

respectively as[35]:

8 ( -( = # ) , (6) 9

4 8 ( 51 − # ) 6. (7) -( = 1 − 9 4

To calculate the permeability of the porous media (in Darcy regime) following equation is used[33]:

10

ACCEPTED MANUSCRIPT = %(

is the mean flow velocity and >& is implemented body force through x-direction.

RI PT

where, 9:

89: %( 89: %( = , (8) ∆< >& ∆=

The characteristics length scale for the fractal geometries can be given as[36]:

√ . (9) -(

SC

=

M AN U

3. Lattice Boltzmann method

In BGK (Bhatnagar-Gross-Krook) model, the evolution of particle distribution function can be written as[37]: (= + A BC, C + BC) −

1 F DE

(=, C) −

(=, C)G + BC #1 −

1 ) > , (10) 2DE

represents volume-averaged single particle distribution function along

TE D

where,

(=, C) =

direction and

is the corresponding equilibrium particle distribution function. In addition, t is time, x refers

to a position in space, A stands for particle velocity in

direction, DE is relaxation time which is

AC C

direction.

EP

related to kinematic viscosity by H = IJ (DE − 0.5)BC, and > is external force component along

For the MRT-LBM scheme, the time evolution equation can be expressed as follows[37]: (= + A BC, C + BC) −

(=, C) = −K L FMN (=, C) − MOP (=, C)G, (11) N

11

ACCEPTED MANUSCRIPT where, MN and MN

OP

are the momentum vector and its corresponding equilibrium momentum

vector, respectively. Q is the identity matrix, K is the transformation matrix and K

is its

1 2 1

SC

−1 −1

1 2  1  1 1 ,  − 1 − 1  0 − 1

−1

−1 0

(12)

M AN U

M D2Q9

1 1 1 1 1 1 1 − 4 − 1 − 1 − 1 − 1 2 2   4 −2 −2 −2 −2 1 1  1 0 −1 0 1 −1 0  = 0 −2 0 2 0 1 −1  0 1 0 −1 1 1 0 0 0 −2 0 2 1 1  1 −1 1 −1 0 0 0 0 0 0 0 0 1 −1 

RI PT

inverse form. The transport matrix for D2Q9 model can be defined as:

1

In Eq. (11), L is relaxation matrix which is a non-negative diagonal matrix consisting of nine

relaxation rates R . For D2Q9 model, the relaxation matrix can be defined as[37]:

TE D

L = STUV(RW , R , R , RX , RY , RZ , R[ , R\ , R] ) (13)

In Eq. (13), RW , RX , and RZ are set to zero to ensure mass and momentum conservations[37]. Also,

EP

R\ and R] are set to 1/DE in order to reproduce the same viscosity at that of the SRT (single

relaxation time) model. Other relaxation rates have no physical meaning for incompressible

AC C

flows and are usually selected from linear stability analysis in the range of 0 < R' < 2 [37].

Meanwhile, for D2Q9 model, the equilibrium momentum vector (MOP ) can also be written

as[37]:

MOP = _`, −2` + 3(a& + ab )/`W , ` − 3(a& + ab )/`W , a& , −a& , ab , −ab , (a&

− ab )/`W , (a& ab )/`W c (14)

where, a& = `d& and ab = `db .

12

ACCEPTED MANUSCRIPT

Also, the external forcing term and equilibrium distribution function in Eq. (10) can be expressed respectively as: 1 1 1 A . d + Y (A . d ) − d. dg (15) 2IJ 2IJ IJ

> = `e

RI PT

= `e f1 +

A .> (16) IJ

the weight coefficient in

SC

where, IJ is the speed of sound related to streaming speed (h = B=/BC) by IJ = h/√3, and e is direction, defined for D2Q9 model as:

M AN U

4 l = 0, j9 1 e = = 1 − 4, (17) k9 j 1 i 36 = 5 − 8.

TE D

Macroscopic flow quantities can be found, using calculated distribution functions as: `=m

EP

`d = m A

=m = mA

, (18) . (19)

AC C

In order to simulate rarefied gas flow with LBM, a correlation between Knudsen number and relaxation time is required. With this respect, Zhang et al.[31] presented the following relation as:

8 (D − 0.5) = n , (20) 34 b

13

ACCEPTED MANUSCRIPT

where,

b

= o/Bp is the number of lattice sites at y-direction (for a micro-channel, o is the

height of the channel). Kn is Knudsen number which is defined as Eq. (1).

RI PT

As mentioned earlier, the gaseous flow restricted to a solid surface (wall) is exposed to constant and uniform shear stresses. In this case, describing the gaseous flow behavior by Eq. (20) in LBM at higher rarefaction effects will result in significant inaccuracy. To overcome this

SC

dilemma, the effective mean free path can be defined to decrease effects of Knudsen layer at high Knudsen numbers. This effect can be redirected to the viscosity (or relaxation time) by defining effective viscosity or effective relaxation time. In this case, we have[28, 38]: ), (21)

M AN U

D ∗ = Dr(

s ∗ = sr(

), (22)

where, D ∗ and s ∗ are the effective relaxation time and effective viscosity, respectively. In

TE D

addition, r function can be defined as[28]: r(

)=

2 CU 4

_√2

X/Y

c (23)

EP

4. Boundary conditions

AC C

In this section, the utilized boundary conditions through the numerical simulation are discussed. These boundary conditions are bounce-back, diffusion reflection, and periodic boundary conditions. In this study, bounce-back boundary condition is used in no-slip flow regime, whereas developed diffusion reflection boundary condition is used in slip and transient flow regimes. These boundary conditions are explained in the following sub-sections. 4.1.Bounce-back boundary condition

14

ACCEPTED MANUSCRIPT

Gas particles are reflected after colliding with solid surfaces. For fluids with low dilution effect, applying the bounce-back boundary condition is the most common method for this reflecting illustration. In this case, the collided particles to the solid surface are bounced back in the

boundary node (=t ) can be expressed as[39]:

u̅ (=t , C)

RI PT

reversed direction (Fig. 3(a)). In this case, the distribution functions after the collision at the

= ' _= , Cc, (24)

SC

where, Iu̅ = −I' , =t is the node position at boundary wall, and = is the fluid node near =t

M AN U

(= = =t − I' BC).

Bounce-back boundary condition presented above can be easily implemented to complex geometries like porous media due to its local property by using LBM in comparison with the classical computation fluid dynamic (CFD) methods.

TE D

4.2.Diffusion reflection boundary condition If the wall is assumed to be extremely rough at the molecular scale, the gas particles will reflect

EP

at some random angle, which are uncorrelated with their incident angle. In this case, the velocities of the reflected particles follow the Maxwellian distribution law, which is called

AC C

diffusion reflection boundary condition (Fig. 3(b))[18, 40-43]. In this study, wall temperature is equal to the fluid temperature, and hence has no significant effect on the fluid flow simulation[40]. To get particles reflected diffusively at time t in the bottom wall node point M (Fig. 4), it must be notified that this point receives particles originating from the fluid nodes as well as ghost nodes. The particles’ distribution functions mixed together and then redirected towards the same nodes due to the discrete characteristics of

15

ACCEPTED MANUSCRIPT

space velocities which are not allowed to follow other directions. Information on the diffusively reflected particles from point M towards fluid nodes is contained in the value of

Z

,

and

[

(post-collision) in their corresponding directions. Due to the mixing process that occurs at point ,

Z

and

[

(the unknown post-collision

RI PT

M, the diffusively reflected distribution functions

quantities) should be dependent on the incoming distribution functions of

Y,

\

and

].

Therefore, conservation of diffusively reflected particle number in point M (for a straight wall)

+

Z

=

[

+

Y

\

+

] . (25)

M AN U

+

SC

requires[31]:

In addition, these diffusively reflected distribution functions ( , Maxwellian distribution as follows[40]:

where,

,

Z

and

[

Z

Z

=

TE D

=

[

[

Z

and

[)

must follow the

, (26)

are equilibrium distribution functions at their corresponding direction

EP

and can be obtained from Eq. (15).

AC C

By solving Eqs. (25) and (26), the unknown distribution functions can be obtained as follows: =

Z

=

[

=

+ + +

Z Z Z [ Z

+

[

+

[

+

[

(

Y

+

\

+

] ). (27)

(

Y

+

\

+

] ), (28)

(

Y

+

\

+

] ). (29)

The similar procedure can be carried out the upper wall of a channel, too. 16

ACCEPTED MANUSCRIPT

Eqs. (25)-(29) were driven for straight walls which can be used in whole domain of a microchannel flow. However, for micro-porous media, which are composed of straight, convex, and concave walls, implementation of boundary condition is rather different and difficult. The most

RI PT

conventional method with high accuracy can be curved boundary conditions. However, no slip boundary conditions for the curved surfaces are presented up to now. One alternative way to account the slip effects on the curved surfaces can be a step-by-step curved surface formation

SC

which resembles behavior of a real curved surface with relatively low accuracy. The accuracy of this method can be increased by increasing the number of solid obstacles of the porous media. In

M AN U

this study, the step-by-step curved surface is used. Regarding these points, the smallest square and circular solid blocks used in Sierpinski carpets are depicted in figures 5 and 6, respectively. Figure 5 illustrates the boundary nodes at the convex wall of the smallest solid obstacle for D2Q9 model. The mixing process that occurs at the convex point of a square solid obstacle must

TE D

satisfy the conservation of diffusively reflected particle number at this point. The unknown postcollision distribution functions of distribution functions of

X, Y

+

Z

\:

+

EP

+

and

,

[

Z,

,

+

]

=

[

and

X

+

Y

]

should be dependent on the incoming

+

\ . (30)

AC C

Also, these diffusively reflected distribution functions must follow the Maxwellian distribution as follows:

=

=

Z

Z

=

[

[

=

]

]

. (31)

By solving Eqs. (30) and (31), the unknown distribution functions can be computed as:

17

ACCEPTED MANUSCRIPT

=

[

=

]

=

Z

+

[

+

]

+

+

Z

+

[

+

]

+

+

+

[

+

]

+

+

+

[

+

]

+

+

Z Z [ Z ]

+

+

(

X

+

Y

+

\ ), (32)

(

X

+

Y

+

\ ), (33)

(

X

+

Y

+

\ ), (34)

(

X

+

Y

+

\ ), (35)

(

X

+

Y

+

\ ). (36)

RI PT

Z

+

SC

=

+

M AN U

=

Z

[

]

For a concave wall, the unknown distribution functions from solid surface toward fluid node ( Z )

write:

TE D

is bounced back at the concave point to conserve mass and momentum (Figure 6). Thus, we can

Z

=

\ . (37)

EP

Like Eq. (24) for bounce-back, the above presented diffusion reflection boundary conditions can be easily implemented to complex geometries like porous media due to its local property.

AC C

4.3.Periodic boundary condition The periodic boundary condition is considered to be the simplest one (in terms of implementation in the program) in the lattice Boltzmann method. In this condition, the distribution function exiting from one side of the domain is assumed to be equal from the opposite side. For example, for a two-dimensional flow field with an input on the west side and an outlet on the east, we have[44]: 18

ACCEPTED MANUSCRIPT ' (='( , C + BC) = ' (=wx* , C), T = 1,5,8 , (38) ' (=wx* , C

+ BC) = ' (='( , C), T = 3,6,7. (39)

RI PT

5. Model verification

In this section, the ability of the developed diffusion reflection boundary condition in LBM is evaluated to simulate the flow field in two different problems and the results are compared with

SC

those of the literature. In the first problem, the fluid flow characteristic in 2D micro-channel is simulated and the velocity profile at different Knudsen numbers is compared with that of the

M AN U

literature. In the second problem, two fractal porous media are selected, simulated, and their results are compared with those of the literature under the same conditions. In all cases, MRTLBM is used and the relation of Knudsen number and relaxation time is corrected by considering effective relaxation time (Eq. 21).

TE D

5.1.Micro-channel flow

Figure 7 illustrates a schematic of the simulated micro-channel geometry in this section. It

EP

should be noted that both the inlet and outlet boundaries are set as periodic boundary conditions, while for the upper and lower walls, the bounce-back and developed diffusion reflection boundary conditions are applied. Also, the flow field is driven by applying a constant body force

AC C

of 10 Z
in the x-axis. Moreover, the length and height of the channel are chosen 100Lu

and 20Lu (Lu stands for lattice unit), respectively. Fig. 8 shows the x-component velocity profile for the described micro-channel flow. The simulation is carried out for five different Knudsen numbers of 0.113, 0.226, 0.451, 0.677, and 1.13 and results are compared with those of Zhang et al.[31] with and without wall function

19

ACCEPTED MANUSCRIPT

(WF) and Ohwada et al.[45]. The x-component velocity is obtained by integration of velocity at

the middle of the channel (= = /2) and is normalized by the mean velocity at this section. As this figure indicates, the results of present study are in good agreement with those of the Zhang et

RI PT

al., especially when these groups applied wall function in their simulation to capture the effect of Knudsen layer more accurately. Also, the results of DSMC simulation are added to these figures to show how our results deviate at high Kn numbers. As shown in these figures, the results of

SC

Ohwada et al. and DSMC are in close agreement with each other, since the main concept behind these two methods are mainly similar. Moreover, to show the superiority of LBM versus Navier-

M AN U

Stokes based solution, the results of Navier-Stokes simulation presented in Zhang et al. is also presented. The results obviously indicate the demerits of N-S solution versus merits of LBM. To further study the use of the present model, variation of the normalized volume flow rate against Knudsen number is presented in Fig. 9 and the results are compared with available

TE D

numerical, theoretical, and experimental data. As figure 9 indicates, there is a minimum in the normalized flow rate in the micro-channel flow in the transition flow regime which is known as Knudsen phenomenon. Knudsen observed the minimum at about

~1. This behavior has been

EP

studied by many scholars both theoretically by Cercignani et al.[46] and Dongari et al.[47] and

AC C

experimentally by Dong [48]. The present model agrees very well with Dongari et al.[47] and Cercignani et al. [46] for Kn less than about 0.1, which corresponds approximately to the end of the slip flow regime. Also, our present model is relatively in good agreement with theoretical and experimental data at transition flow regime (0.1 <

< 10). It should be noted that this good

agreement is achieved without incorporating any kind of complex and wall adjustable slip models like wall function used in Zhang et al. [31]. 5.2.Fractal porous media 20

ACCEPTED MANUSCRIPT

In this section, two fractal geometries with porosities of 0.79 and 0.702 are selected and simulated at Darcy regime. With this regards, the permeability for each case at different Knudsen numbers of zero, 0.0369 and 0.0539 are calculated and compared with results of Zheng et al.[49]

RI PT

(Table 1). As mentioned earlier, the developed diffusion reflection boundary condition has a high capability to simulate the effect of rarefied gaseous flow in micro-porous media for both slip and transition flow regimes when the effective relaxation time is used, in which the results of Table 1

SC

indicate this fact.

M AN U

6. Results and discussions

In this study, pore-scale simulation of rarified gas flow in square- and circular-based Sierpinski carpets is carried out, using lattice Boltzmann method. For this purpose, the first three stages of these two-dimensional fractal geometries at slip and transition flow regimes are considered. The

TE D

effects of key parameters such as Knudsen number and porosity on the flow structures and permeability are investigated. At last, the relationship between the permeability with Knudsen number and porosity in the slip and transition flow regimes is investigated and expressed with a

EP

high coefficient of determination.

AC C

6.1.Results of square-based Sierpinski carpet To study the effects of slippage and porosity on the flow patterns, vortex structures, and streamlines along with the velocity contours (in the x-direction) for the square-based Sierpinski carpet in the first, second, and third stages, Figs. (10-12) are presented, respectively. Porosity for the first, second, and third stages are 0.888, 0.79, and 0.702, respectively, which can be calculated using Eq. (5). Simulations have been carried out for 10 [
21

constant body

ACCEPTED MANUSCRIPT

force applied along the longitudinal axis and in six different Knudsen numbers of 0, 0.07, 0.1, 0.3, 0.6, and 1.

RI PT

According to Fig. (10), one can say that the velocity value for the first stage of the square-based Sierpinski carpet is reduced as Knudsen number increases. The rate of velocity reduction in the slip flow regime is greater than that of the transition one. According to Eq. (20), as Knudsen number increases, the relaxation time increases, and hence the kinematic viscosity of the fluid

SC

increases. Increasing the Knudsen number increases slippage effect on the solid-fluid interface,

M AN U

which will increase the velocity at the interface[50]. However, augmentation of the kinematic viscosity shows a downward trend in the velocity, which is dominant on the rough surfaces. It is essentially due to the reduction of frictional force on the solid surfaces in high Knudsen numbers, which results in a reduction of the mass flow rate in regions far from the solid surface[50]. These findings are also consistent with the results of the other researchers[46, 47, 51]. As stated by

TE D

Cercignani et al.[46] in 2004, for solid surfaces that follow the Maxwellian distribution (such as the solid surface with diffusion properties), the mass flow rate decreases with the increase of the Knudsen number up to the commencement of the molecular dynamics regime (

~1). After this

EP

point, depending on the physical conditions of the problems, the mass flow rate can be increased

AC C

or remained constant. This reduction in the mass flow rate at low Knudsen numbers (slip flow regime) and its augmentation in the high Knudsen numbers (between transition and molecular dynamics flow regimes) cause Knudsen paradox[47]. This phenomenon which was first reported by Knudsen[52], occurs when the convective mass flow rate starts to be dominated by the diffusive mass flow rate [47]. Also, as Fig. 10 illustrates it can be said that the size of generated vortices behind the solid body due to the separation phenomenon is reduced as Knudsen number increases. So that, the vortex 22

ACCEPTED MANUSCRIPT

size is minimized at the Knudsen number of 0.3, after which the surface of the solid object acts as a streamline. This result is also consistent with the results observed by Tang et al.[50].

RI PT

In order to study more on the effect of Knudsen number and porosity on the flow pattern in fractal porous media, Fig. 11 is presented. This figure presents the velocity contours and streamlines at different Knudsen numbers in the second stage of the square-based Sierpinski carpet. Similar to the first stage (Fig. 10), the velocity decreases with the increase of the Knudsen

SC

number, so that the reduction rate of velocity in the slip flow regime is higher than that of the

M AN U

transition flow regime.

Moreover, according to Fig. 11 it can be said that the size of the produced vortices due to the separation phenomenon behind the solid object is decreased as Knudsen number increases (similar to the first stage). However, as predicted, the size of these vortices is smaller than that of

TE D

the vortices created in the first stage. Therefore, the size of the vortices is decreased as porosity decreases, where the size reduction rate at the second stage is much lower than that of the first stage one. So that, the presence of the small vortices in Kn=1 is observable yet.

EP

Finally, for further study on the effects of the fractal geometries, porosity, and Knudsen number on vortex structures, flow patterns, velocity contours and streamlines for the third stage of the

AC C

square-based Sierpinski carpet is illustrated in Fig. 12. One of the significant implications of this figure is that the size of the vortices at low porosities is considerably smaller than that of the high porosities, which is mainly due to the low flow momentum and high number of fractal geometries along the flow path. Also, according to Fig. 12, it can be said that the size of the produced vortices is very smaller than that of the previous stages that are not visible. Therefore, with the increase of the Knudsen

23

ACCEPTED MANUSCRIPT

number and decrease of the porosity, the size of the flow vortices is decreased and eventually disappeared due to the augmentation of the kinematic viscosity and fractal geometry of the porous medium. This will also increase the permeability. In this case (at high Knudsen numbers

RI PT

and low porosities), the flow regime is in the range of the Darcy regime, which no flow separation is observed and the effect of the Knudsen number on the flow patterns is negligible. Fig. 13 shows the variation of the permeability with the pressure gradient along the x-direction

SC

(body force) at different porosities and Knudsen numbers for the square-based Sierpinski carpet.

M AN U

For this purpose, the permeability is obtained for various body forces using Eq. (8). As this figure indicates, at high porosities, the permeability is decreased as body force increases. At high porosities, the reduction rate of the permeability strongly depends on the Knudsen number; so that in low rarefaction effects, the reduction rate of the permeability is very high. These results are consistent with the findings of Zheng et al.[49] and Tang et al.[50] for fractal micro-porous

TE D

media. On the other hand, in low porosities, the fluid flow behavior approaches to the Darcy regime and gets away from the Darcy-Forchimer regime. As a result, the permeability is almost constant with any variation of the body force in the x-direction in this regime. Another

EP

significant indication of Fig. 13 is that increasing the porosity will result in the augmentation of

AC C

the permeability which is not linear. This variation will further be discussed in Fig. 14. Fig. 14 shows the permeability variations with the porosity and Knudsen number for the squarebased Sierpinski carpet in the Darcy regime. According to this figure, it can be concluded that increasing porosity and Knudsen number increases permeability, as pointed out before. Throughout this study, it is found that the permeability varies logarithmically with Knudsen number, while it varies exponentially with the porosity. To show this fact accurately, the correlation between the permeability with the Knudsen number and porosity with a 24

ACCEPTED MANUSCRIPT

determination coefficient of 98.19% for the square-based Sierpinski carpet can be predicted in terms of Eq. (40). In this case, we have (0 <

< 1, zA ≪ 1, U S 0.702 ≤ - ≤ 0.888):

K = a0 × e b0ε + c 0 × Ln ( K n ) + d 0 × e b0ε × Ln ( Kn ) + e 0 ,

(40)

RI PT

where, UW , }W , IW , SW , and AW are constant coefficients and can be expressed with confident range of 95% as follows:

SC

a0 = 3 .853 × 10 −7 ( − 2 .022 × 10 −6 , 2 .793 × 10 −6 ) ,

c 0 = − 0 .2898 ( − 7 .633 , 7 .053) ,

M AN U

b 0 = 22 .84 (15 .85 , 29 .82) ,

(41)

d 0 = 3 .184 × 10 −9 ( − 2 .48 × 10 −8 , 3 .116 × 10 −8 ) ,

TE D

e 0 = 16 .25 (1 .364 , 31 .13).

6.2.Results of circular-based Sierpinski carpet Due to the ability of the developed program to simulate the gaseous rarefaction effects for

EP

different geometries (with first order accuracy), the circular-based Sierpinski carpet in slip and

AC C

transition flow regimes is considered, simulated and compared with the square-based Sierpinski carpet under the same conditions. Based on this fact, Fig. 15 shows the evolution of velocity contours, streamlines, and the size of the vortices for the first stage of the circular-based Sierpinski carpet (ε=0.912) with different Knudsen numbers. According to this figure and similar to Fig. 10, it can be said that the velocity magnitude decreases with the increase of the Knudsen number, where the reduction rate of velocity in the slip flow regime is considerably higher than

25

ACCEPTED MANUSCRIPT

that of the transition one. As mentioned earlier, this reduction rate is pronounced up to Kn=0.3, where it is not considerable thereafter.

RI PT

As another significant implication of Fig. 15, the size of the generated vortices due to the separation phenomenon behind the circular cylinder is decreased as Knudsen number increases (like the square shapes). However, the size of these vortices is reached to their minimum magnitude at lower Knudsen numbers (Kn=0.1), where the streamlines fit the body thereafter.

SC

This is mainly due to the geometry characteristics of the circular cylinder versus of the square

M AN U

one.

To study further the effects of Knudsen number and porosity on the circular-based Sierpinski carpet, the velocity magnitude and streamlines for the second (ε=0.835) and third (ε=0.766) stages of this fractal porous medium is presented in Figs. 16 and 17, respectively. As shown in

TE D

these figures, similar results are observed for these fractal geometries. To discuss further on quantitative variation of the velocity due to the rarefaction effects on both square- and circular-based Sierpinski carpets, Fig. 18 is presented. This figure is depicted for the third stage of the Sierpinski carpets at six different Kn numbers (

~0, 0.07, 0.1, 0.3, 0.6, and 1)

~0. As shown in figures 18(a) and 18(b), the quantitative reduction of velocity at the

AC C

1 and

EP

and ~/ = 1 . The velocity profiles are normalized with the mean velocity integrated at ~/ =

slip flow regime is considerably higher than that of the transition flow regime, particularly at peak points (at maximum velocities). To show detail variation of the velocity profiles at transition flow regime, the figures are magnified at the peak point, and presented in a clear frame.

26

ACCEPTED MANUSCRIPT

Similar to Fig. 13 for the square-based Sierpinski carpet, Fig. 19 is plotted to show the variation of the permeability with the pressure gradient along the x-direction (body force) at different porosities and Knudsen numbers for the circular-based Sierpinski carpet. Likewise, the

RI PT

permeability is decreased at high porosities as body force increases. At high porosities, the reduction rate of the permeability strongly depends on the Knudsen number. Therefore, the reduction rate of the permeability is very high at low rarefaction effects. These results are in

SC

close agreement with those of the square-based Sierpinski carpet, too.

M AN U

Fig. 20 shows the variation of the permeability with porosity and Knudsen number for the circular-based Sierpinski carpet in the Darcy regime. According to this figure, it can be concluded that like the square-based Sierpinski carpet, the permeability varies logarithmically with Knudsen number and exponentially with the porosity. With this regard, the correlation between the permeability with the Knudsen number and porosity (with a determination

TE D

coefficient of 96.54%) can be predicted as follows (0 < 0.912):

K = a00 × e b 00ε + c 00 × Ln ( K n ) + d 00 × e b 00ε × Ln ( K n ) + e 00 ,

< 1, zA ≪ 1, U S 0.766 ≤ - ≤

(42)

EP

where, UWW , }WW , IWW , SWW , and AWW are constant coefficients and can be expressed with confident

AC C

range of 95% as follows:

a00 = 4 .946 × 10 −9 ( − 4 .745 × 10 −8 , 5 .734 × 10 −8 ) , b 00 = 27 .21(15 .67 , 38 .75) ,

c 00 = 1 .416 ( − 11 .29 ,14 .12) ,

(43)

d 00 = 1 .579 × 10 −10 ( − 1 .554 × 10 −9 , 1 .87 × 10 −9 ) ,

27

ACCEPTED MANUSCRIPT

e 00 = 28 .15 (1 .229 , 55 .07).

7. Concluding remarks

RI PT

Pore-scale simulation of the rarefied gas flow inside the square- and circular-based Sierpinski carpets at different porosities and Knudsen numbers are conducted, using lattice Boltzmann method. For this purpose, a simple methodology for diffusion reflection boundary condition is

SC

developed and validated for two valuable benchmarks, namely micro-channel flow and fractal porous media, where the results were in good agreement with literature. In our main problem

M AN U

(Sierpinski carpets), the effects of different key parameters such as pressure gradient in the xdirection (body force), porosity and Knudsen number on the velocity contours and streamlines, velocity profiles, flow structures, and permeability are investigated. Also, a correlation between the permeability with porosity and Knudsen number in the Darcy regime ( Re << 1 ) is presented.

•

TE D

Some concluding remarks are as follows:

Velocity magnitude is reduced as Knudsen number increases, where the reduction rate of velocity in the slip flow regime is greater than that of the transition one. Permeability varies logarithmically with Knudsen number and exponentially with the porosity.

EP

•

At high porosities, the permeability is decreased as body force increases.

•

At high porosities and low rarefaction effects, the reduction rate of the permeability with

AC C

•

body force is considerable (Darcy-Forchimer regime). •

At low porosities, the fluid flow behavior approaches that of the Darcy regime and gets away from the Darcy-Forchimer regime. As a result, the permeability is almost constant with any variation of the body force in the x-direction in the Darcy regime. 28

ACCEPTED MANUSCRIPT

•

Increasing the porosity and Knudsen number will increase the permeability in the Darcy regime.

[5] [6] [7] [8]

[9] [10]

[11]

[12]

[13] [14] [15]

SC

M AN U

[4]

TE D

[3]

EP

[2]

Q. Zheng, J. Xu, B. Yang, and B. Yu, "Research on the effective gas diffusion coefficient in dry porous media embedded with a fractal-like tree network," Physica A: Statistical Mechanics and its Applications, vol. 392, pp. 1557-1566, 2013. A. Hosa, A. Curtis, and R. Wood, "Calibrating Lattice Boltzmann flow simulations and estimating uncertainty in the permeability of complex porous media," Advances in Water Resources, vol. 94, pp. 60-74, 2016. B. B. Mandelbrot, "Geometry of Turbulence: Intermittency," The fractal geometry of nature, vol. 173, pp. 97-105, 1983. G. R. Rakhshandehroo, M. Shaghaghian, A. Keshavarzi, and N. Talebbeydokhti, "Temporal variation of velocity components in a turbulent open channel flow: Identification of fractal dimensions," Applied Mathematical Modelling, vol. 33, pp. 38153824, 2009. W. Cao, Z. Cai, and Z. Tang, "Fractal structure of lunar topography: An interpretation of topographic characteristics," Geomorphology, vol. 238, pp. 112-118, 2015. G. Lei, P. Dong, S. Mo, S. Gai, and Z. Wu, "a Novel Fractal Model for Two-Phase Relative Permeability in Porous Media," Fractals, vol. 23, p. 1550017, 2015. C. G. Jacquin and P. Adler, "Fractal porous media II: geometry of porous geological structures," Transport in Porous Media, vol. 2, pp. 571-596, 1987. D. Aniszewska, "Fractal characteristics of defects growth in porous ceramics modeled with the Movable Cellular Automata," Applied Mathematical Modelling, vol. 39, pp. 2409-2415, 2015. A. J. Katz and A. Thompson, "Fractal sandstone pores: implications for conductivity and pore formation," Physical Review Letters, vol. 54, p. 1325, 1985. B. Yu, M. Zou, and Y. Feng, "Permeability of fractal porous media by Monte Carlo simulations," International journal of heat and mass transfer, vol. 48, pp. 2787-2794, 2005. R. Pitchumani and B. Ramakrishnan, "A fractal geometry model for evaluating permeabilities of porous preforms used in liquid composite molding," International journal of heat and mass transfer, vol. 42, pp. 2219-2232, 1999. B. Yu, "Comments on “A fractal geometry model for evaluating permeabilities of porous preforms used in liquid composite molding”," International journal of heat and mass transfer, vol. 44, pp. 2787-2789, 2001. B. Yu, L. J. Lee, and H. Cao, "A fractal in‐plane permeability model for fabrics," Polymer Composites, vol. 23, pp. 201-221, 2002. B. Yu, J. Li, and D. Zhang, "A fractal trans-plane permeability model for textile fabrics," International communications in heat and mass transfer, vol. 30, pp. 127-138, 2003. B. Yu and W. Liu, "Fractal analysis of permeabilities for porous media," AIChE journal, vol. 50, pp. 46-57, 2004.

AC C

[1]

RI PT

References

29

ACCEPTED MANUSCRIPT

[22]

[23] [24] [25] [26] [27]

[28]

[29]

[30]

[31]

[32]

RI PT

[21]

SC

[20]

M AN U

[19]

TE D

[18]

EP

[17]

F. Meng, H. Zhang, Y. Li, X. Zhang, and F. Yang, "Application of fractal permeation model to investigate membrane fouling in membrane bioreactor," Journal of Membrane Science, vol. 262, pp. 107-116, 2005. M. Othman and Z. Helwani, "Simulated fractal permeability for porous membranes," Applied Mathematical Modelling, vol. 34, pp. 2452-2464, 2010. M. Gad-el-Hak, "The fluid mechanics of microdevices-the Freeman scholar lecture," Transactions-American Society of Mechanical Engineers Journal of FLUIDS Engineering, vol. 121, pp. 5-33, 1999. R. N. Moghaddam and M. Jamiolahmady, "Slip flow in porous media," Fuel, vol. 173, pp. 298-310, 2016. W.-M. Zhang, G. Meng, and X. Wei, "A review on slip models for gas microflows," Microfluidics and nanofluidics, vol. 13, pp. 845-882, 2012. Y. Zhang, R. Qin, Y. Sun, R. W. Barber, and D. Emerson, "Gas flow in microchannels–a lattice Boltzmann method approach," Journal of Statistical Physics, vol. 121, pp. 257267, 2005. C. L. Bailey, R. W. Barber, and D. R. Emerson, "Is it safe to use Navier–Stokes for gas microflows," in European Congress on Computational Methods in Applied Sciences and Engineering, 2004. X. Nie, G. D. Doolen, and S. Chen, "Lattice-Boltzmann simulations of fluid flows in MEMS," Journal of statistical physics, vol. 107, pp. 279-289, 2002. S. Succi, "Mesoscopic modeling of slip motion at fluid-solid interfaces with heterogeneous catalysis," Physical review letters, vol. 89, p. 064502, 2002. C. Lim, C. Shu, X. Niu, and Y. Chew, "Application of lattice Boltzmann method to simulate microchannel flows," Physics of fluids, vol. 14, pp. 2299-2308, 2002. M. Sbragaglia and S. Succi, "Analytical calculation of slip flow in lattice Boltzmann models with kinetic boundary conditions," Physics of Fluids, vol. 17, p. 093602, 2005. P. Taheri, M. Torrilhon, and H. Struchtrup, "Couette and Poiseuille microflows: analytical solutions for regularized 13-moment equations," Physics of Fluids, vol. 21, p. 017102, 2009. I.-W. Park, M.-S. Shin, S.-J. Byun, and J.-Y. Yoon, "Simulation of gas flow in a microchannel by lattice Boltzmann method," in Fluid Machinery and Fluid Mechanics, ed: Springer, 2009, pp. 195-200. Z. Guo, T. Zhao, and Y. Shi, "Physical symmetry, spatial accuracy, and relaxation time of the lattice Boltzmann equation for microgas flows," Journal of Applied physics, vol. 99, p. 074903, 2006. F. Verhaeghe, L.-S. Luo, and B. Blanpain, "Lattice Boltzmann modeling of microchannel flow in slip flow regime," Journal of Computational Physics, vol. 228, pp. 147-157, 2009. Y.-H. Zhang, X.-J. Gu, R. W. Barber, and D. R. Emerson, "Capturing Knudsen layer phenomena using a lattice Boltzmann model," Physical Review E, vol. 74, p. 046704, 2006. A. Isfahani and A. Soleimani, "Thermal Lattice Boltzmann Simulation of Rarefied Gas Flows in Nanochannels for Wide Range of Knudsen Number," in Advanced Materials Research, 2012, pp. 5313-5317.

AC C

[16]

30

ACCEPTED MANUSCRIPT

[38]

[39] [40]

[41]

[42] [43]

[44] [45]

[46] [47]

[48] [49]

[50]

RI PT

SC

[37]

M AN U

[36]

TE D

[35]

EP

[34]

Z.-L. Chen, N.-T. Wang, L. Sun, X.-H. Tan, and S. Deng, "Prediction method for permeability of porous media with tortuosity effect based on an intermingled fractal units model," International Journal of Engineering Science, vol. 121, pp. 83-90, 2017. B. Yu and J. Li, "Some fractal characters of porous media," Fractals, vol. 9, pp. 365-372, 2001. A. E. Khabbazi, J. Hinebaugh, and A. Bazylak, "Analytical tortuosity–porosity correlations for Sierpinski carpet fractal geometries," Chaos, Solitons & Fractals, vol. 78, pp. 124-133, 2015. M. Kaviany, "Fluid Mechanics," in Principles of Heat Transfer in Porous Media, ed: Springer, 1991, pp. 15-113. M. Salimi, M. Taeibi-Rahni, and F. Jam, "Heat transfer analysis of a porously covered heated square cylinder, using a hybrid Navier–Stokes–lattice Boltzmann numerical method," International Journal of Thermal Sciences, vol. 91, pp. 59-75, 2015. Y. Yuan and S. Rahman, "Extended application of lattice Boltzmann method to rarefied gas flow in micro-channels," Physica A: Statistical Mechanics and its Applications, vol. 463, pp. 25-36, 2016. Z. Guo and C. Shu, Lattice Boltzmann method and its applications in engineering: World Scientific, 2013. V. Sofonea and R. F. Sekerka, "Boundary conditions for the upwind finite difference Lattice Boltzmann model: Evidence of slip velocity in micro-channel flow," Journal of Computational Physics, vol. 207, pp. 639-659, 2005. W. G. Vincenti and C. H. Kruger, "Introduction to physical gas dynamics," Introduction to physical gas dynamics, by Vincenti, Walter Guido; Kruger, Charles H. New York, Wiley [1965], 1965. C. Cercignani, Rarefied gas dynamics: from basic concepts to actual calculations vol. 21: Cambridge University Press, 2000. Y. Sone, "Theoretical and Numerical Analyses of the Boltzmann Equation-Theory and Analysis of Rarefied Gas Flows-Part I," Lecture Notes, Department of Aeronautics and Astronautics, Graduate School of Engineering, Kyoto University, Kyoto, 1998. A. A. Mohamad, Lattice Boltzmann method: fundamentals and engineering applications with computer codes: Springer Science & Business Media, 2011. T. Ohwada, Y. Sone, and K. Aoki, "Numerical analysis of the Poiseuille and thermal transpiration flows between two parallel plates on the basis of the Boltzmann equation for hard‐sphere molecules," Physics of Fluids A: Fluid Dynamics, vol. 1, pp. 2042-2049, 1989. C. Cercignani, M. Lampis, and S. Lorenzani, "Variational approach to gas flows in microchannels," Physics of Fluids, vol. 16, pp. 3426-3437, 2004. N. Dongari, A. Sharma, and F. Durst, "Pressure-driven diffusive gas flows in microchannels: from the Knudsen to the continuum regimes," Microfluidics and nanofluidics, vol. 6, pp. 679-692, 2009. W. Dong, "University of California report no," UCRL-3353, 1956. J. Zheng, W. Zhang, G. Zhang, Y. Yu, and S. Wang, "Effect of porous structure on rarefied gas flow in porous medium constructed by fractal geometry," Journal of Natural Gas Science and Engineering, vol. 34, pp. 1446-1452, 2016. G. Tang, W. Tao, and Y. He, "Gas slippage effect on microscale porous flow using the lattice Boltzmann method," Physical Review E, vol. 72, p. 056301, 2005.

AC C

[33]

31

ACCEPTED MANUSCRIPT

RI PT SC M AN U TE D EP

[52]

Z. Deng, Y. Chen, and C. Shao, "Gas flow through rough microchannels in the transition flow regime," Physical Review E, vol. 93, p. 013128, 2016. M. Knudsen, "Die Gesetze der Molekularströmung und der inneren Reibungsströmung der Gase durch Röhren," Annalen der Physik, vol. 333, pp. 75-130, 1909.

AC C

[51]

32

ACCEPTED MANUSCRIPT

Tables:

ε = 0.79

Kn

RI PT

Table 1 Comparison of the permeability value in Darcy regime between present study and Zheng et al. for different Knudsen numbers and porosities.

ε = 0.702

Present work

Zheng et al.

Present work

Kn → 0

4.4 × 10 − 15

4.6 × 10 − 15

3.91 × 10 − 16

3.84 × 10 − 16

Kn = 0.0369

9.36 × 10 − 15

9.94 × 10 − 15

3.14 × 10 − 15

3.41 × 10 − 15

Kn = 0.0539

12.86 × 10 − 15

13.04 × 10 − 15

5.34 × 10 − 15

5.43 × 10 − 15

AC C

EP

TE D

M AN U

SC

Zheng et al.

ACCEPTED MANUSCRIPT

Figure 1 Generation process of different stages of the square-based Sierpinski carpet: (a) first stage, (b) second stage, and (c) third stage.

RI PT

Figure 2 Generation process of different stages of the circular-based Sierpinski carpet: (a) first stage, (b) second stage, and (c) third stage.

Figure 3 Boundary conditions for the distribution functions: (a) bounce-back and (b) diffusion

SC

reflection. Figure 4 Boundary nodes at the bottom wall for the D2Q9 model.

M AN U

Figure 5 Boundary nodes at the convex wall of the smallest solid obstacle for the D2Q9 model. Figure 6 Boundary nodes at the concave wall of the smallest solid obstacle for the D2Q9 model. Figure 7 Schematic diagram of a micro-channel geometry.

TE D

Figure 8 Comparison of the x-component velocity profile with the results of other studies for the flow inside the micro-channel at different Knudsen numbers. Figure 9 Variation of normalized volume flow rate against Knudsen number and comparison

EP

with available numerical and experimental data.

AC C

Figure 10 X-component velocity contours and streamlines around the square-based Sierpinski carpet for the first stage at various Knudsen numbers. Figure 11 X-component velocity contours and streamlines around the square-based Sierpinski carpet for the second stage at various Knudsen numbers. Figure 12 X-component velocity contours and streamlines around the square-based Sierpinski carpet for the third stage at various Knudsen numbers.

ACCEPTED MANUSCRIPT

Figure 13 Variation of the permeability versus body force at various porosities and Knudsen numbers for the square-based Sierpinski carpet.

based Sierpinski carpet.

RI PT

Figure 14 Variation of the permeability versus porosity and Knudsen number for the square-

carpet at the first stage for various Knudsen numbers.

SC

Figure 15 X-component velocity contours and streamlines around the circular-based Sierpinski

Figure 16 X-component velocity contours and streamlines around the circular-based Sierpinski

M AN U

carpet at the second stage for various Knudsen numbers.

Figure 17 X-component velocity contours and streamlines around the circular-based Sierpinski carpet at the third stage for various Knudsen numbers.

Figure 18 X-component velocity profile at X/L=1 for square- and circular-based Sierpinski

TE D

carpet at 3rd stage and different Knudsen numbers. Figure 19 Variation of the permeability versus body force at various porosities and Knudsen

EP

numbers for the circular-based Sierpinski carpet. Figure 20 Variation of the permeability versus porosity and Knudsen number for the circular-

AC C

based Sierpinski carpet.

(b)

(c)

AC C

(a)

EP

TE D

M AN U

Figure 1

SC

(a)

RI PT

ACCEPTED MANUSCRIPT

(b) Figure 2

(c)

(a)

SC

RI PT

ACCEPTED MANUSCRIPT

(b)

AC C

EP

TE D

M AN U

Figure 3

Figure 4

TE D

M AN U

SC

RI PT

ACCEPTED MANUSCRIPT

AC C

EP

Figure 5

TE D

M AN U

SC

RI PT

ACCEPTED MANUSCRIPT

AC C

EP

Figure 6

M AN U

SC

RI PT

ACCEPTED MANUSCRIPT

AC C

EP

TE D

Figure 7

(a)

(b)

AC C

EP

TE D

(c)

M AN U

SC

RI PT

ACCEPTED MANUSCRIPT

(e) Figure 8

(d)

M AN U

SC

RI PT

ACCEPTED MANUSCRIPT

AC C

EP

TE D

Figure 9

(a)

(b)

SC

RI PT

ACCEPTED MANUSCRIPT

(d)

AC C

EP

(e)

TE D

M AN U

(c)

Figure 10

(f)

SC

RI PT

ACCEPTED MANUSCRIPT

(b)

(d)

AC C

EP

(c)

TE D

M AN U

(a)

(e)

(f) Figure 11

SC

RI PT

ACCEPTED MANUSCRIPT

(b)

(d)

AC C

EP

(c)

TE D

M AN U

(a)

(e)

(f) Figure 12

AC C EP

Permeability (Lu2)

TE D M AN U SC

RI PT

ACCEPTED MANUSCRIPT

Figure 13

(a)

ACCEPTED MANUSCRIPT

First Stage Square-based Sierpinski Carpet 270

K= 4.1773ln(Kn) + 268.15 R² = 0.9868

RI PT

266 264 262 260

SC

Permeability, K (Lu2)

268

258 256 0.2

0.4 0.6 0.8 Knudsen Number, Kn

1

M AN U

0

(b)

Second Stage Square-based Sierpinski Carpet

TE D

24

K = 1.5287ln(Kn) + 22.083 R² = 0.9738

22

EP

21 20 19

AC C

Permeability, K (Lu2)

23

18 17

0

0.2

0.4

0.6

0.8

Knudsen Number, Kn

(c)

1

1.2

1.2

ACCEPTED MANUSCRIPT

Third Stage Square-based Sierpinski Carpet K= 0.566ln(Kn) + 3.2174 R² = 0.9847

RI PT

3

2.5

2

1.5 0

0.2

0.4

0.6

SC

Permeability, K (Lu2)

3.5

0.8

1

M AN U

Knudsen Number, Kn

(d)

AC C

EP

TE D

Figure 14

(a)

(b)

1.2

SC

RI PT

ACCEPTED MANUSCRIPT

(d)

(e)

TE D

M AN U

(c)

AC C

EP

Figure 15

(f)

ACCEPTED MANUSCRIPT

(b)

SC

RI PT

(a)

(d)

AC C

EP

(e)

TE D

M AN U

(c)

Figure 16

(f)

SC

RI PT

ACCEPTED MANUSCRIPT

(b)

(d)

AC C

EP

(c)

TE D

M AN U

(a)

(e)

(f) Figure 17

M AN U

SC

RI PT

ACCEPTED MANUSCRIPT

AC C

EP

TE D

(a) Square-based Sierpinski carpet at 3rd stage

(a) Circular-based Sierpinski carpet at 3rd stage Figure 18

M AN U

SC

RI PT

ACCEPTED MANUSCRIPT

AC C

EP

TE D

Figure 19

(a)

ACCEPTED MANUSCRIPT

First Stage Circular-based Sierpinski Carpet

330

RI PT

Permeability, K (Lu2)

340

K = 12.538ln(Kn) + 337.11 R² = 0.9819

320 310

290 0

0.2

0.4

0.6

SC

300

0.8

1.2

M AN U

Knudsen Number, Kn

1

(b)

Second Stage Circular-based Sierpinski Carpet 30

TE D

K= 3.3103ln(Kn) + 27.989 R² = 0.9694

26 24 22

AC C

20

EP

Permeability, K (Lu2)

28

18

0

0.2

0.4

0.6

Knudsen Number, Kn

(c)

0.8

1

1.2

ACCEPTED MANUSCRIPT

Third Stage Circular-based Sierpinski Carpet 5.5

K = 1.0342ln(Kn) + 4.8706 R² = 0.9845

4.5

RI PT

Permeability, K(Lu2)

5

4 3.5 3

SC

2.5 2 0.2

0.4 0.6 0.8 Knudsen Number, Kn

M AN U

0

(d)

AC C

EP

TE D

Figure 20

1

1.2

ACCEPTED MANUSCRIPT

Highlights •

A novel methodology in implementation of diffusion reflection boundary condition for

•

RI PT

rarefied gaseous flow is developed. Proposed diffusion reflection boundary condition is validated by two valuable benchmarks.

Gaseous rarefaction effects for two Sierpinski carpets with fractal geometries are

SC

•

investigated.

M AN U

A general correlation between permeability with porosity and Knudsen number for both

EP

TE D

Sierpinski carpets at Darcy regime is presented.

AC C

•