Hydrodynamic Entrance Length for Laminar Flow in Microchannels with Rectangular Cross Section

This work presents a detailed numerical investigation on the required development length (L = L/B) in laminar Newtonian fluid flow in microchannels with rectangular cross section and different aspect ratios (AR). The advent of new microfluidic technologies shifted the practical Reynolds numbers (Re) to the range of unitary (and even lower) orders of magnitude, i.e., creeping flow conditions. Therefore, accurate estimations of L at Re ≤ O(1) are important for microsystem design. At such low Reynolds numbers, in which inertial forces are less dominant than viscous forces, flow characteristics become necessarily different from those at the macroscale where Re is typically much larger. A judicious choice of mesh refinement and adequate numerical methods allowed obtaining accurate results and a general correlation for estimating L, valid in the ranges 0 ≤ Re ≤ 2000 and 0.1 ≤ AR ≤ 1, thus covering applications in both macro and microfluidics.


Introduction
Over the past two decades, many microfluidic systems have been designed and built for a wide range of applications. These systems are essentially formed by microchannels and chambers whose geometries have a strong influence on the flow characteristics. Improving the design of such systems requires understanding in detail the fluid flow characteristics. One of the flow characteristics that still needs to be better quantified is the microchannel entrance length, required to achieve fully developed flow, which is an important design parameter especially when it represents a significant fraction of the total microchannel length. However, entrance effects are typically ignored by researchers who usually assume the laminar flow to be fully developed [1].
Different cross-sections (e.g., rectangular, trapezoidal, triangular, circular and elliptical) have been used by numerous researchers with the purpose of exploring the flow behavior in microchannels [2][3][4][5][6][7][8][9]. One important fact that should be highlighted is that the flow field in rectangular channels is more complex than in circular channels, because of the additional dependence on the aspect ratio (AR) of the cross-section, as shown in Figure 1.
When a fluid enters a rectangular channel, the initial uniform velocity profile is gradually redistributed, accelerating around the centerplane(s) and de-accelerating near the walls. The fluid will reach a location after which the velocity profile no longer changes in all directions, and under such conditions the flow is considered to be fully developed. Estimating this entrance length, also called development length, is one of the classical problems in fluid mechanics with practical use in microsystems design and studies of transition to turbulence [10]. Here, following the criteria proposed by Shah and London [11], the normalized development length, L = L/B, is defined as the non-dimensional distance measured from the entrance of the microchannel, where the velocity profile is uniform, up to the axial location where the flow reaches 99% of the corresponding fully developed velocity profile. Many analytical, experimental and numerical works on the development length problem have been published over the years [11][12][13][14].
Most of the previous works proposed a linear functional relation between the development length and Reynolds number, L = C 1 + C 2 Re, in which C 1 represents the asymptotic limit in the creeping flow regime, Re → 0. The proposed values for C 1 in the literature are scattered, with the linear correlation hardly predicting the development length accurately for low Reynolds number flows (see Table 1). Table 1. Existing correlations to estimate L = f (Re) for low Re channel flows between parallel plates (Re = ρUL * /η, with L * representing the channel width). num. and exp. stand for numerical and experimental work, respectively.

Remark 1.
It should be noted that, to present a correlation where the results for the different AR can be comparable, the dimension B was used as the characteristic length scale in all dimensionless development lengths. This means that this correlation cannot be compared with the one from Durst et al. [17], for planar channels, AR → 0. We define Re = ρUB/η and Re H = ρUH/η, thus Re H = AR × Re. Note also that the Reynolds based on the hydraulic diameter (D h ) is given by For a square duct, B = H and Re D h = Re.
The nonlinear behaviors of the entrance length with respect to the Reynolds number for pipe and parallel plate flows were first investigated by Atkinson et al. [15], using a linear combination of creeping flow with boundary-layer type solutions, and Chen [16], using the integral momentum method to derive approximate solutions for the development length. The resulting functional relations for L in planar channel flows are listed in Table 1, along with the appropriate fitting constants. These nonlinear correlations indicate that the development length has a rational relationship with the Reynolds number when the Reynolds number is small, with L approaching a constant value as Re → 0.
When the Reynolds number increases, the constant and rational portions become less dominant and the entrance length variation with Re becomes linear.
Sadri and Floryan [21] presented a linear correlation for the relation of L(Re) ranging from Re = 0.01 to 2200. They relied on a numerical method based on the use of the stream function and vorticity in the flow governing equations, with a compact fourth-order finite-difference scheme. The values of their correlation overpredict the entry length at low Re due to the omission of the upstream flow development and at high values of Re due to the omission of flow separation effects (they decoupled the upstream and downstream regions of channel entry). Sadri and Floryan also compared different definitions of the entry length and studied their influence on the obtained development lengths of the flow.
More recently, Durst et al. [17] improved the linear functional prediction accuracy by superimposing diffusion and convection together in their model while suggesting a method to determine the normalized entrance length in laminar pipe and planar channel flows. It should be noted that, in the works previously cited, 3D effects were neglected.
Recent experimental fluid flow analysis, focused in microchannels with rectangular and square cross-sections [18,19], also present nonlinear correlations for L as a function of Re. Their conclusion is that the development lengths were found to be shorter than classical correlations found in [17]. The explanation for this difference is based on the effect of the aspect ratio, the different off-center velocity maxima of the inlet profiles and due to the pre-development of the axial velocity. Recently, Lobo and Chatterjee [20] presented a numerical study on the development of flow in a square mini-channel, taking into account the effect of flow oscillations. They also performed a steady-state study and proposed a correlation for the development length for a square microchannel.
In this work, we present a detailed numerical investigation on the development length (L = L/B) in rectangular microchannels by considering the laminar regime (including the limit of creeping flow conditions) and the dependence on the aspect ratio, AR.
It should be remarked that, although some works dedicate particular attention to the effects of surface roughness, wettability and the presence of gaseous layers (thus multiphase flows) leading to interfacial slip typical of some specific fluids and particular conditions [22], here we assume a single phase fluid and that the quality of manufacturing is sufficiently good to ensure smooth surfaces (as in channels made of PDMS through soft lithography [23] or made from fused-silica glass [24] so that the no-slip boundary condition remains valid).
The remainder of this work is organized as follows. In the next section, we present the governing equations, the numerical method and computational meshes used in the simulations. In Section 3, we present the results and discussion. The paper ends with the main conclusions.

Governing Equations, Numerical Method and Computational Meshes
To predict the developing flow field within the microchannel, we assume that the flow is three-dimensional, laminar, isothermal, incompressible and steady. The governing equations for such conditions are those expressing conservation of mass, and the momentum equation, where u is the velocity vector, p is the pressure, t is the time, ρ is the fluid density and τ is the Newtonian extra stress tensor, defined as: with D representing the symmetric rate of strain tensor and η the dynamic viscosity. An in-house finite-volume method code (FVM) was used in this numerical study (for more details, see [25,26]). In this code, all diffusive terms are discretized using second-order central differences, and the advective term in the momentum equation is discretized with the CUBISTA high-resolution scheme [25], which is formally of third-order accuracy. The term ∇ · τ in Equation (2) is usually represented as η∇ 2 u in Newtonian fluid flows, but our numerical method is also applicable for more complex fluids and uses the stress tensor notation. The channel geometry is represented in Figure 1. A uniform velocity, U, is imposed at the inlet (y/B = 0) of the microchannel, where all other velocity components and deviatoric stresses are set to zero. Vanishing streamwise gradients are applied to all variables at the outlet plane, except for the pressure, which is linearly extrapolated to the outlet from the two nearest upstream cells. The classical no-slip condition is imposed at the channel walls.
The iterative time stepping procedure was stopped whenever the L 2 − norm ∑ n i=1 |x i | 2 of the residuals vector of the system of equations was less than a tolerance of 10 −8 , which corresponds to steady-state flow. The calculations for the creeping flow limiting case, Re → 0, were carried out by dropping out the convective term of the momentum equation.
In order to achieve mesh-independent results, the full physical domain depicted in Figure 1 was discretized into eight different sets of computational domains, with consecutively refined meshes. A detailed summary of those meshes is presented in Table 2. To avoid an effect of the computational domain on the calculated development length, the microchannel was made significantly longer than the value of L, so that, at the end, the length of the microchannel domain varied between four and twenty times the calculated development length. In turn, this depended also on the Reynolds number.  For the x direction, the cells were uniformly distributed with the total number of cells, NC, shown in Table 2. The present mesh refinement strategy was designed in order to provide a consistent estimate of uncertainty of our simulations using Richardson's extrapolation technique [27]. This method allows us to estimate the mesh independent extrapolated development length (L ext ) as L ext = L M4 + (L M4 − L M3 )/(2 p − 1) based on our numerical simulations in meshes M3 and M4, for the sake of example, where p is the order of accuracy of the computations. Assuming a second-order accuracy, p = 2, L ext can be estimated as L ext = L M4 + (L M4 − L M3 )/3. This technique also allows us to estimate the accuracy of a specific calculation, defined as where L i is the length computed numerically in mesh M i . In addition to the meshdependency study, another parameter used for estimating the accuracy in this work is the velocity relative error (e u,M i ) for mesh M i , defined as, where u max,num is the predicted maximum velocity at the outlet plane and u max,theo is the corresponding fully developed theoretical maximum velocity value computed using the analytical solution [28]. It should be remarked that preliminary tests with different meshes where conducted initially, in order to select the most appropriate meshes, without compromising the convergence accuracy and the computational time needed to perform the simulations.

The Development Length
In order to estimate L, simulations were conducted from small Re (creeping flow conditions) up to Re = 2000, covering the laminar regime. Note that for high AR the simulations were performed for lower values of Re, so that the flow stays laminar (see, e.g., [18]).
The computational domain size of the microchannel for the 3D flow field is shown in Table 2. A square section (AR = 1) and three additional rectangular sections with AR = 0.5, AR = 0.25 and AR = 0.1 were investigated.
The results obtained for L ext considering different AR are presented in Table 3, for a number of representative Re values. It should be remarked that the errors are below 3.5%, within the asymptotic range of convergence (see Appendix A, Tables A1-A4 for a  The data shown in Table 3 can be better understood by looking at Figure 2. A decrease of the development length with the decrease of AR is observed for the whole range of Re considered. For creeping flow, we note a more pronounced decrease of L as AR deviates from AR = 1. For low Re, the difference in the development length for the different AR is less significant when compared with the results obtained for high Re. For high Re, convection starts to play an important role, and, to understand this flow behavior, the local velocity profiles must be analyzed. Note that we also compare our results with the experimental data of Lee et al. [18] for AR = 0.4 and 0.38. The experimental results are qualitatively in agreement with our numerical predictions.   [18]. The dashed lines show the fit obtained with the correlation given in Equation (6). We also show a comparison with the correlation proposed in [20], and a good agreement is observed (except for low Re where the correlation proposed by [20] fails). The correlation proposed in [29] also shows good agreement with the correlation proposed in [20] (in [20], several correlations are developed to account for different duct shapes).   Comparing Figure 3a,c, we observe that for the lower aspect ratio the overshoot is more pronounced, possibly due to the restrictive no slip boundary condition of the top and bottom walls that are much closer to the center than the side walls. Figure 3b shows the evolution of the non-dimensional axial velocity component at several axial positions along the z−direction for the case of AR = 0.1 and Re = 0.001 (the results are qualitatively similar for other AR). In this case, the left and right walls are far from the development region, and we have an almost 2D case, where a more parabolic profile is obtained.

Flow Dynamics
The overshoots are generated as a result of the abrupt local fluid deceleration near the channel wall at the inlet section due to the no-slip boundary condition. Since the fluid near the symmetry axis is not accelerated immediately, whereas the fluid near the wall is stationary as soon as it enters the inlet region, to satisfy the continuity equation, velocity overshoots are formed near the wall, close to the inlet section. This behavior was also reported by Durst et al. [17] for Newtonian fluids in 2D channels.
We can also observe that there is a dependence on the development of the axial velocity with Re, with the increase of inertia leading to longer axial development lengths to achieve fully developed flow (see also Figure 2). The non-dimensional axial velocity component presents a local maximum that becomes more pronounced and closer to the channel walls as the Reynolds number increases.
By comparing the evolution of the non-dimensional axial velocity component at several axial positions along the z−direction for the case of AR = 0.1 and Re = {0.001, 100}, we observe that near the entrance the case of Re = 0.001 shows a more pronounced development, while the combination of inertia/diffusion effects requires higher L to stabilize the flow field in the case of Re = 100 (see Figures 3b and 4c).
The evolution of the non-dimensional streamwise velocity along the center of the duct is shown in Figure 5 for creeping flow conditions and Re = 100 for AR = 1. The development of the axial velocity profiles is not linear, following a logistic-type function behavior. Initially, the fluid evolves slowly as a plug (before the the momentum diffusion from the wall reaches the center) and then quickly evolves to its fully developed value. This occurs for both values of Re considered, but, as expected, achieving the fully developed regime takes higher L when convection plays an important role.

Entrance Length Correlation
In order to obtain a nonlinear correlation between L and Re, we performed a fit to the numerical results shown in Table 3. The following correlation, valid in the ranges 0 ≤ Re ≤ 2000 and 0.1 ≤ AR ≤ 1, was obtained: The correlation follows the same structure of the one proposed by Durst et al. [17], but now C 1 , C 2 and C 3 are functions of AR. For the range of AR studied in this work, we obtained the following functions: An overall average relative error of 2.2% and a maximum relative error of 6.0% were obtained, and the good quality of fit is shown in Figure 2 (dashed lines).

Conclusions
An extensive and systematic numerical study was carried out regarding the entrance length, L, for Newtonian fluids in 3D microchannels with rectangular cross section. An adequate choice of meshes, with consistent mesh refinement, and accurate numerical methods allowed the prediction of accurate values of L, and a unified nonlinear correlation was proposed for estimating L, valid in laminar regime for 0 ≤ Re ≤ 2000 and 0.1 ≤ AR ≤ 1. This correlation allows an accurate prediction of entrance length and represents a novel tool for the design and fabrication of efficient microfluidic devices.