Thermal Lattice BoltzmannModel for Nonisothermal Gas Flow in a Two-Dimensional Microchannel

In this paper, the (ermal Lattice Boltzmann Method (TLBM) is used for the simulation of a gas microflow. A 2D heated microchannel flow driven by a constant inlet velocity profile Uin 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. (e rarefaction effects described by the Knudsen number, on the velocity and temperature profiles are evaluated. (e aim of this study is to prove the efficiency of the TLBMmethod 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). (e results also show an interesting sensitivity of velocity and temperature profiles with the rarefaction degree and the imposed temperature gradient of the walls.


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][5][6]. is method is used to simulate different types of flows: the Rayleigh-Benard convection [7], micro-Poiseuille flow [8][9][10][11][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 r f 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 bounceback 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. e LBM approximation is a linear discretization of the Boltzmann equation. e 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. e miniaturization of devices, MEMS technology, leads to the microsystems where the gas particles mean free path is comparable to the system characteristic length. erefore, 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. e 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. e 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][9][10][11][27][28][29][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. e 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.

Statement of Problem
In the present simulation, we studied forced convection inside a long microchannel in the case of the body force is zero. e 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 direction z [32]. e 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 T p (T p � T H /p) to cold one T C , so T p varies from T H to T C (p is an integer ≥ 1). e variation of walls temperature in the second case is as follows:

Lattice Boltzmann Method.
A square grid and D2Q9 model is used for both distribution function density f and temperature g. e governing equations for these distribution functions are written according to the BGK model as follows: 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 e macroscopic density ρ, velocity u, and internal energy by unit of mass e are obtained by the following relations: where e � DRT/2, in which D is the number of physical dimensions (equal to 2 in the current work).

Flow Boundary
Conditions. e 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 (f 1 , f 5 , f 8 ). Using equations (6a) and (6b), these function densities are calculated as follows: 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]: where u slip x,y�0 and u slip x,y�H are the slip velocities at the bottom and top walls. e walls are considered at rest (u x,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 functions (f 4 , f 7 , f 8 ) are calculated as follows: where the gas mean free path is defined as λ � KnH and H is the total number of lattice nodes in the vertical direction.

Temperature Boundary Conditions.
At the inlet boundary, the unknown energy distribution functions (g 1 , g 5 , g 8 ) are obtained as follows: e 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 C jump is the temperature jump coefficient defined as follows: In which, ϕ represents the thermal accommodation coefficient and is assumed to be unity, c 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: At the top wall, the unknown distribution functions (g 4 , g 7 , g 8 ) are calculated as follows:

Nusselt Number Calculation.
Heat transfer characteristic of the flow can be determined using the Nusselt number Nu, which is the ratio of convective and conductive heat transfers. e Nusselt number along the horizontal axis is calculated as follows: e bulk temperature at a cross-section is calculated as follows: e average Nusselt number is calculated as follows:

Finite Difference Method
e advection-diffusion equation in 2D is expressed by 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 omitted (v � 0) and the following theoretical velocity profile u across a microchannel is used in the finite difference method [8]: where σ is the momentum accommodation coefficient. e domain is discretized into equal segments (Δx � Δy) and the finite difference approach (FDM) is used, and equation (18) becomes where n is the number of the time step Δt and (1/τ g ) � ((2αΔt)/Δx 2 ). e term (1 − (1/τ g )) must be positive, which implies that Δt ≤ (Δx 2 /(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.

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 to T C ; 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 − T C )/(T H − T C )).
us, the temperatures T H , T C , and T p 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.

Mesh Independence Study.
e channel aspect ratio is fixed in this study to be AR � 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.

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.2 and (y/L) > 0.8) [6,[8][9][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. is is due to the interactions of particles which become more important (Figure 4). e 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. e 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.

Numerical Results and Analysis.
e second case is evaluated to investigate the effect of the nonisothermal walls on the gas flow. e 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. e 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. e 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 Mathematical Problems in Engineering     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. e temperature jump and the    (Figure 9). By increasing θ p , the profile of the temperature along the vertical axis near the center of the channel (x/L) � 0.25 and (x/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)). is 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. is 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   [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). e 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). e temperature jump is calculated at the center of the wall as follows: 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).
In 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 bounceback 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.

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. e explored Knudsen numbers correspond to the slip regime. e 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