Analysis of the slip flow in the hydrodynamic entrance region of a 2D microchannel

Two-dimensional developing flow in the entrance of a microchannel has been studied numerically. Due to its nature, a microchannel can be used in tight space applications and the length of channel can get very small values. Furthermore, if the hydrodynamic development length of flow in microchannel has comparably the same value with the channel length, the channel entrance parameters play critical role on the flow performance of a microscale channel. Lattice Boltzmann Method (LBM) was considered for studying and simulating the developing slip flows through a rectangular microchannel. A unique computational code for this study was developed by using LBM. The slip velocity boundary condition along with Knudsen number values in the slip flow regime was used for this model. The bounce-back boundary condition was considered at the top and bottom walls of the microchannel. The effects of the Reyn-olds numbers (1-100) and Knudsen numbers (0.001, 0.01, 0.1) on the hydrodynamic entrance length has been examined. The numerical results have been compared with the previous studies in the literature and the similarities have been found satisfactory. The entrance length is found to be increasing with the increase of Reynolds and Knudsen numbers. A correlation for hydrodynamic development length as a function of Knudsen and Reynolds numbers was obtained by using the data extracted from LBM simulations performed in this study.


INTRODUCTION
As the demand on microchannels increases gradually from the need for use in new technologies, developments in manufacturing technologies of microscale channels increase as well. Therefore, the research of the effect of various cross sections and lengths of channels on the flow characteristics of fluid in microscale channels gains more importance. The velocity profiles are found to have fast changes from uniform to parabolic geometry in this short channel lengths.
Microscale channels can be seen in different application areas in Micro-Electro-Mechanical-Systems (MEMS) [1][2][3], where their sizes are in micron scale. The attractive side of MEMS is due to their miniaturization and installation on tight spaces in microelectronics. The increased developments found in microscale devices have been carried on these devices to address the most challenging tight space applications. Further enhancements made in the production capability of the micro scale devices allow these devices to be available in smaller sizes.
The microscale devices can be used in gas and/or liquid application areas. The applications can be seen in defense, transportation, medicine, computer and space industry [4,5]. Some of the MEMS have internal flow passages that include fluid flows in these micro channels, such structures can be seen in lab-on-a-chip devices, which generally found application in automotive, aerospace, medical and military areas [6]. The behavior of the flows in such micro scale systems is often different from what we encountered in daily life. It may be necessary to find different study technics for an efficient modelling of fluid flows in microscale channels. The flow modelling is important to make efficient designs by considering its cost and performance especially in high volume productions. Before performing the experiments and manufacturing the prototype products, the simulation of the system can be time and cost effective. Comparing the macro scale devices, conducting a test with a microscale device is not easy and does not always give reliable results. Due to dimensions in 0.1-10 μm sizes, it is not efficiently tested comparing to macro-scale devices and it is generally quite difficult to conduct a test and get reliable results for microscale devices.
These systems have generally seen with gaseous flows and the behavior of the gaseous flow defines the performance of a microchannel system [7,8]. The ratio of the gas mean free path (λ) to the characteristic length of flow field (L) is known as the Knudsen number (Kn) a dimensionless parameter. Kn number classifies the flows according to that continuum or free molecular flow. And with this classification, a suitable solution approach can be selected for the flow modelling. The mean free path (λ) is the distance that a gas particle must travel between two consecutive collisions with other particles. Different mathematical approaches have been developed to study the behavior of fluid flows in different flow regimes. The general Navier-Stokes equation can be used to study the fluid flow and the non-slip boundary condition can be applied on wall boundaries when Kn ≤ 0.001. As the size of the flow domain is reduced, Kn increases further and the rarefaction affects the flow field as a result the velocity slip is going to be experienced on boundaries of channel [9]. In the slip flow regime where Knudsen number is equal or between the values of 0.001 and 0.1, the Navier-Stokes equation can still be used to solve the considered flow domain along with the slip conditions that must be taken into account on wall boundaries. In the transition (0.1 ≤ Kn ≤ 10) and free molecule (Kn > 10) flow regimes, the continuum assumption cannot be used and some other two approaches, Molecular Dynamics (MD) or Direct Simulation Monte Carlo (DSMC) methods must be used to solve the flow domain [10,11]. Molecular Dynamics may not be a good choice for the simulations because of the restriction in both length and time scales due to computational reasons [12].
The Boltzmann Equation (BE) can cover all of the three flow regimes mentioned above. The BE can be used to describe the velocity distribution function along with molecular transport and collisions. The solution of BE is not easy so the additional approaches can be used for the BE to approximate. One of these approaches, Lattice Boltzmann Method (LBM), used widely especially in microchannel flows and is gaining popularity nowadays. In LBM, the discrete-velocity Boltzmann equation is solved with the velocities on a discrete lattice [13].
The lattice Boltzmann method (LBM) has been found as an alternative and useful tool for studying the micro scale flows in channels [14]. Some of its advantages over the other solution methods are; the comparable computational cost of LBM over the Navier-Stokes solvers [15], being faster than other two solution methods (MD and DSMC) [16], and the application capability in a wide flow regimes compared to Navier-Stokes solvers [17][18][19]. Implementation of slip boundary condition is an important feature of LBM for microflows. There are different boundary conditions used for the model [20]; but the Maxwell slip model can be found as the most used one. In this model, velocity slip is described along with Knudsen number, velocity gradient and accommodation coefficient [21].
A wide variety of fluid flow cases in microchannels have been studied by using LBM with adequate results. Arkilic et al. [22,23] conducted experimental and analytical studies with gaseous flow in microchannels. Nie at al. [24] studied on two dimensional micro channel flows by using LBM with bounce back boundary condition. Isfahani and et al. [25][26][27] studied the flow characteristics of the nanofluids in porous media experimentally and numerically. Tang at al. [28] used a boundary condition based on the kinetic theory to determine a gaseous slip flow. Zhang at al. [20] considered the accommodation coefficient for flow surface interactions. Shirani and Jafari [29] used the bounce-back boundary condition. Myong et al. [30] studied the flow by using Langmuir slip model. Kim at al. [31] used the Langmuir slip model. Chen et al. [32] studied the compressible Navier-Stokes equations by using the velocity slip boundary condition. Renksizbulut et al. [33] studied the behavior of flow in the entrance region of rectangular microchannels.
In this study, the characteristics of fluid flow in the hydrodynamic entrance region has been evaluated for a laminar flow in slip flow regime in a two dimensional microscale channel. The single relaxation time approximation of the LBM is used for flow model [34]. As a result, the effect of the Reynolds and Knudsen numbers on the hydrodynamic development length of micro scale channel has been investigated.

NUMERICAL METHOD
The forced flow is moving between two parallel plates as shown in Figure 1. The channel thickness (H) and length (L) are the characteristic length scales as seen on geometry. The thermal properties of the fluid are assumed to be constant. In Maxwell slip model, the boundary velocity depends on Knudsen number; for this reason we used different Knudsen numbers of 0.1, 0.01, 0.001 and repeated these values for the simulations in the slip flow regime.

The Lattice Boltzmann Method (LBM)
The Boltzmann transport equation describes the distribution of particles by means of a distribution function. The equation is related to the total rate of change of the distribution function being equal to the collision rate. Since the collision term is complex, this equation is not easy to solve. One of the approach to this problem is Bhatnagar-Grossi Krook (BGK) model along with simplifying the collision operator [35]. The Equation 1-10 [36] below provide calculation of distribution functions as well as the macroscopic parameters. After the application of this model under the assumption of no external forces exist, Equation (1) is obtained: The equation above is discretized in Lattice Boltzmann Method (LBM) along i th direction. It's a linear partial differential equation where the left and right hand sides represent the streaming and the collision, respectively. Here f represents the distribution function, f eq is the local equilibrium distribution function, c is the unit vector of the particle velocity, and τ is the relaxation factor. The collision frequency shown in Equation (2) is defined as the inverse of the relaxation factor: (2) In LBM, particles can move in prescribed directions in lattice domain that is related to the chosen velocity space discretization. In this study, the most frequent approach for 2D problems, that is D2Q9 model is preferred [36].
In D2Q9 model, the particles can move in horizontal, vertical or crosswise directions to reach the neighbor lattice or they can stay stationary at the center ( Figure 2).
Collision, streaming, calculation of the microscopic values at boundaries and calculation of the macroscopic quantities are categorized in four main steps in a LBM solution loop.
The Equation (1) can be discretized in space and time as below [8]: The collision step is calculated as: (4) In this equation, is the distribution function after the collision that is calculated by using values. The streaming step is calculated as: The equation above is the representation of the arriving of f i (r,t) to the next lattice lattice velocity c i after each time step represented by Δt.
The next step is the update of the microscopic values that is the distribution functions at boundaries. The details of this step will be explained in the following sections. The final step is the calculation of the macroscopic values; velocity components, density and pressure shown in Equation (6,7,8,9).
The local equilibrium distribution function is calculated as: (9) In the equation above, w i is the weighting factor and it is equal to 4/9 for i = 9, 1/9 for i=1,2,3,4 and 1/36 for i=5,6,7,8. c i is the unit vector and c s is speed of sound in lattice units. In D2Q9 approach, c s is . The selection of relaxation factor is important in LBM. The relaxation factor τ and the kinematic viscosity of the fluid must fulfill the relation shown in Equation (10). (10) In this study, τ -1/2 ≪ 1 is considered for the stability reasons and the lattice spacing Δx and the time steps Δt are kept unity since these values are frequently used in LBM studies in open literature.

Boundary Conditions
There are four boundaries in the flow domain to be described in the written code. Uniform inlet velocity and outflow boundary conditions are applied at the inlet and outlet boundaries, respectively. Since the study deals with slip flow regime, special boundary conditions must be applied at the walls instead of the standard ones.
The dimensionless Maxwell first-order slip boundary condition is formulated as, [37]: In the equations above, the accommodation coefficient σ is chosen as 1 by assuming diffuse reflection that is widely used in gas-solid boundaries [38]. u w values are 0 since walls are stationary. Kn is calculated as λ/H. The derivatives can be discretized by using second order scheme. The corresponding top and bottom boundary conditions as well as the inlet and outlet boundary conditions are shown as the following equations.
Bottom Wall Slip Boundary Condition:

Simulation Procedure
A flow chart of the developed solver is simply given below in Figure 3. First, input parameters such as Knudsen number, Re number, the inlet velocity in lattice scale, vertical and horizontal node numbers are defined. Then the grid is generated by dividing the domain into small cells each having a length of unity in both directions. After defining initial conditions, main cycle starts. f values are computed by collision and streaming steps. The application of microscopic conditions proceeds with the calculation of macroscopic parameters. The algorithm goes back to the collision step until the relative variation value of all the parameters with respect to the last time step reaches 10 -5 .

Validation of the Code
After a grid independency study, 510×51 domain is con-sidered to be sufficient to represent the flow for this study and is used for all simulations performed in this paper.
As the lack of experimental result matching suit-ably with this study, the microchannel flow simulation is compared to the results of the work carried out by Karniadakis, et al. [12]. Their flow model seems to be an alternative for the comparison purposes and has been found as highly accurate with the study done by Ahmed and Beskok [39].
In Equation (31), b is a constant that reflects the accuracy of the method and can take values of 0, -1 and -2. The simulation for validation shown in Figure 4 is performed for Re 10 and compared to the analytical formula in Equation 31 where Re effect is not considered. The figure shows that the curve obtained with Equation 31 with b=-1 agrees well with the results obtained in this study by using LBM for micro-flow in the slip regime.
Since the flow domain is modeled as a 2D channel, hydraulic diameter D h is calculated as twice of the height H of the channel. In order to show the results on a non-dimensional scale, y axis is represented as 2y/D h to obtain maximum value at channel upper wall as 1. Horizontal axis is also nondimensionalized by D h .

Developing Velocity Profiles
The non-dimensional velocity profiles with varying Knudsen and Reynolds numbers in the entrance region of a microchannel duct are calculated and shown in Figure  5. The coordinates in the flow direction are non dimensionalised by the hydraulic diameter (x/D h ), whereas the horizontal velocity at the centerline of the channel is normalized by the horizontal velocity at the channel centerline for fully developed condition (U CL/ U CL,FD ). The plots indicate that the velocity profiles change their shape along the channel from uniform to parabola where the minimum and maximum velocities are observed at the channel wall and core as expected, respectively. As the uniform velocity profile enters the channel, it soon after develops, the velocity profile changes due to viscous friction at the wall. As the Kn number increases from 0.001 to 0.1, the velocity near the wall gets a higher value due to increase in the slip velocity. Thus, the centerline velocity (maximum velocity) at fully developed flow condition decreases to satisfy continuity. Hence, the ratio of the inlet velocity to the maximum velocity at the fully developed conditions increases. It is also seen from the plots that as the flow approaches the fully developed region, U CL/ U CL,FD approaches unity since the maximum velocity at the channel center increases along the channel, as expected. For each Re, when Kn increases, the flow regime causes the profile to reach fully developed conditions at a closer position to the outlet of the channel which means that L e /D h is longer. This phenomenon is thought to be directly related to the increasing slip velocity at walls-where slip velocity is proportional to Kn number given in Equation 11 and Equation 12-that decreases the speed of the development of the velocity profile. Moreover, for a constant Kn value, as the Re number increases, the development of the velocity profiles takes longer again which results in a higher value of L e /D h . As Re increases, the ratio of the viscous forces propagated from the walls to the inertial forces decrease. Thus, the penetration of these viscous effects into the channel is slowed down. Consequently, hydrodynamic entrance length increases.
The non-dimensional velocity profiles at the channel exit for various Kn numbers are shown in Figure 6. The coordinates in the flow direction are non dimensionalised by the hydraulic diameter (2y/D h ), whereas the horizontal velocity at the exit of the channel is normalized by the mean velocity at the channel exit for fully developed condition (U FD / Umean,FD ). It is already known that for laminar flow, nondimensional velocity profile for fully developed conditions is independent of Re number.
The results are obtained for all Re numbers in this study and obtained the same profile for the same Kn number. Thus, the results for Re 50 as a representative of other Re numbers in this work is chosen and shown in Figure 6. The effect of the slip is clearly seen on this Figure. It is clear from the velocity profiles that the wall slip increases and the core velocity decreases along with increasing Kn number. The core velocity profile becomes flatter with increasing Kn number. The slip values at walls and the profiles in Kn 0.001 and 0.01 are close together and the slip at walls are small. One should note that Kn 0.001 is just at the limit between slip and continuum regions, so a low slip value is expected. Since the flow regime based on Kn number is given on a logarithmic scale. On the other hand, Kn 0.1 is at the limit of transition flow that shows considerably a higher velocity slip at wall that cause a lower maximum velocity at the center.

Developing Flow Streamlines
In Figure 7 and 8 velocity vectors and contours non-dimensionalized by the inlet velocity are presented with varying Kn numbers at Re=1 and Re=100. The results show that at a given Re number, as the Kn number increases, the velocity magnitude near the wall increases due to increase in slip velocity. The velocity profile uniformity increases and the magnitude of maximum velocity at the center line decreases with the increase of Kn number. The comparison of Kn 0.001 and Kn 0.01 results in Figure 7 show minor difference that is in accordance with the Figure 6. On the other hand, blue regions with low velocity disappears in Kn 0.1 due to high slip velocity at walls. Besides, it is remarked easily that red regions are disappeared as well due to lower maximum velocity at the center and thus a more uniform profile. These comments are valid also for Figure 8. The notable difference compared to Figure 7 is that relatively high Re number causes the profile to develop later that is in accordance with the results shown in Figure 5.

Hydrodynamic Development Length
The hydrodynamic entrance length can be considered as the position along a channel or pipe where the maximum velocity in the related velocity profile reaches 99% of the corresponding fully developed value. The micro-scale channels are generally restricted in size to fit in a given space, especially in electro mechanic systems. Due to its sizes at micro meter scale, the channel entrance effects must be taken into account, particularly in short channels. In some micro channel applications, the flow may not reach to the fully developing regime due to the short length of the  channel. Therefore, the hydrodynamic developing length can be a parameter to examine the channel entrance effects during the design of micro scale channels.
There are several experimental and numerical works carried out for the hydrodynamic entrance length in the literature. Some of those are: Atkinson et al. [40] proposed equation below for the hydrodynamic entrance length: Chen [41] proposed the equation below: (33) Durst et al. [42] studied the development lengths in channel flows both numerically and experimentally. Their proposal for the entrance length: Han [43] proposed below equation for square microchannels: (35) Wiginton and Dalton [44] suggested below equation for an aspect ratio of 1: Renksizbulut and Niazmand [55] obtained a correlation for an aspect ratio between 0.5 <ε < 2 and Reynolds number 10 < Re < 1000: Ahmad and Hassan [46] suggested an empirical correlation for the channel flows to obtain the dimensionless entry length given as equation below: (38) Galvis et al. [47] suggested a correlation to obtain the entrance length for microchannels with aspect ratio between 0.2 < ε < 1 and for Re < 50: (39) The comparison of entrance length data, extracted from this study and the correlations obtained in previous studies above for microchannel flows is showed in Figure 9. As it is seen there is no single correlations but are many different correlations for the entrance length estimations in microchannels. In our present study, a new correlation is obtained for the estimation of the entrance length in rectangular microchannels involving Knudsen numbers between 0.001 to 0.1 and Reynolds numbers 1 to 100.
One can observe from this comparison that, when the Reynolds and Knudsen numbers increase the hydrodynamic entry length (L e /D h ) also increases. At a low Re number, the developing length is also small and the flow is fully developed through a long part of the channel. Based on the present results, the below relationship developed previously by Ma et al. [48] for the entrance length calculation in rectangular microchannels is considered here: (40) where a, b, c are constants, s(Kn) and f(ε) are for describing the relation with L e /D h . As it is a 2D study then aspect ratio relation in this case f(ε)=1.
For the curve fitting, the values of a, b, c are obtained from below equation: Finally, for this study, the coefficients are obtained to give the equation below: (42)

CONCLUSION
In this study, slip flow in micro-scale channel has been investigated by developing a code in LBM. The results of numerical simulations have been used for examining the entrance length and its dependence on Kn and Re numbers in a two dimensional micro scale channel. The fluid velocity profiles, contours and hydrodynamic development length are evaluated by changing Kn within the slip flow regime limits and Re between 1 and 100. Hydraulic development length has been compared with the literature and comparable results are obtained; L e /D h is found to be increasing with the increase of Re and Kn. For each Re, when Kn increases, the flow regime causes the profile to reach fully developed conditions at a closer position to the outlet of the channel that means a longer L e /D h value. This is related to the slip velocity at walls to be proportional to Kn number. When Kn increases, increasing slip velocity decreases the speed of the development of the velocity profile. Moreover, for a constant Kn value, as the Re number increases, the penetration of these viscous effects into the channel is slowed down and the development of the velocity profiles takes longer again which results in a higher value of L e /D h . The slip values at walls and the profiles in Kn 0.001 and 0.01 are found close together and the slip at walls are found small as expected due to relatively lower Kn values. Kn 0.1 which is at the limit of transition flow shows considerably a higher velocity slip at wall. This causes a lower maximum velocity at the center and the velocity profiles to become flatter at the channel exit. A formula for hydrodynamic development length as a function of Kn and Re numbers is obtained by using data extracted from LBM simulations performed. One can conclude that LBM is a relatively easy, fast, simple and promising method as an alternative to conventional modelling techniques and can be powerful to resolve complex problems where applicable.

ACKNOWLEDGMENTS
The authors would like to thank TÜBİTAK Rail Transport Technologies Institute for the infrastructure used for this work.

AUTHORSHIP CONTRIBUTIONS
Authors equally contributed to this work.

DATA AVAILABILITY STATEMENT
The authors confirm that the data that supports the findings of this study are available within the article. Raw data that support the finding of this study are available from the corresponding author, upon reasonable request.

CONFLICT OF INTEREST
The author declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

ETHICS
There are no ethical issues with the publication of this manuscript.