# Thermal Lattice Boltzmann Model for Nonisothermal Gas Flow in a Two-Dimensional Microchannel

Authors: Youssef Elguennouni; Mohamed Hssikou; Jamal Baliti; Mohammed Alaoui
Journal: Mathematical Problems in Engineering
 (2020)

---

## Abstract

In this paper, the Thermal Lattice Boltzmann Method (TLBM) is used for the simulation of a gas microflow. A 2D heated microchannel flow driven by a constant inlet velocity profileUin and nonisothermal walls is investigated numerically. Two cases of micro-Poiseuille flow are considered in the present study. In the first case, the temperature of the walls is kept uniform, equal to zero; therefore, the gas is driven along the channel under the inlet parameters of velocity and temperature. However, in the second one, the gas flow is also induced by the effect of temperature decreasing applied on the walls. For consistent results, velocity slip and temperature jump boundary conditions are used to capture the nonequilibrium effects near the walls. The rarefaction effects described by the Knudsen number, on the velocity and temperature profiles are evaluated. The aim of this study is to prove the efficiency of the TLBM method to simulate Poiseuille flow in case of nonisothermal walls, based on the average value of the Nusselt number and by comparing the results obtained from the TLBM with those obtained using the Finite Difference Method (FDM). The results also show an interesting sensitivity of velocity and temperature profiles with the rarefaction degree and the imposed temperature gradient of the walls.


---

## Body

## 1. Introduction

Microdevice technology has shown an interesting growth in the last few decades. In order to understand the behavior and the physics of gas flows in such micro-electro-mechanical systems (MEMS) better, researchers are focused on several approaches. Kinetically, such flows are governed by the Boltzmann equation, in which the solution is better approximated by direct simulation of Monte Carlo (DSMC) [1] and dynamic molecular (MD) [2]. To save the computation time, other alternatives have been used mainly in the slip regime such as moment equations [3]. To combine the advantages of both approaches, the lattice Boltzmann method (LBM) becomes recently a powerful tool of such applications, and this approach is a hybrid method which combines the kinetic description given by the Boltzmann equation and the classical computational fluid dynamics (CFD), mesoscopic approach. Several attempts are made by scientists to improve the LBM approach and extend its ability to simulate more complex geometry flows [4–6]. This method is used to simulate different types of flows: the Rayleigh-Benard convection [7], micro-Poiseuille flow [8–12], micro-Couette flow [11, 12], Lid-driven cavity [4, 12], etc.However, the study of such flow needs a good choice and implementation of the boundary conditions (BC), which is a crucial step in the LBM simulation. In this context, different BC are tested in the literature for different problems. Nie et al. [13] used the bounce-back boundary conditions, which is compared with the DSMC method [14]. Lim et al. [15] employed specular reflection and a second-order extrapolation scheme to capture the slip effect on the interaction gas surface. To test a hybrid BC, Tang et al. [16] defined a reflection coefficient rf to combine the bounce-back and the specular reflection boundary conditions [12] and also proposed a thermal boundary condition for a double-population thermal lattice Boltzmann method equation [17]. Lee and Lin [18] proposed a wall equilibrium boundary condition with a Knudsen number to capture the slip velocity. A slip boundary condition based on the Maxwellian-scattering kernel was proposed by Zhang et al. [19, 20] using a tangential momentum accommodation coefficient (TMAC). To simulate isothermal two-dimensional microchannel flows Niu et al. [21, 22] have used the diffuse scattering boundary condition [12]. Zou and He [23] proposed a new method to specify boundary conditions, where the conditions are constructed in the consistency with wall boundary condition, based on the idea of bounce-back of the nonequilibrium distribution function. Another investigation was carried out by Tian et al. [24, 25] who used the first order of Maxwell slip boundary conditions.The LBM approximation is a linear discretization of the Boltzmann equation. The collision term in this equation is, often, approximated according to the Bhatnagar, Gross, and Krook (BGK) [26]. Gas flows in microfluidic devices are characterized by their rarefaction degree given by the Knudsen number Kn=λ/l (λ is the mean free path of the gas molecules and l is the characteristic length of the system). According to its value, four regimes of the gas flow are considered: the continuum regime Kn≲0.001, the slip regime 0.001≲Kn≲0.1, the transition regime 0.1≲Kn≲10, and the free molecular regime Kn≳10. As the Knudsen number increases, rarefaction effects become more important and the continuum approach breaks down because of nonequilibrium effects. The miniaturization of devices, MEMS technology, leads to the microsystems where the gas particles mean free path is comparable to the system characteristic length. Therefore, accurate approaches are needed to better understand the physics of such microflows, which usually belong to slip or early transition regimes. In this context, with the computer performances increasing, several numerical methods are adopted as alternatives to study such flows whose analytical description is generally more complicated. The Lattice Boltzmann method is one of the most successful numerical approaches used for gas flow simulation, unlike others which need a long time of running (DSMC, for example). In this model, the governing equations of mass and momentum conservation are satisfied at each lattice nodes.The thermal lattice Boltzmann method is used to simulate two- and three-dimensional microchannel flows by using velocity slip and temperature jump boundary conditions [24, 25] and has been studied in several papers [8–11, 27–30]. In this paper, the TLBM and FDM are used to simulate micro-Poiseuille flow, in the slip regime, which is usually encountered in the MEMS devices [31]. According to the top and bottom walls temperature, two cases are considered; in the first case, duct walls are assumed isothermal; however, nonisothermal walls are imposed in the second case. The aim of this study is to prove the efficiency of the TLBM method to simulate micro-Poiseuille flow even in case of nonisothermal walls by comparing velocity and temperature profiles with the finite difference method.


## 2. Statement of Problem

In the present simulation, we studied forced convection inside a long microchannel in the case of the body force is zero. The microchannel walls are stationary and heated such that the temperature gradient is constant. A rectangular microchannel is considered with a large cross-sectional aspect ratio, and the flow is independent of directionz [32]. The flow is fully developed between two long and parallel plates; the characteristic length scale in this problem is the channel height H (Figure 1). According to the top and bottom walls temperature, two cases are considered in this study. In the first case, the duct walls are isothermal; therefore, the inlet parameters of temperature and velocity induce the gas. However, nonisothermal walls are imposed in the second case; the bottom and upper sides are linearly heated from the hot value TpTp=TH/p to cold one TC, so Tp varies from TH to TC (p is an integer ≥1). The variation of walls temperature in the second case is as follows:(1a)Tp=THp,p≥1,(1b)Tx=Tp−Tp−TC×xL.Figure 1
Studied domain configuration. (a) Case 1: walls with uniform temperatureTC and (b) case 2: walls with a linear decreasing temperature from Tp to TC.
(a)
(b)


## 3. Thermal LB Model

### 3.1. Lattice Boltzmann Method

A square grid and D2Q9 model is used for both distribution function densityf and temperature g. The governing equations for these distribution functions are written according to the BGK model as follows:(2a)fkx+ckΔt,t+Δt=fkx,t1−1τf+1τffkeqx,t,(2b)gkx+ckΔt,t+Δt=gkx,t1−1τg+1τggkeqx,t.In which,τf and τg represent the relaxation times of density and internal energy functions, respectively, and are related to the kinematic viscosity v and thermal diffusivity α by(3)τf=νcs2+0.5,τg=αcs2+0.5,where cs is the velocity of sound, related to the lattice speed, temperature, and gas constant by cs=c/3=RT.The relationship between the kinematic viscosityv and the coefficient of thermal diffusion α is determined by ν/α=Pr, where Pr is the Prandtl number.At the equilibrium state, the Maxwell distribution functionsfkeq and gkeq are written in the Taylor expansion as follows:(4a)fkeqi,j=wkρi,j1+3ck⋅uc2+92ck⋅u2c4−32u2c2,(4b)gkeqi,j=wkTi,j1+3ck⋅uc2+92ck⋅u2c4−32u2c2,where the weight factors wk are w0=4/9, w1−4=1/9, and w5−8=1/36 and c is the lattice speed given by c=Δx/Δt=Δy/Δt.The discrete scheme D2Q9 is shown in Figure2(5)ck=0,0,k=0,ccosk−1π2,sink−1π2,k=1−4,2ccosk−5π2+π4,sink−5π2+π4,k=5−8.Figure 2
D2Q9 velocities.The macroscopic densityρ, velocity u, and internal energy by unit of mass e are obtained by the following relations:(6a)ρ=∑k=08fk,(6b)ρu=∑k=08fkck,(6c)ρe=∑k=08gk,where e=DRT/2, in which D is the number of physical dimensions (equal to 2 in the current work).


### 3.2. Boundary Conditions

#### 3.2.1. Flow Boundary Conditions

The boundary conditions proposed by Zou and He [23] are used at the inlet, in which the normal velocity component is assumed to be equal to zero and the density is to be determined. After the streaming step, the unknown distribution functions are at the inlet boundary f1,f5,f8. Using equations (6a) and (6b), these function densities are calculated as follows:(7a)ρin=f0+f2+f4+2f3+f6+f71−Uin,(7b)f1=f3+23ρinUin,(7c)f5=f7−12f2−f4+16ρinUin,(7d)f8=f6+12f2−f4+16ρinUin.A simple extrapolation scheme is used for the velocity at the outlet.To capture velocity slip phenomenon, the Maxwell first order of slip boundary condition is applied in its dimensionless form, near the longitudinal walls [24]:(8a)ux,y=0slip=ux,y=0−ux,wall=σKn∂u∂yy=0,(8b)ux,y=Hslip=ux,wall−ux,y=H=σKn∂u∂yy=H,where ux,y=0slip and ux,y=Hslip are the slip velocities at the bottom and top walls. The walls are considered at rest ux,wall=0 and σ is the momentum accommodation coefficient which is assumed to be equal unity, to simulate completely diffuse reflection [33].At the top wall, the unknown distribution functionsf4,f7,f8 are calculated as follows:(9a)ux,y=H=λ4ux,H−1−ux,H−22+3λ,(9b)ρtop=f0+f1+f3+2f2+f6+f5,(9c)f4=f2,(9d)f7=ρtop1−ux,y=H−f0+f2+f42−f3+f6,(9e)f8=ρtop1+ux,y=H−f0+f2+f42−f1+f5,where the gas mean free path is defined as λ=KnH and H is the total number of lattice nodes in the vertical direction.


#### 3.2.2. Temperature Boundary Conditions

At the inlet boundary, the unknown energy distribution functionsg1,g5,g8 are obtained as follows:(10a)g1=Tinletw1+w3−g3,(10b)g5=Tinletw5+w7−g7,(10c)g8=Tinletw6+w8−g6.The outlet boundary condition for energy distributions is implemented by using an extrapolation scheme [34]. Similarly, the temperature jump boundary conditions caused by the rarefaction effects at the walls are given by(11a)Tx,y=0jump=Tx,y=0−Tx,wall=ϕKn2γγ+1Pr∂T∂yy=0,(11b)Tx,y=Hjump=Tx,wall−Tx,y=H=ϕKn2γγ+1Pr∂T∂yy=H.C
jump is the temperature jump coefficient defined as follows:(12)Cjump=ϕ2γγ+1Prλ=kλ.In which,ϕ represents the thermal accommodation coefficient and is assumed to be unity, γ is the specific heat ratio, and Pr is the Prandtl number.Temperature jump boundary conditions are written at the bottom and top walls as follows:(13a)Tx,y=0=Cjump4Tx,1−Tx,2+2Tx,wall2+3Cjump,(13b)Tx,y=H=Cjump4Tx,H−1−Tx,H−2+2Tx,wall2+3Cjump.At the top wall, the unknown distribution functionsg4,g7,g8 are calculated as follows:(14a)Tx,y=H=Cjump4Tx,H−1−Tx,H−2+2Tx,wall2+3Cjump,(14b)g4=Tx,y=Hw2+w4−g2,(14c)g7=Tx,y=Hw5+w7−g5,(14d)g8=Tx,y=Hw6+w8−g6.


### 3.3. Nusselt Number Calculation

Heat transfer characteristic of the flow can be determined using the Nusselt numberNu, which is the ratio of convective and conductive heat transfers. The Nusselt number along the horizontal axis is calculated as follows:(15)Nux=2H∂T/∂ywallTwall−Tbulk.The bulk temperature at a cross-section is calculated as follows:(16)Tbulk=∫SρuTdS∫SρudS.The average Nusselt number is calculated as follows:(17)Nuavg=∫0LNuxdxL.


## 3.1. Lattice Boltzmann Method

A square grid and D2Q9 model is used for both distribution function densityf and temperature g. The governing equations for these distribution functions are written according to the BGK model as follows:(2a)fkx+ckΔt,t+Δt=fkx,t1−1τf+1τffkeqx,t,(2b)gkx+ckΔt,t+Δt=gkx,t1−1τg+1τggkeqx,t.In which,τf and τg represent the relaxation times of density and internal energy functions, respectively, and are related to the kinematic viscosity v and thermal diffusivity α by(3)τf=νcs2+0.5,τg=αcs2+0.5,where cs is the velocity of sound, related to the lattice speed, temperature, and gas constant by cs=c/3=RT.The relationship between the kinematic viscosityv and the coefficient of thermal diffusion α is determined by ν/α=Pr, where Pr is the Prandtl number.At the equilibrium state, the Maxwell distribution functionsfkeq and gkeq are written in the Taylor expansion as follows:(4a)fkeqi,j=wkρi,j1+3ck⋅uc2+92ck⋅u2c4−32u2c2,(4b)gkeqi,j=wkTi,j1+3ck⋅uc2+92ck⋅u2c4−32u2c2,where the weight factors wk are w0=4/9, w1−4=1/9, and w5−8=1/36 and c is the lattice speed given by c=Δx/Δt=Δy/Δt.The discrete scheme D2Q9 is shown in Figure2(5)ck=0,0,k=0,ccosk−1π2,sink−1π2,k=1−4,2ccosk−5π2+π4,sink−5π2+π4,k=5−8.Figure 2
D2Q9 velocities.The macroscopic densityρ, velocity u, and internal energy by unit of mass e are obtained by the following relations:(6a)ρ=∑k=08fk,(6b)ρu=∑k=08fkck,(6c)ρe=∑k=08gk,where e=DRT/2, in which D is the number of physical dimensions (equal to 2 in the current work).


## 3.2. Boundary Conditions

### 3.2.1. Flow Boundary Conditions

The boundary conditions proposed by Zou and He [23] are used at the inlet, in which the normal velocity component is assumed to be equal to zero and the density is to be determined. After the streaming step, the unknown distribution functions are at the inlet boundary f1,f5,f8. Using equations (6a) and (6b), these function densities are calculated as follows:(7a)ρin=f0+f2+f4+2f3+f6+f71−Uin,(7b)f1=f3+23ρinUin,(7c)f5=f7−12f2−f4+16ρinUin,(7d)f8=f6+12f2−f4+16ρinUin.A simple extrapolation scheme is used for the velocity at the outlet.To capture velocity slip phenomenon, the Maxwell first order of slip boundary condition is applied in its dimensionless form, near the longitudinal walls [24]:(8a)ux,y=0slip=ux,y=0−ux,wall=σKn∂u∂yy=0,(8b)ux,y=Hslip=ux,wall−ux,y=H=σKn∂u∂yy=H,where ux,y=0slip and ux,y=Hslip are the slip velocities at the bottom and top walls. The walls are considered at rest ux,wall=0 and σ is the momentum accommodation coefficient which is assumed to be equal unity, to simulate completely diffuse reflection [33].At the top wall, the unknown distribution functionsf4,f7,f8 are calculated as follows:(9a)ux,y=H=λ4ux,H−1−ux,H−22+3λ,(9b)ρtop=f0+f1+f3+2f2+f6+f5,(9c)f4=f2,(9d)f7=ρtop1−ux,y=H−f0+f2+f42−f3+f6,(9e)f8=ρtop1+ux,y=H−f0+f2+f42−f1+f5,where the gas mean free path is defined as λ=KnH and H is the total number of lattice nodes in the vertical direction.


### 3.2.2. Temperature Boundary Conditions

At the inlet boundary, the unknown energy distribution functionsg1,g5,g8 are obtained as follows:(10a)g1=Tinletw1+w3−g3,(10b)g5=Tinletw5+w7−g7,(10c)g8=Tinletw6+w8−g6.The outlet boundary condition for energy distributions is implemented by using an extrapolation scheme [34]. Similarly, the temperature jump boundary conditions caused by the rarefaction effects at the walls are given by(11a)Tx,y=0jump=Tx,y=0−Tx,wall=ϕKn2γγ+1Pr∂T∂yy=0,(11b)Tx,y=Hjump=Tx,wall−Tx,y=H=ϕKn2γγ+1Pr∂T∂yy=H.C
jump is the temperature jump coefficient defined as follows:(12)Cjump=ϕ2γγ+1Prλ=kλ.In which,ϕ represents the thermal accommodation coefficient and is assumed to be unity, γ is the specific heat ratio, and Pr is the Prandtl number.Temperature jump boundary conditions are written at the bottom and top walls as follows:(13a)Tx,y=0=Cjump4Tx,1−Tx,2+2Tx,wall2+3Cjump,(13b)Tx,y=H=Cjump4Tx,H−1−Tx,H−2+2Tx,wall2+3Cjump.At the top wall, the unknown distribution functionsg4,g7,g8 are calculated as follows:(14a)Tx,y=H=Cjump4Tx,H−1−Tx,H−2+2Tx,wall2+3Cjump,(14b)g4=Tx,y=Hw2+w4−g2,(14c)g7=Tx,y=Hw5+w7−g5,(14d)g8=Tx,y=Hw6+w8−g6.


## 3.2.1. Flow Boundary Conditions

The boundary conditions proposed by Zou and He [23] are used at the inlet, in which the normal velocity component is assumed to be equal to zero and the density is to be determined. After the streaming step, the unknown distribution functions are at the inlet boundary f1,f5,f8. Using equations (6a) and (6b), these function densities are calculated as follows:(7a)ρin=f0+f2+f4+2f3+f6+f71−Uin,(7b)f1=f3+23ρinUin,(7c)f5=f7−12f2−f4+16ρinUin,(7d)f8=f6+12f2−f4+16ρinUin.A simple extrapolation scheme is used for the velocity at the outlet.To capture velocity slip phenomenon, the Maxwell first order of slip boundary condition is applied in its dimensionless form, near the longitudinal walls [24]:(8a)ux,y=0slip=ux,y=0−ux,wall=σKn∂u∂yy=0,(8b)ux,y=Hslip=ux,wall−ux,y=H=σKn∂u∂yy=H,where ux,y=0slip and ux,y=Hslip are the slip velocities at the bottom and top walls. The walls are considered at rest ux,wall=0 and σ is the momentum accommodation coefficient which is assumed to be equal unity, to simulate completely diffuse reflection [33].At the top wall, the unknown distribution functionsf4,f7,f8 are calculated as follows:(9a)ux,y=H=λ4ux,H−1−ux,H−22+3λ,(9b)ρtop=f0+f1+f3+2f2+f6+f5,(9c)f4=f2,(9d)f7=ρtop1−ux,y=H−f0+f2+f42−f3+f6,(9e)f8=ρtop1+ux,y=H−f0+f2+f42−f1+f5,where the gas mean free path is defined as λ=KnH and H is the total number of lattice nodes in the vertical direction.


## 3.2.2. Temperature Boundary Conditions

At the inlet boundary, the unknown energy distribution functionsg1,g5,g8 are obtained as follows:(10a)g1=Tinletw1+w3−g3,(10b)g5=Tinletw5+w7−g7,(10c)g8=Tinletw6+w8−g6.The outlet boundary condition for energy distributions is implemented by using an extrapolation scheme [34]. Similarly, the temperature jump boundary conditions caused by the rarefaction effects at the walls are given by(11a)Tx,y=0jump=Tx,y=0−Tx,wall=ϕKn2γγ+1Pr∂T∂yy=0,(11b)Tx,y=Hjump=Tx,wall−Tx,y=H=ϕKn2γγ+1Pr∂T∂yy=H.C
jump is the temperature jump coefficient defined as follows:(12)Cjump=ϕ2γγ+1Prλ=kλ.In which,ϕ represents the thermal accommodation coefficient and is assumed to be unity, γ is the specific heat ratio, and Pr is the Prandtl number.Temperature jump boundary conditions are written at the bottom and top walls as follows:(13a)Tx,y=0=Cjump4Tx,1−Tx,2+2Tx,wall2+3Cjump,(13b)Tx,y=H=Cjump4Tx,H−1−Tx,H−2+2Tx,wall2+3Cjump.At the top wall, the unknown distribution functionsg4,g7,g8 are calculated as follows:(14a)Tx,y=H=Cjump4Tx,H−1−Tx,H−2+2Tx,wall2+3Cjump,(14b)g4=Tx,y=Hw2+w4−g2,(14c)g7=Tx,y=Hw5+w7−g5,(14d)g8=Tx,y=Hw6+w8−g6.


## 3.3. Nusselt Number Calculation

Heat transfer characteristic of the flow can be determined using the Nusselt numberNu, which is the ratio of convective and conductive heat transfers. The Nusselt number along the horizontal axis is calculated as follows:(15)Nux=2H∂T/∂ywallTwall−Tbulk.The bulk temperature at a cross-section is calculated as follows:(16)Tbulk=∫SρuTdS∫SρudS.The average Nusselt number is calculated as follows:(17)Nuavg=∫0LNuxdxL.


## 4. Finite Difference Method

The advection-diffusion equation in 2D is expressed by(18)∂T∂t+u∂T∂x+v∂T∂y=α∂T2∂x2+∂T2∂y2,where α is the coefficient of thermal diffusion and u,v are the components of velocity.For a steady Poiseuille flow, the vertical velocity component is omittedv=0 and the following theoretical velocity profile u across a microchannel is used in the finite difference method [8]:(19)uy=6Umeany/H−y2/H2+2−σ/σKn1+62−σ/σKn,where σ is the momentum accommodation coefficient.The domain is discretized into equal segmentsΔx=Δy and the finite difference approach (FDM) is used, and equation (18) becomes(20)Tn+1i,j=Tni,j1−1τg+1τgTni+1,j+Tni−1,j+Tni,j+1+Tni,j−12−uΔtTni,j−Tni−1,jΔx,where n is the number of the time step Δt and 1/τg=2αΔt/Δx2.The term1−1/τg must be positive, which implies that Δt≤Δx2/2α.Temperature jump boundary conditions used in TLBM approach are also used in the FDM approach at the bottom and the top walls. In the inlet, the temperature is taken equal to unity and a simple extrapolation can be used at the outlet.


## 5. Results

In the present study, all simulations are performed using a developed Fortran code. A gas flow between two parallel plates at rest has been simulated using the thermal lattice Boltzmann method and finite difference approaches. Two cases of micro-Poiseuille flow are treated; in the first case, the temperature of the upper and lower plates is taken equal toTC; however, a linear decreasing of plate temperature is imposed in the second case. To describe the results of simulation easily, it is more convenient to use the following normalization for temperature: θ=T−TC/TH−TC. Thus, the temperatures TH, TC, and Tp become θH=1, θC=0, and θp=1/p. In this study, we take Pr=0.71, k = 1.667, and Reynolds number is fixed at Re=10.


### 5.1. Mesh Independence Study

The channel aspect ratio is fixed in this study to beAR=4. For Kn=0.05 and θp=0.5, Table 1 shows the effect of the mesh on the average Nusselt number, while Table 2 shows the effect of the mesh on the horizontal velocity and the temperature near the inlet x/L=0.04,y/H=0.5, (not to mention that θp does not affect the velocity). In this study, a mesh of 200 × 50 has been used for numerical investigation because it gives a stable solution.Table 1
Effect of the mesh on the average Nusselt number forKn=0.05 and θp=0.5.
Mesh
180× 45
200× 50
220× 55
Nu
6.78356
6.83014
6.86398Table 2
Effect of the mesh on the velocity and temperature.
Mesh
180× 45
200× 50
220× 55
u
/
U
in
1.05
1.04
1.04
θ
0.974
0.977
0.977


### 5.2. Numerical Validation

To ensure the validation of the present model, the first case is treated and the results found have been compared with the FD method and the results presented in Zarita and Hachemi [8]. To observe the rarefaction effect on the gas flow behavior, the velocity profiles for Kn = 0.01 and 0.08 (Figure 3(a)) and the temperature profiles for Kn = 0.01, 0.05 and 0.08 (Figure 3(b)) are plotted. Figures show the effect of Kn on the velocity slip and temperature jump at the centerline of the microchannel. Good agreement between the analytical solution and the numerical one obtained by TLBM for velocity is observed. In agreement with kinetic theory, it is shown that velocity slip is sensitive, inversion of velocity profiles, to the rarefaction degree near the longitudinal walls, the Knudsen layer y/L<0.2andy/L>0.8 [6, 8–10, 22]. To compare the transient time of both solutions given by FDM and TLBM approaches, the temperature profile is plotted as a function of time steps number for Kn=0.01−0.08. For a given Kn, the FDM and the TLBM require the same time to reach the steady state, and with increasing Kn, the temperature value in the equilibrium state at the bulk increases. This is due to the interactions of particles which become more important (Figure 4). The convergence of the Nusselt number to a constant value is reached as soon as the velocity and the temperature become constant (Figure 5). Both methods give the same temperature contours for Kn=0.01 and Kn=0.08. The higher values are obtained near the inlet and when Kn increases the values of the temperature increases at the outlet of the microchannel (see Figure 6). So, it is observed that the rarefaction effect increases the convection from the hot inlet stream for both approaches.Figure 3
(a) Velocity profiles and (b) temperature profiles along the vertical centerline axis according toKn.
(a)
(b)Figure 4
Evolution of the temperature at the center of the microchannel.Figure 5
Nusselt number variation along the microchannel according toKn.Figure 6
Normalised temperature isolines in the microchannel: (a) TLBM-Kn = 0.01, (b) TLBM-Kn = 0.08, (c) FDM-Kn = 0.01, and (d) FDM-Kn = 0.08.
(a)
(b)
(c)
(d)


### 5.3. Numerical Results and Analysis

The second case is evaluated to investigate the effect of the nonisothermal walls on the gas flow. The walls are not isothermal because of the effect of the conduction of the inlet temperatureθinlet=1, and this effect is assumed to be linear in this study. To compare the transient time of both solutions given by FDM and TLBM approaches for θp=0.1 and θp=0.8, the temperature profile at the microchannel center is plotted as a function of time steps number for Kn=0.01−0.08 (Figure 7). By comparing Figures 7(a) and 7(b), the effect of the Knudsen number decreases when θp has a large value. The temperature reaches the steady state after approximately 7500 number time steps, and its value in the equilibrium state increases by increasing Kn and θp. Figure 7(c) shows the effect of θp on the temperature value, which increases with increasing θp, at the center of microchannel for Kn=0.05. The temperature jump increases with Kn and θp [32] (see Table 3). By increasing the rarefaction degree Kn, its effect on the temperature profile along the horizontal axis near the wall y/H=0.01 for θp=0.1 and x/L<0.4 is important (Figures 8(a) and 8(c)), whereas its effect vanishes when θp approaches to unity and its profile becomes almost linear (Figures 8(b) and 8(d)). After that, inlet excitation, both approaches predict a decrease in temperature to reach the cold value at the outlet. In the vertical direction, the temperature profiles obtained by both approaches of TLBM and FDM, for different values of Kn number, are plotted as a function of y-coordinate. In addition to the good agreement between the results, the temperature is sensitive to the rarefaction degree Kn. The temperature jump and the amplitude of its curve near the channel inlet x/L=0.01andx/L=0.05 along the vertical axis increase as Kn increases (Figure 9). By increasing θp, the profile of the temperature along the vertical axis near the center of the channel x/L=0.25andx/L=0.5 begins to lose its parabolic appearance and take a horizontal shape even for the low values of Kn (Figures 9(b) and 9(d)). This is since the gas becomes more rarefied when θp has a great value. Both cases give the same velocity profiles since equation (6b) which gives the velocity is independent of the temperature (see Table 4). To compare the transient time of the solution given by the TLBM approach, the evolution of the horizontal velocity at the center of microchannel is plotted as a function of time stepsnumber for Kn=0.01−0.08. This approach reaches the steady state after approximately 2000 time steps and the effect of the Knudsen number on the velocity is seen because by increasing Kn, the velocity decreases at the center of the microchannel (Figure 10(a)). Using the slip boundary condition, the so-called velocity slip phenomenon which appears under rarefaction effects, for Kn=0.01−0.08, is captured by the TLB method (Figures 10(b)–10(d)) [6, 22]. Near the entrance region, the flow is developed fast; therefore, the Knudsen number has a significant effect on the Nusselt number (see Figure 11). The effect of nonisothermal walls for θp = 0.5 on the temperature in the microchannel is clearly visible and the values of the temperature contours are strictly greater than 0. It is observed that the rarefaction effect increases the convection from the hot inlet stream for both methods, and when Kn increases, the small values of the temperature vanish at the outlet of the microchannel (see Figure 12).Figure 7
Evolution of the temperature at the center of the microchannel as a function of time: (a) TLBM and FDM forθp = 0.1, (b) TLBM and FDM for θp = 0.8, and (c) TLBM for Kn = 0.05.
(a)
(b)
(c)Table 3
Effect ofKn and θp on wall temperature jump.
θ
p
0 (case 1)
0.1
0.5
0.8
1
Kn
0.01
0.00844
0.00888
0.01063
0.01195
0.01283
0.03
0.02754
0.02854
0.03254
0.03554
0.03753
0.05
0.04902
0.05020
0.05491
0.05844
0.06079
0.08
0.08417
0.08504
0.08851
0.09112
0.09285Figure 8
Temperature profile along the horizontal axis according toθp and Kn obtained using TLB and FD methods for (a) θp=0.1 and Kn=0.01, (b) θp=0.8 and Kn=0.01, (c) θp=0.1 and Kn=0.08, and (d) θp=0.8 and Kn=0.08.
(a)
(b)
(c)
(d)Figure 9
Temperature profile along the vertical axis according toθp and Kn obtained using TLB and FD methods for (a) θp=0.1 and Kn=0.01, (b) θp=0.8 and Kn=0.01, (c) θp=0.1 and Kn=0.08, and (d) θp=0.8 and Kn=0.08.
(a)
(b)
(c)
(d)Table 4
Velocity and temperature values atx/L=0,y/H=0 for Kn=0.05 obtained using the TLB method.
θ
p
0 (case 1)
0.5
u
/
U
in
0.78947
0.78947
θ
0.78445
0.89175Figure 10
(a) Evolution of the velocity at the center of the microchannel as a function of time and velocity profiles for (b)Kn=0.01, (c) Kn=0.05, and (d) Kn=0.08.
(a)
(b)
(c)
(d)Figure 11
Nusselt number variation along the microchannel for (a)θp=0.8 and (b) Kn=0.05.
(a)
(b)Figure 12
Normalised temperature isolines in the microchannel forθp=0.5 (a) TLBM-Kn=0.01, (b) TLBM-Kn=0.08, (c) FDM-Kn=0.01, and (d) FDM-Kn=0.08.
(a)
(b)
(c)
(d)The temperature jump is calculated at the center of the wall as follows:(21)θjump=θxL=0.5−θwall.For both cases, by increasing the Knudsen number value, the average value of the Nusselt number decreases and in parallel whenθp increases its value increases, but by increasing this value to θp = 1, the results become insignificant for that number (Table 5).Table 5
Value of the average Nusslet number at differentθp and at different Kn.
θ
p
0 (case 1)
0.1
0.5
0.8
1
Kn
0.01
8.34289
8.46053
8.70021
8.84516
6.87401
0.03
7.35054
7.43502
7.62395
7.73609
6.50464
0.05
6.62209
6.68289
6.83014
6.91603
6.00751
0.08
5.76999
5.80805
5.90945
5.96658
5.31223In the TLBM approach, the flow boundary conditions proposed by Zou and He [23] are used at the inlet, and these conditions are constructed in consistence with the wall boundary condition based on an idea of bounce-back of nonequilibrium distribution functions given by second-order accurate of equations (7a)–(7d) and for energy distribution function equations (10a)–(10c) are used. However, in the FD method, the analytical solution of theoretical velocity (equation (19)) is used while the inlet temperature is taken equal to unity, that is why the methods are quite different for y/H=0.01 and y/H=0 at the inlet x/L=0 (see Figure 8 and Table 4).To sum up, the TLBM method is able to simulate micro-Poiseuille flow in the case of nonisothermal walls, and it can capture the slip velocity and jump temperature at the wall even at the inlet.


## 5.1. Mesh Independence Study

The channel aspect ratio is fixed in this study to beAR=4. For Kn=0.05 and θp=0.5, Table 1 shows the effect of the mesh on the average Nusselt number, while Table 2 shows the effect of the mesh on the horizontal velocity and the temperature near the inlet x/L=0.04,y/H=0.5, (not to mention that θp does not affect the velocity). In this study, a mesh of 200 × 50 has been used for numerical investigation because it gives a stable solution.Table 1
Effect of the mesh on the average Nusselt number forKn=0.05 and θp=0.5.
Mesh
180× 45
200× 50
220× 55
Nu
6.78356
6.83014
6.86398Table 2
Effect of the mesh on the velocity and temperature.
Mesh
180× 45
200× 50
220× 55
u
/
U
in
1.05
1.04
1.04
θ
0.974
0.977
0.977


## 5.2. Numerical Validation

To ensure the validation of the present model, the first case is treated and the results found have been compared with the FD method and the results presented in Zarita and Hachemi [8]. To observe the rarefaction effect on the gas flow behavior, the velocity profiles for Kn = 0.01 and 0.08 (Figure 3(a)) and the temperature profiles for Kn = 0.01, 0.05 and 0.08 (Figure 3(b)) are plotted. Figures show the effect of Kn on the velocity slip and temperature jump at the centerline of the microchannel. Good agreement between the analytical solution and the numerical one obtained by TLBM for velocity is observed. In agreement with kinetic theory, it is shown that velocity slip is sensitive, inversion of velocity profiles, to the rarefaction degree near the longitudinal walls, the Knudsen layer y/L<0.2andy/L>0.8 [6, 8–10, 22]. To compare the transient time of both solutions given by FDM and TLBM approaches, the temperature profile is plotted as a function of time steps number for Kn=0.01−0.08. For a given Kn, the FDM and the TLBM require the same time to reach the steady state, and with increasing Kn, the temperature value in the equilibrium state at the bulk increases. This is due to the interactions of particles which become more important (Figure 4). The convergence of the Nusselt number to a constant value is reached as soon as the velocity and the temperature become constant (Figure 5). Both methods give the same temperature contours for Kn=0.01 and Kn=0.08. The higher values are obtained near the inlet and when Kn increases the values of the temperature increases at the outlet of the microchannel (see Figure 6). So, it is observed that the rarefaction effect increases the convection from the hot inlet stream for both approaches.Figure 3
(a) Velocity profiles and (b) temperature profiles along the vertical centerline axis according toKn.
(a)
(b)Figure 4
Evolution of the temperature at the center of the microchannel.Figure 5
Nusselt number variation along the microchannel according toKn.Figure 6
Normalised temperature isolines in the microchannel: (a) TLBM-Kn = 0.01, (b) TLBM-Kn = 0.08, (c) FDM-Kn = 0.01, and (d) FDM-Kn = 0.08.
(a)
(b)
(c)
(d)


## 5.3. Numerical Results and Analysis

The second case is evaluated to investigate the effect of the nonisothermal walls on the gas flow. The walls are not isothermal because of the effect of the conduction of the inlet temperatureθinlet=1, and this effect is assumed to be linear in this study. To compare the transient time of both solutions given by FDM and TLBM approaches for θp=0.1 and θp=0.8, the temperature profile at the microchannel center is plotted as a function of time steps number for Kn=0.01−0.08 (Figure 7). By comparing Figures 7(a) and 7(b), the effect of the Knudsen number decreases when θp has a large value. The temperature reaches the steady state after approximately 7500 number time steps, and its value in the equilibrium state increases by increasing Kn and θp. Figure 7(c) shows the effect of θp on the temperature value, which increases with increasing θp, at the center of microchannel for Kn=0.05. The temperature jump increases with Kn and θp [32] (see Table 3). By increasing the rarefaction degree Kn, its effect on the temperature profile along the horizontal axis near the wall y/H=0.01 for θp=0.1 and x/L<0.4 is important (Figures 8(a) and 8(c)), whereas its effect vanishes when θp approaches to unity and its profile becomes almost linear (Figures 8(b) and 8(d)). After that, inlet excitation, both approaches predict a decrease in temperature to reach the cold value at the outlet. In the vertical direction, the temperature profiles obtained by both approaches of TLBM and FDM, for different values of Kn number, are plotted as a function of y-coordinate. In addition to the good agreement between the results, the temperature is sensitive to the rarefaction degree Kn. The temperature jump and the amplitude of its curve near the channel inlet x/L=0.01andx/L=0.05 along the vertical axis increase as Kn increases (Figure 9). By increasing θp, the profile of the temperature along the vertical axis near the center of the channel x/L=0.25andx/L=0.5 begins to lose its parabolic appearance and take a horizontal shape even for the low values of Kn (Figures 9(b) and 9(d)). This is since the gas becomes more rarefied when θp has a great value. Both cases give the same velocity profiles since equation (6b) which gives the velocity is independent of the temperature (see Table 4). To compare the transient time of the solution given by the TLBM approach, the evolution of the horizontal velocity at the center of microchannel is plotted as a function of time stepsnumber for Kn=0.01−0.08. This approach reaches the steady state after approximately 2000 time steps and the effect of the Knudsen number on the velocity is seen because by increasing Kn, the velocity decreases at the center of the microchannel (Figure 10(a)). Using the slip boundary condition, the so-called velocity slip phenomenon which appears under rarefaction effects, for Kn=0.01−0.08, is captured by the TLB method (Figures 10(b)–10(d)) [6, 22]. Near the entrance region, the flow is developed fast; therefore, the Knudsen number has a significant effect on the Nusselt number (see Figure 11). The effect of nonisothermal walls for θp = 0.5 on the temperature in the microchannel is clearly visible and the values of the temperature contours are strictly greater than 0. It is observed that the rarefaction effect increases the convection from the hot inlet stream for both methods, and when Kn increases, the small values of the temperature vanish at the outlet of the microchannel (see Figure 12).Figure 7
Evolution of the temperature at the center of the microchannel as a function of time: (a) TLBM and FDM forθp = 0.1, (b) TLBM and FDM for θp = 0.8, and (c) TLBM for Kn = 0.05.
(a)
(b)
(c)Table 3
Effect ofKn and θp on wall temperature jump.
θ
p
0 (case 1)
0.1
0.5
0.8
1
Kn
0.01
0.00844
0.00888
0.01063
0.01195
0.01283
0.03
0.02754
0.02854
0.03254
0.03554
0.03753
0.05
0.04902
0.05020
0.05491
0.05844
0.06079
0.08
0.08417
0.08504
0.08851
0.09112
0.09285Figure 8
Temperature profile along the horizontal axis according toθp and Kn obtained using TLB and FD methods for (a) θp=0.1 and Kn=0.01, (b) θp=0.8 and Kn=0.01, (c) θp=0.1 and Kn=0.08, and (d) θp=0.8 and Kn=0.08.
(a)
(b)
(c)
(d)Figure 9
Temperature profile along the vertical axis according toθp and Kn obtained using TLB and FD methods for (a) θp=0.1 and Kn=0.01, (b) θp=0.8 and Kn=0.01, (c) θp=0.1 and Kn=0.08, and (d) θp=0.8 and Kn=0.08.
(a)
(b)
(c)
(d)Table 4
Velocity and temperature values atx/L=0,y/H=0 for Kn=0.05 obtained using the TLB method.
θ
p
0 (case 1)
0.5
u
/
U
in
0.78947
0.78947
θ
0.78445
0.89175Figure 10
(a) Evolution of the velocity at the center of the microchannel as a function of time and velocity profiles for (b)Kn=0.01, (c) Kn=0.05, and (d) Kn=0.08.
(a)
(b)
(c)
(d)Figure 11
Nusselt number variation along the microchannel for (a)θp=0.8 and (b) Kn=0.05.
(a)
(b)Figure 12
Normalised temperature isolines in the microchannel forθp=0.5 (a) TLBM-Kn=0.01, (b) TLBM-Kn=0.08, (c) FDM-Kn=0.01, and (d) FDM-Kn=0.08.
(a)
(b)
(c)
(d)The temperature jump is calculated at the center of the wall as follows:(21)θjump=θxL=0.5−θwall.For both cases, by increasing the Knudsen number value, the average value of the Nusselt number decreases and in parallel whenθp increases its value increases, but by increasing this value to θp = 1, the results become insignificant for that number (Table 5).Table 5
Value of the average Nusslet number at differentθp and at different Kn.
θ
p
0 (case 1)
0.1
0.5
0.8
1
Kn
0.01
8.34289
8.46053
8.70021
8.84516
6.87401
0.03
7.35054
7.43502
7.62395
7.73609
6.50464
0.05
6.62209
6.68289
6.83014
6.91603
6.00751
0.08
5.76999
5.80805
5.90945
5.96658
5.31223In the TLBM approach, the flow boundary conditions proposed by Zou and He [23] are used at the inlet, and these conditions are constructed in consistence with the wall boundary condition based on an idea of bounce-back of nonequilibrium distribution functions given by second-order accurate of equations (7a)–(7d) and for energy distribution function equations (10a)–(10c) are used. However, in the FD method, the analytical solution of theoretical velocity (equation (19)) is used while the inlet temperature is taken equal to unity, that is why the methods are quite different for y/H=0.01 and y/H=0 at the inlet x/L=0 (see Figure 8 and Table 4).To sum up, the TLBM method is able to simulate micro-Poiseuille flow in the case of nonisothermal walls, and it can capture the slip velocity and jump temperature at the wall even at the inlet.


## 6. Conclusion

In this study, the efficiency of the TLBM method to simulate micro-Poiseuille flow in the case of nonisothermal walls is demonstrated. Slip and jump boundary conditions (SJBC) are used to capture the nonequilibrium effect near the walls. The explored Knudsen numbers correspond to the slip regime. The Nusselt number which depends on the velocity and temperature profiles proves the effectiveness of the results. Good agreement is obtained between TLBM and FDM results. So, regarding its fast convergence, it shows the ability of TLBM to describe the velocity slip and temperature jump as a good alternative that can be used to describe the gas microflows usually encountered in the micro-electro-mechanical systems (MEMS) and nano-electro-mechanical systems (NEMS) devices. Unlike the kinetic methods, such as DSMC, which need a long time to achieve satisfactory results, the TLBM is more promising for such flows.



---

*Source: 8638494-2020-03-19.xml*