1 Introduction

Laminar-to-turbulent flow transition in microchannels (MCs) is useful to enhance mixing and heat transfer in microsystems. Due to the large pressure drop linked to the small dimensions of MCs, microdevices usually work in laminar regime, especially in the presence of liquid working fluids. On the contrary, due to lower density, Reynolds numbers larger than 2000 can be easily reached for gas microflows with an acceptable pressure drop in channels having hydraulic diameters less than 0.5 mm. Since MCs are used as elementary building blocks of micro-heat exchangers and micro-heat-sinks, it is essential to ascertain under which conditions, the laminar-to-turbulent flow transition inside such geometries can be expected. After the pioneering work of Tuckerman and Pease (1981), most of the earlier experimental groups involved in the analysis of internal flows within microchannels reported an anticipated laminar-to-turbulent flow transition in MCs using liquid (Peiyi and Little 1983; Peng et al. 1994) as well as gas flows (Choi 1991), whereas a few groups like (Harms et al. 1999; Adams et al. 1999) confirmed the agreement to the typical rules of the laminar-to-turbulent transition in conventional tubes using liquid flows. Recent reviews like those of (Asadi et al. 2014; Dixit and Ghosh 2015) seem to give a final answer to the question by highlighting how experimental data obtained in the recent decade confirm that single-phase flows in MCs follow conventional laws. This conclusion is in agreement with the observation of Morini (2004b), according to which macro-theory holds even for MCs if scaling effects can be considered negligible. Morini and Yang (2013) individuated the main scaling effects that need to be verified before interpreting experimental results which can be responsible for deviation from conventional correlations. About the rules of the laminar-to-turbulent transition in MCs, Morini (2004a) gathered previously published experimental data for flow transition in MCs and compared it with the predictions obtained by using the results validated for conventional channels. The analysis on rectangular MCs showed that a series of published values of critical Reynolds number (\({\text {Re}}_\mathrm{c}\)) can be explained by using the Obot–Jones model. Obot (1988) demonstrated that the \({\text {Re}}_\mathrm{c}\) is a function of the cross-sectional geometry of the channel and, for rectangular channels, it tends to increase if the channel aspect ratio (\(\alpha\):height to width) is decreased (shallow channels).

To address the issue of the anticipated laminar-to-turbulent transition in MCs observed by pioneering groups, microparticle image velocimetry (\(\mu\)PIV) studies were conducted by Li and Olsen (2006a) using water as working fluid. They tested four MCs with almost the same hydraulic diameter (\(D_\mathrm{h}\)) and having \(\alpha\) values between 0.2 and 1. They found that \({\text {Re}}_\mathrm{c}\) increased from 1715 to 2315 by decreasing \(\alpha\) from 1 to 0.2. Similar experimentation was also performed by Wibel and Ehrhard (2009) who reported that an increase in \({\text {Re}}_\mathrm{c}\) from 1600 to 2200 is observed for a decrease in \(\alpha\) from 1 to 0.2. Li and Olsen (2006a) and Wibel and Ehrhard (2009) tested a set of rectangular MCs having a geometry which enabled the investigation of the effect of aspect ratio on the laminar-to-turbulent transition by keeping a constant hydraulic diameter. Both studies showed a good agreement with macro-scale predictions for liquid flows. More recently (Kim 2016) presented experimental results of liquid flow and heat transfer for 10 chips containing 10 parallel MCs each. The critical Reynolds number increased from 1700 to 2400 with a decrease in the aspect ratio from 1.0 to 0.25. MC dimensions were changed such that \(\alpha\) and \(D_\mathrm{h}\) varied simultaneously.

Due to the large difference in terms of density between working fluid and seeding particles, optical experimental techniques like \(\mu\)PIV are not recommended for gas flows in MCs and therefore a pressure drop study (Moody chart) remains one of the most efficient tools to observe turbulent transition in microflows. Such analysis for gas flows in 11 different microtubes with \(D_\mathrm{h}=125{-}180\ \upmu {\text {m}}\) has been presented by Morini et al. (2007, 2009). By comparing their results with previously published gas flow pressure drop studies, it was shown that if friction factor is calculated by taking into account compressibility effects and minor losses, \({\text {Re}}_\mathrm{c}\) was in between 1800 and 2000, in complete disagreement with the earlier transitions reported by other groups (Peiyi and Little 1983; Kandlikar et al. 2005). An analytical study, based on an energy gradient that deals with the current subject, is reported by Wang et al. (2012). Based on their analysis, it was shown that in a rectangular channel the shorter side dictates the transition. Moreover, a correlation was also provided to predict \({\text {Re}}_\mathrm{c}\) as a function of \(\alpha\) for rectangular MCs.

While reporting the transitional characteristics of microflows, much less attention has been paid to the inlet characteristics of the MCs. It has been shown by Ghajar and Tam (1994) and Tam and Ghajar (1997) that for conventional sized pipes, inlet shape has a significant effect on the transition and hence heat transfer. Similar study was repeated more than 20 years later for rectangular MCs by Dirker et al. (2014) where they employed three different inlet manifold shapes namely sudden contraction, bellmouth and swirl inlet. Results showed that sudden contraction shows the most delayed transition at around Reynolds number (Re) of 2000 whereas swirl inlet type showed the transition as early as Re of 1500. Studies by both of the above-mentioned groups employed liquid as a working fluid. Although effects of different inlet shapes using isothermal air flow in a conventional sized duct have recently been reported by Moruz et al. (2019) no such analyses however, experimental or numerical, are up to now available for gas flow inside MCs. In order to fill this gap, in this paper a detailed analysis of the laminar-to-turbulent transition in rectangular microchannels in the presence of gas flows will be presented. The experimental values of the \({\text {Re}}_\mathrm{c}\) are compared with the predictions obtained by a RANS-based transitional turbulence model, as proposed by Menter et al. (2006) and later modified by Abraham et al. (2008). This is done to assess the applicability of such a model for the prediction of gas flow behavior in transitional regime in MCs. This model is further used to complement experimental investigations on the effect of the aspect ratio of MCs on the laminar-to-turbulent flow transition.

2 Experimental setup and data reduction

2.1 Testing apparatus

The experimental test bench and 3D model of MC assembly used in current work are shown in Fig. 1. Pressurized nitrogen gas is allowed to enter MC assembly perpendicular to the axial direction of MC and leaves perpendicularly as well through the outlet manifold. The total pressure drop along the length of MC assembly is measured using a differential pressure transducer (Validyne®  DP15). The differential pressure transducer has an interchangeable sensing element which allows for accurate measurements over the whole range of encountered pressure drop. Atmospheric pressure is measured using an absolute pressure sensor (Validyne®  AP42). A K-type, calibrated thermocouple is used to measure the temperature at the entrance of the inlet manifold. Thermocouple voltage and the amplified voltage signals from pressure sensors are fed to the internal multiplexer board of Agilent 39470A which is used as a switch. Voltages from the switch are measured using Agilent 34420A with a \(7-1/2\) digits resolution and are finally stored in a PC using a Labview® program. Sensors’ uncertainties are given in Table 2.

Fig. 1
figure 1

Experimental setup (a) and 3D model of MC assembly (b)

Full size image

2.2 Channels fabrication and characterization

MCs were fabricated by milling a PMMA plastic sheet (5 mm thick) using Roland® MDX40A. Modeling and the subsequent CNC toolpath generation were carried out in Autodesk® Fusion360. A flat end mill with two flutes was used to make MC slot with spindle rotational and translational speed of 10,000 and 300 mm/min, respectively. A constant depth of cut of 100 \(\upmu\)m was used. As shown in Fig. 1b, the lower PMMA sheet containing milled MC was sealed using a top cover of the same material with an O-ring in between. Both top and bottom plastics were compressed, between two 5-mm-thick aluminum plates to ensure leak tightness, using 16 peripheral M6 bolts. The test section was checked for leakage by applying a pressure of almost 10 bars at the inlet connector with outlet connector being closed. Later, the inlet connector was also closed and pressure inside the assembly was then monitored for at least 5 h to spot any major leakage. The test section was considered leakage-free if a decrease rate of the measured pressure was observed to be less than 3% per hour. Such a test was repeated before initiating each experimental campaign. Local static pressure was measured at 3 axial locations (\(x= 0.58\,{\text {L}}, 0.72\,{\text {L}}\) and 0.87L) using small holes (200–250 \(\upmu\)m) drilled in the top plastic. A single differential pressure transducer was used for local pressure readings thanks to a solenoid valve-switching assembly. The differential pressure between each pressure port and the atmosphere was measured. Moreover, the absolute atmospheric pressure is also continuously recorded. The absolute pressure of each cross section is simply obtained by adding the measured absolute atmospheric pressure and the measured differential pressure.

A total of 11 MCs were manufactured for the experimental analysis. At first, four MCs microfabricated with different inlet manifold shapes, namely sudden contraction (SC), V-shaped (VS), rounded entrance (RE) and bellmouth (BM) were used to study the effect of inlet shape on the laminar-to-turbulent flow transition. In Fig. 2, the PMMA chips containing the four investigated MCs with different inlet shapes are shown. On the contrary, a set of 7 MCs (MC1–MC7; see Table 1), having a RE inlet manifold (see Fig. 1b), were fabricated with the aim to study the effect of aspect ratio (\(\alpha =h/w\)) of the MC on the laminar-to-turbulent flow transition. Dimensions of MCs were determined using two methodologies. At first, the realized MCs were observed under an inverted microscope (Nikon Eclipse TE 2000-U). Due to the open channel, it was possible to mark the top surface of the chip and then focusing it using the microscope. The depth of channel cross section was then found out by measuring the total distance traveled by the lens (in \(1~\upmu\)m steps) to focus the bottom face of the channel that could be observed with ease due to milling signature of the blade as shown in Fig. 3. On the other hand, the width of the MC was measured using pixel counts for images of bottom and top surfaces. The length of one pixel was ascertained by comparing a known fixed length (width of a commercial etched MC) against the total number of pixels in the image. The width and depth of the MC were finally measured using this methodology at five locations along the length of the MC. These values were finally averaged to find out the mean width and height of the realized MC to calculate the average hydraulic diameter.

Fig. 2
figure 2

MCs with different inlet shapes: schematics (a) and fabricated (b)

Full size image
Fig. 3
figure 3

Microscopic view of the top (a) and bottom (b) surfaces of fabricated MC. \(W_\mathrm{t}\) and \(W_\mathrm{b}\) denote widths of top and bottom surfaces, respectively

Full size image

To validate the dimensions of the MC evaluated using a microscope, 7 MCs (MC1–MC7; see Table 1) were also observed using optical profiler (Talysurf CCI). Another reason to use optical profilometer was to estimate the average surface roughness of the MC bottom wall with the selected machining parameters. Surface roughness, width and height of the channels were measured similarly at five axial locations and an average of three readings was taken at each location. A sample of images and extracted data at one specific cross section for a representative MC is shown in Fig. 4. Finally, representative values for each individual location were averaged to find out the mean height and width. The maximum difference between the calculated \(D_\mathrm{h}\) using the microscope and profilometer for all the seven microchannels was less than \(10\%\). Dimensions of the other 4 MCs (with different inlet shapes) are measured using the inverted microscope. In this case, only the surface roughness measurements are not performed and rather due to the same machining parameters selected for the milling process, its value is considered equal to the average value of the measured roughness associated with the other 7 MCs. The measured average width (w), height (h), length (L) and other derived geometrical characteristics of all the MCs are summarized in Table 1.

Fig. 4
figure 4

3D reconstruction of the MC (a), average dimensions (b) and surface roughness (c) of MC4 using optical profilometer

Full size image
Table 1 Geometrical characteristics of the manufactured MCs (Roughness values indicated for channels SC, VS, RE and BM are not directly measured.)
Full size table

2.3 Data reduction

For a unidirectional flow of ideal gas, average Fanning friction factor between inlet ‘in’ and outlet ‘out’ of a MC with length L can be obtained by the following expression for adiabatic compressible flow, as indicated by Hong et al. (2019):

$$\begin{aligned} f_{f,av}= & {} \frac{D_\mathrm{h}}{L}\left[ \left( -2\ln \frac{p_\mathrm{in}}{p_\mathrm{out}}+2\ln \frac{T_\mathrm{in}}{T_\mathrm{out}}\right) \right. \nonumber \\&\left. -\left( \frac{1}{\dot{G}^2R\left( T_\mathrm{in}+\frac{{u_\mathrm{in}}^2}{2c_p}\right) }\right) \right. \nonumber \\&\times \left( \frac{{p_\mathrm{out}}^2-{p_\mathrm{in}}^2}{2}+\frac{B^2}{2}\ln \frac{p_\mathrm{out}+\sqrt{p_\mathrm{out}^2+B^2}}{p_\mathrm{in}+\sqrt{p_\mathrm{in}^2+B^2}}\right. \nonumber \\&\left. \left. +\frac{1}{2}\left( p_\mathrm{out}\sqrt{p_\mathrm{out}^2+B^2}-p_\mathrm{in}\sqrt{p_\mathrm{in}^2+B^2}\right) \right) \right], \end{aligned}$$
(1)

where \(B^2=4\beta \times \frac{{\dot{G}}^2R^2}{2c_p}\times \left( T_\mathrm{in}+\frac{u_\mathrm{in}^2}{2c_p}\right)\) and \(\beta\), the kinetic energy recovery coefficient, is taken as 2 for laminar and 1 for turbulent flows. p and T denote cross-sectional average pressure and temperature of gas and \({\dot{G}}\) is mass flux (\({\dot{G}}=\frac{{\dot{m}}}{A}\)). R and \(c_p\) denote the specific gas constant and heat capacity of the fluid at constant pressure, respectively. \({\dot{m}}\) represents the mass flow rate, A is the cross-sectional area of the MC and \(u_\mathrm{in}\) denotes the average fluid velocity at the inlet. Hydraulic diameter of a rectangular MC is defined as:

$$\begin{aligned} D_\mathrm{h}=\frac{2wh}{w+h} \end{aligned}$$
(2)

Using the measured mass flux (\({\dot{G}}\)) and evaluated viscosity (\(\mu\)) at inlet temperature, Reynolds number is calculated with the following equation:

$$\begin{aligned} {\text {Re}}=\frac{{\dot{G}}D_\mathrm{h}}{\mu }. \end{aligned}$$
(3)

In addition, under the hypothesis of adiabatic compressible flow, the energy balance for one dimensional Fanno flow, between inlet ‘in’ and any other cross section at a distance ‘x’ from the inlet of the MC yields to the following quadratic equation for the estimation of average cross-sectional temperature, as demonstrated by Kawashima and Asako (2014):

$$\begin{aligned} \beta \left( \frac{\rho _\mathrm{in} ^2 u_\mathrm{in} ^2R^2}{2c_p p_x ^2}\right) T_x ^2 + T_x - \left( T_\mathrm{in} + \frac{u_\mathrm{in} ^2}{2c_p}\right) =0 \end{aligned}$$
(4)

Finally, knowing the average pressure and temperature of a specific cross section, average density and velocity of compressible gas can be obtained using ideal gas and continuity equations, respectively. An initial estimate of the gas velocity at MC inlet is made using the measured mass flow rate and density of the gas at the entrance of inlet reducer/connector (see Fig. 1b). MC inlet properties are then calculated iteratively using the methodology outlined in Rehman et al. (2019); Hong et al. (2016). Once inlet conditions are determined, the temperature at MC outlet can be easily evaluated using Eq. 4, given the cross-sectional pressure is known. In current work, a hypothesis of perfect expansion is considered which essentially assumes atmospheric pressure at the MC outlet (\(p_\mathrm{out}=p_\mathrm{atm}\)). This hypothesis holds if the gas flow is far from chocking condition. Using numerical simulations, Hong et al. (2018) showed that the choked flow is expected in the laminar and transitional flow regimes for microtubes having \(D_\mathrm{h}<100~\upmu\)m. Since the smallest \(D_\mathrm{h}\) among the currently experimented channels is 220 \(\upmu\)m, perfect expansion of the gas can safely be assumed in transitional and early turbulent regimes for all the MCs during experimental data reduction. Moreover, during the experimental runs the maximum Mach number reached in the transitional regime at the last pressure port (\(x/L=0.87\)) for the MC with the smallest \(D_\mathrm{h}\) is 0.37 which is way below the sonic or supersonic choking expected in MC gas flows (Kawashima et al. 2016; Hong et al. 2018; Rehman et al. 2019). Using the outlet temperature evaluated by Eq. 4 and under perfect expansion assumption, the average Fanning friction factor of MC can finally be calculated using Eq. 1. In this analysis, not only the average value of the friction factor along the whole length of MC is obtained by the experimental tests but also the semi-local friction factor values are evaluated using the pressure taps created along the length of MCs. The semi-local friction factor (\(f_\mathrm{SL}\)) is defined as the average value of the friction factor between two local pressure taps. To calculate \(f_\mathrm{SL}\), a simplified formulation of Eq. 1 is used between two adjacent ports, ‘a’ and ‘b’ as given by Kawashima and Asako (2014):

$$\begin{aligned} f_\mathrm{SL}=\frac{D_\mathrm{h}}{L}\left[ {\frac{p_{a} ^2 -p_{b} ^2 }{RT_\mathrm{av} {\dot{G}}^2}-2\ln \left( \frac{p_{a}}{p_{b}}\right) +2\ln \left( \frac{T_{a}}{T_{b}}\right) }\right], \end{aligned}$$
(5)

where \(T_\mathrm{av}\) is the average temperature of the gas between the pressure ports ‘a’ and ‘b’. In order to assess the accuracy of the experimental measurements, the uncertainty associated with each instrument used in the experimental campaigns is reported in Table 2. In a detailed study, Yang (2012), based on the theory of error propagation, highlighted that uncertainty on the evaluated friction factor is highly influenced by the uncertainty on hydraulic diameter. In the present case, the uncertainty on \(D_\mathrm{h}\), evaluated by using the optical profilometer, is of the order of \(1.5 \%\). The uncertainty associated with the measurement of L is around \(0.4\%\). Therefore, the uncertainty on \(f_\mathrm{f}\) (with Eq. 1) and Re (with Eq. 3) is estimated to be equal to \(\pm 10\%\) and \(\pm 3\%\), respectively. All these uncertainties are summarized in Table 3.

Table 2 Typical uncertainties of instruments used
Full size table
Table 3 Typical uncertainties of experimentally deduced parameters based on theory of propagation
Full size table

3 Numerical modeling

From Table 1, it is possible to observe how the channels MC1–MC7 are characterized by different aspect ratio values and hydraulic diameters. To individuate the effects of these two parameters (\(\alpha\) and \(D_\mathrm{h}\)) separately on the flow transition, a numerical parametric study was performed in parallel to the experimental tests. This was done using a commercial CFD code (CFX-ANSYS 18.1) in which \(\gamma -{\text {Re}}_{\theta }\) SST transitional turbulence model was adopted for the analysis of transition flows. This intermittency-based model due to Menter et al. (2006) was originally proposed for the analysis of laminar-to-turbulent flow transition in external flows. Later, Abraham et al. (2008) proposed modifications of two model coefficients and showed that Menter’s intermittency model can also be proficiently employed to predict the laminar-to-turbulent transition in internal flows. In this work, the modified formulation of Menter’s model, as proposed by Abraham et al. (2008), was used.

For the numerical analysis of the effect of inlet geometry shape on the laminar-to-turbulent flow transition, exact replicas of the tested manifold shapes (i.e., SC, RE, VS, BM) were modeled. The computational model accounted for all the geometrical features of the experimental set-up as shown in Fig. 5. The computational domain included the microchannel, inlet (of respective shapes) and outlet manifolds, and two inlet/outlet reducers. Inlet reducer connected the test section assembly to the piping from the mass flow meter whereas outlet reducer vented the gas exiting from MC to the atmosphere. In these set of simulations, importance has been given to the joining part of the manifold and MC where the mesh is significantly refined. A structured mesh of \(40\times 40\) was used for cross section of MC and manifolds whereas the number of divisions along the lengths of MC and manifolds were equal to 100 and 40 each (both inlet and outlet), respectively. A dense mesh in the cross section was selected in order to guarantee that the non-dimensional geometric distance of the first node element from the MC wall (\(y^+\)) was lower than 1 for all the simulated cases. For all the simulated cases, mesh orthogonality stayed between 0.95 and 1. Height of the reducer (\(H_\mathrm{in}\), see Fig. 5) was 30 mm with an internal diameter of 4 mm. The height of the inlet manifold was kept the same as the height of the simulated MC (h). Irrespective of the inlet shape of the MCs, all the simulated MCs had a sudden expansion outlet where gas, after passing through the outlet reducer, was vented to the atmosphere. Two identical reducers are present at the inlet and outlet of the assembly; they were modeled in the computational domain as they might influence the total pressure drop, and hence friction factor evaluation, due to the presence of 90\(^{\circ }\) flow bends. Boundary conditions were placed at the inlet/outlet of the reducers.

Fig. 5
figure 5

Computational model used for numerical analysis of bellmouth (BM)-shaped MC

Full size image

For the analysis of the effect of aspect ratio on the flow transition, a simpler computational model was used instead where inlet and outlet manifolds were not modeled. However, in order to capture any underexpansion for highly turbulent flow, a diverging nozzle was modeled at the outlet of the MC. To complement the experimental campaign by extending the range of aspect ratios investigated, in current numerical simulations \(\alpha\) was varied from 0.1 to 1 by keeping a constant value of \(D_\mathrm{h}\). Three \(D_\mathrm{h}\) values (\(100\, \upmu {\text {m}}, \ 300 \,\upmu {\text {m}}, \ 500\, \upmu {\text {m}}\)) were considered and for each fixed value of \(D_\mathrm{h}\), four \(\alpha\) values (0.1, 0.25, 0.5 and 1) were simulated. A total of 12 MCs were considered; the geometrical dimensions of the rectangular MCs modeled numerically with the aim to complement the experimental campaign are summarized in Table 4. The length of the MC was kept constant in all the cases (\(L=100\) mm), which provided a ratio \(\frac{L}{D_\mathrm{h}} \ge 200\) even for the largest \(D_\mathrm{h}\) considered. This ensured that there were no entrance effects present in the friction factors evaluation as highlighted by Cheng et al. (2017) and Barlak et al. (2011). A similar meshing strategy, as described earlier for the MCs with different inlet shapes, was adopted. As the total length of the MC in this case was greater than earlier explained computational models (ones with different inlet shapes), the number of elements along the length in these simulations was increased to 200.

Table 4 Geometric dimensions of simulated MCs for the analysis of aspect ratio
Full size table

3.1 Solution procedure

In order to obtain the desired Re, the mass flow rate calculated using Eq. (3) was imposed at the inlet of MC. Ideal nitrogen gas entered the reducer/manifold and was vented to the atmosphere using an orthogonal outlet reducer. The value of turbulence intensity (TI) specified at the inlet was 5% for all considered Re, the same value adopted by Abraham et al. (2008). In the current work, the inlet boundary condition was set at the entrance of the domain (reducer entrance) for the cases where the inlet manifold was modeled (VS, RE, BM, SC). In all the other cases (MC8-MC10), it was applied directly at the inlet of MC, in agreement with the works due to Abraham et al. (2008) and Minkowycz et al. (2009).

Steady-state RANS simulations were performed for all turbulent cases. High-resolution turbulence numerics were employed with a higher-order advection scheme available in ANSYS-CFX (ANSYS 2018). Pseudo-time marching was done using a physical timestep of 0.01 s. A convergence criteria of \(10^{-6}\) for RMS residuals of governing equations was chosen while monitor points for pressure and velocity at the MC inlet and outlet were also observed during successive iterations. In case where residuals stayed higher than convergence criteria, the solution was deemed converged if monitor points did not show any variation for 200 consecutive iterations. Reference pressure of 101 kPa was used for the simulation and all the other pressures were defined with respect to this reference pressure. Due to small measured surface roughness (around 1 \(\upmu\)m), walls of the MC were treated as smooth in numerical model. A no slip boundary condition was applied at the walls. Energy equation was activated using total energy option available in CFX which adopts energy equation without any simplifications in governing equations solution. Kinematic viscosity dependence on gas temperature was taken into account using the Sutherland’s law.

$$\begin{aligned} \mu = \mu _{\mathrm{ref}} \ \left( \frac{T}{T_{\mathrm{ref}}}\right) ^{\frac{3}{2}}\left( \frac{T_{\mathrm{ref}}+S}{T+S}\right) , \end{aligned}$$
(6)

where constants for nitrogen gas are:

$$\begin{aligned} \mu _{\mathrm{ref}}= & {} {1.7812\times 10^{-5}}\,{\mathrm{Pa}\ \mathrm{s}} \\ T_{\mathrm{ref}}= & {} {298.15}\,{\mathrm{K}}\\ S= & {} {111}\,{\mathrm{K}} \end{aligned}$$

Further details of boundary conditions are summarized in Table 5. To estimate the average friction factor from numerical results two virtual cross sections, defined at 0.05% and 99.95% of the MC length, were treated as the inlet and outlet of MC, respectively. The results obtained in correspondence of these cross sections were further analyzed in MATLAB to evaluate numerical friction factors by using Eq. 5. It is to be noted that in the numerical model atmospheric pressure was specified at the outlet of the whole computational domain and not at MC outlet, as it was the assumed in experimental data reduction due to unavailability of the measured pressure at the MC outlet.

Table 5 Boundary conditions
Full size table

3.2 Grid independence

A mesh independence study was performed for channel 9d (\(D_\mathrm{h}={300}\,{\upmu {\text {m}}}\) and \(\alpha =1\)). Four different meshes were utilized where numbers of elements were systematically increased as shown in Table 6. For each mesh, three values of Re (1000, 2400, 10,000) were simulated. In the laminar and fully turbulent flow regimes, the difference in terms of average friction factor (\(f_\mathrm{f}\)) between different mesh setups was insignificant whereas it changed considerably in the transitional flow regime. Relative error in average \(f_\mathrm{f}\) at \({\text {Re}} =2400\) was 25% between Coarse and Medium and between Medium and Fine mesh. When this comparison was made between Fine and Very Fine mesh, however, the difference in terms of \(f_\mathrm{f}\) was reduced to only 2.5%. Therefore, considering a trade-off between computational time and solution accuracy, a fine mesh with 320,000 elements was used for further analysis. All the MCs were meshed with similar strategy by putting an emphasis on the near wall region such that \(y^+\le 1\) and to have at least 10 number of elements in the boundary layer for all considered cases.

Table 6 Types of mesh for grid independence study generated for MC9d
Full size table

4 Results and discussion

Laminar-to-turbulent flow transition is established using the friction factor trend as a function of the Reynolds number. Critical Reynolds number (\({\text {Re}}_{\mathrm{c}}\)) is defined as the Reynolds number in correspondence of which the friction factor attains its first minimum and then starts to increase. Following this definition, \({\text {Re}}_\mathrm{c}\) becomes dependent on the experimental conditions and thus it is expressed by means of a range of Reynolds numbers and not as an absolute value. Such range is individuated during the post-processing of the experimental and CFD results by using MATLAB. Theoretical values of friction factor in the laminar regime for the tested rectangular MCs are obtained by means of the Shah and London correlation (Shah and London 1978; Sahar et al. 2017):

$$\begin{aligned}&f_{\mathrm{R},{\mathrm{lam}}}= \frac{96}{{\text {Re}}}\left( 1- 1.3553\alpha \right. \nonumber \\&\left. \quad + 1.9467{\alpha }^2 - 1.7012{\alpha }^3 + 0.9564{\alpha }^4 - 0.2537{\alpha }^5\right) \end{aligned}$$
(7)

For turbulent regime, the following Blasius law can be used to compare the frictional factor results from experimental and numerical analysis:

$$\begin{aligned} f_\mathrm{B}= 0.3164{\text {Re}}^{-0.25}. \end{aligned}$$
(8)

4.1 Effect of inlet manifold

Due to their reduced axial lengths the four MCs, used for the analysis of the effect of the inlet manifold shape on the laminar-to-turbulent transition (SC, VS, RE, BM), do not have pressure taps along their length. For this reason, only average values of \(f_\mathrm{f}\) are evaluated both experimentally and numerically in this case. By simply observing the trend of the friction factor as a function of the Reynolds number shown by Fig. 6a, it is well evident that the geometry of the inlet is able to change significantly the extension of the transition region between laminar and turbulent regimes. To highlight the difference of \({\text {Re}}_\mathrm{c}\) among various inlet configurations, a comparison of experimental \(f_\mathrm{f}\) is given in Fig. 6a. It can be clearly seen that the BM inlet shows the most delayed transition (\({\text {Re}}_\mathrm{c}=~3470\pm 104\)) compared to the other three MCs. This result is in agreement with the evidences of Ghajar and Tam (1994) and Tam and Ghajar (1997) on circular tubes of 15.8 mm and 19.1 mm using water as working fluid. Tam and Ghajar (1997) explained this result by observing that a smoother and gradual reduction of cross section at the inlet of the channel is able to maintain the laminar state longer than other configurations damping turbulent perturbations at MC inlet. Also Minkowycz et al. (2009) showed, by using \(\gamma -{\text {Re}}_{\theta }\) transitional turbulence model, that if the parabolic velocity profile is maintained at the channel inlet with \(TI=1\%\), an abrupt flow transition from the laminar to turbulent happens at very high Re. More recent results obtained by Dirker et al. (2014) seem to contradict these findings where the authors observed an earlier transition in the presence of a BM inlet (\({\text {Re}}_\mathrm{c}=1800\)) for a minichannel having a hydraulic diameter of 1.05 mm and a delayed transition with a SC inlet (\({\text {Re}}_\mathrm{c}=2000\)). It is to be noted that all these studies dealt with liquid flows whereas more recently (Moruz et al. 2019) used air to analyze the effects of inlet shapes on flow transition in a shallow rectangular channel with \(\alpha\) equals to 0.05. Similar to current results, they also observed that BM inlet shape causes the laminar flow to persist for longer range of Re (\({\text {Re}}_\mathrm{c}=3800\)) compared to other inlet shapes that were experimented. In this case, as it will be demonstrated later, the laminar regime is also sustained longer due to the very low value of the aspect ratio of the channel. The other three inlet geometries (VS, SC, RE), do not show significant difference in terms of \({\text {Re}}_\mathrm{c}\).

In order to better observe the transition, a zoomed view of the transitional region is shown in Fig. 6b. Transition occurs for SC, VS and RE inlet shapes in correspondence of a value of \({\text {Re}}_\mathrm{c}\) equals to 2460 ± 73, 2340 ± 70 and 2090 ± 62, respectively. The numerical results obtained for the MC having a BM inlet are in good accordance with experimental data showing a larger critical Reynolds number with respect to the other inlet geometries with \({\text {Re}}_\mathrm{c}\) being 13% lower than experimentally observed value. In case of SC inlet geometry, numerical results indicate the \({\text {Re}}_\mathrm{c}\) of \(\sim\) 2780 which is again 13% higher than experimental values. For the cases of VS and RE inlet geometries, numerical results do not exhibit significant difference in \({\text {Re}}_\mathrm{c}\) (\(\sim\) 2290 and \(\sim\) 2090). For an easy comparison, all the critical Reynolds numbers obtained in this campaign are reported in Table 7 along with the relevant results appeared in the literature for similar geometries.

Fig. 6
figure 6

Comparison of experimental \(f_\mathrm{f}\) curves for all the tested inlet shapes to establish \({\text {Re}}_\mathrm{c}\) (a) and zoomed view in the range of \(1500< {\text {Re}} < 5000\) (b)

Full size image
Table 7 Numerical and experimental \({\text {Re}}_\mathrm{c}\) for different inlet manifold shapes
Full size table

As underlined in Table 7, for the channel having a RE inlet geometry, two relative minimum points of friction factor are observed before the flow becomes fully turbulent. This singular feature can be observed in Fig. 6b where a zoomed view of the transition region is shown. A mention of such experimental behavior can also be found in Dirker et al. (2014) where for BM inlet shape, authors observed two sudden changes of \(f_\mathrm{f}\) curve before transitioning to fully turbulent flow. In order to establish the extent of transitional region in terms of Reynolds values (\(\Delta {\text {Re}}_{\mathrm{LT}}\)), the difference in Re between the onset of transition (\({\text {Re}}_\mathrm{c}\)) and the onset of fully turbulent regime (i.e., in correspondence of the Reynolds number for which the slope of \(f_\mathrm{f}\) becomes parallel to the Blasius law) (\({\text {Re}}_\mathrm{t}\)), is determined for all the four investigated inlet shapes (see Table 8). From both experimental and numerical results, it can be seen that the length of transitional regime is also higher for BM. This is followed by RE and SC. A discrepancy is observed between numerical and experimental results obtained for VS inlet geometry where experimental results show somewhat abrupt transition whereas numerical results show smooth transition to fully turbulent flow.

Table 8 Numerical and experimental \({\text {Re}}_\mathrm{c}\)
Full size table

4.2 Effect of aspect ratio

4.2.1 Experimental investigations

Experimental results in terms of average friction factors between the inlet and the outlet of the microchannel (\(f_\mathrm{f}\)) are shown in Fig. 7 for MC3 and MC4. Laminar \(f_\mathrm{f}\) is in agreement with the correlation proposed by Hong et al. (2008) (dotted line) for compressible laminar flows at large Reynolds numbers, especially for MC4 having a smaller \(D_\mathrm{h}\). For MC4, average \(f_\mathrm{f}\) is slightly lower than the Blasius correlation. In order to better understand the behavior of the average friction factor, in Fig. 7, the trend of the semi-local value of the friction factor \({f_\mathrm{SL}}\), obtained by using two axial pressure ports close to the microchannel outlet (\(x/L=0.72\) and 0.87), is also shown. It is evident that semi-local \({f_\mathrm{SL}}\) follows well the Blasius correlation in fully turbulent regime, exhibiting larger values with respect to the average friction factor. It is worth mentioning how the experimental average and semi-local \(f_\mathrm{f}\) differ in transitional regime but their trends exhibit the same \({\text {Re}}_{\mathrm{c}}\). This means that the \({\text {Re}}_{\mathrm{c}}\) can be extracted from both the curves (average or semi-local \(f_\mathrm{f}\)). A comparison of the values of the semi-local friction factors measured between the last two pressure ports (\(x/L=0.72--0.87\)) for all the tested MCs is shown in Fig. 8. It is evident that flow transition happens in all the channels when Re is in the range of 1863–2772 in agreement with the classical range mentioned in textbooks for laminar-to-turbulent flow transition in conventional channels (Morini et al. 2007, 2009). For sake of completeness, since experiments were performed at a discrete set of Re, a range of critical Reynolds number is given in Table 9. The results of Table 9 highlight how the transition is delayed in rectangular channels having values of \(\alpha\) lower than 0.5 (i.e., MC4). For channels MC5 and MC7, since \(\alpha\) is almost similar, \({\text {Re}}_{\mathrm{c}}\) is also in the same range. These channels have the smallest \({\text {Re}}_\mathrm{c}\) among all the channels experimented. Last two channels in Table 9 are again showing an almost constant range of transitional Re. The experimental results of Table 9 are obtained with a series of 7 microchannels characterized by different aspect ratio values and diameters. In order to check the dependence of \({\text {Re}}_\mathrm{c}\) on \(\alpha\) and \(D_\mathrm{h}\) separately, a numerical parametric study is performed and results are discussed in the next section.

Fig. 7
figure 7

Average and local friction factor for MC3 (a) and MC4 (b)

Full size image
Fig. 8
figure 8

Comparison of semi-local \(f_\mathrm{f}\) (\(x/L=0.72{-}0.87\)) for all the tested MCs (a) and a zoomed region for Re between 1500 and 4500 (b)

Full size image
Table 9 Experimental values of \({\text {Re}}_\mathrm{c}\) for tested MCs having different aspect ratios
Full size table

4.2.2 Numerical parametric study

In order to complement the experimental tests, a series of numerical simulations has been conducted by varying the gas mass flow rate at the inlet of each MC in such a way that Re is varied between 400 and 20000. For a fixed value of \(D_\mathrm{h}\), the value of \(\alpha\) is varied in order to observe the effect of this single parameter on the laminar-to-turbulent transition. As an example, in Fig. 9, the trend of the average friction factor as a function of the Reynolds number is shown for a series of rectangular MCs having the same hydraulic diameter (100 \(\upmu\)m) and different aspect ratios. It is quite evident that with a decrease in \(\alpha\), flow transition is being delayed for a constant value of \(D_\mathrm{h}\). \({\text {Re}}_\mathrm{c}\) increases from 2000 to 4000 by decreasing \(\alpha\) from 1 to 0.1. This trend is in agreement with the experimental observations (made with \(\mu\)PIV and water as working fluid) of Li and Olsen (2006a) and Wibel and Ehrhard (2006) which noted a decrease in \({\text {Re}}_\mathrm{c}\) for an increase in \(\alpha\) for liquid flows in rectangular MCs. From the results of all the simulated MCs (see Table 4), it can be observed, especially for \(\alpha\) in the range of 0.1–0.25, that laminar breakdown and transition to turbulence is delayed. Similar trends are obtained for MC9 and MC10 which, for the sake of brevity are not reported here. It is worth mentioning here that the deviation of the average friction factor from the Shah and London correlation (see Eq. 7) when Reynolds is larger than 1500 in Fig. 9 must not be confused with an early transition typical of microchannels, as concluded in the past by many authors. This deviation is due to the effects linked to the gas compressibility as deeply discussed by Hong et al. (2008), Hong et al. (2012), Hong et al. (2016) and it is stronger for smaller microchannels.

Fig. 9
figure 9

Average \(f_\mathrm{f}\) curves for MCs with \(D_\mathrm{h} = {100}\,{\upmu m}\) and different aspect ratios (MC8)

Full size image

Numerical data can also be analyzed to check if \(D_\mathrm{h}\) plays a significant role on the onset of transition. In Fig. 10a, the effect of the hydraulic diameter, for a series of rectangular microchannels having the same aspect ratio (\(\alpha =\)1), on the trend of the semi-local friction factor in the transition region is shown. A decrease in \({\text {Re}}_\mathrm{c}\) can be observed in Fig. 10a when \(D_\mathrm{h}\) is decreased. This effect is pertinent to the gas flows due to the compressibility (stronger for smaller microchannels) and it is not present for liquid flows. For this reason, when tested with liquid flow, experimental analysis of Li and Olsen (2006b) yielded no significant influence of \(D_\mathrm{h}\) value (200–640 \(\upmu\)m) on the transition where \({\text {Re}}_\mathrm{c}\) ranged between 1718 and 1885 for square MCs. It is evident from Fig. 10b that the Mach number close to the outlet (\(x/L=0.98\)) of MC in the transitional regime is the highest for the smallest MC (\(D_\mathrm{h}=100~\upmu {\text {m}}\)). At the transitional Reynolds number (\({\text {Re}}_\mathrm{c}\)), Mach number (at \(x/L=0.98\)) for the MC with \(D_\mathrm{h}=100~\upmu {\text {m}}\) is 0.68 and it decreases to 0.4 and 0.28 for \(D_\mathrm{h}\) values of 300 \(\upmu\)m and 500 \(\upmu\)m, respectively. This shows how compressibility promotes the laminar-to-turbulent transition for gas flows through small microchannels having a similar aspect ratio. Also the length of the transition regime is only marginally influenced by the value assumed by the hydraulic diameter if \(D_\mathrm{h}>100~\upmu {\text {m}}\).

Fig. 10
figure 10

MCs with \(\alpha =1\): \(f_\mathrm{SL}\) curve (\(x/L=0.5-0.75\)) (a) and cross-sectional Mach number close to the MC outlet (\(x/L=0.98\)) (b)

Full size image

4.2.3 Comparison with Obot–Jones model

For non-circular ducts, Jones (1976) noted that friction factor monotonically increases with an increase in \(\alpha\) and therefore the hydraulic diameter of non-circular ducts should be corrected with a geometric factor (\({\phi ^*}_\mathrm{G}\)) by using a laminar equivalent hydraulic diameter in order to conserve the hydraulic similitude laws. Obot (1988) later followed the work of Jones and concluded that critical Reynolds number for a non-circular duct of an arbitrary shape ‘G’ (\({\text {Re}}_{\mathrm{cr,G}}\)) can also be associated with that of circular geometry (\({\text {Re}}_{\mathrm{cr,o}}\)) using the geometric factor as follows:

$$\begin{aligned} {\text {Re}}_{\mathrm{cr,G}}=\frac{{\text {Re}}_{\mathrm{cr,o}}}{{\phi ^*}_G}, \end{aligned}$$
(9)

where this geometric function can be written as:

$$\begin{aligned} {{\phi ^*}_G}=\frac{64}{(f {\text {Re}})_{G}}. \end{aligned}$$
(10)

The product of \((fRe)_G\) for a rectangular channel can be obtained by using the Shah and London correlation (Eq. 7) for rectangular ducts. For conventional sized commercial circular tube with a finite roughness (\(\epsilon\)), an onset of transition can be estimated as suggested by Idelchick (1986):

$$\begin{aligned} {\text {Re}}_{\mathrm{cr,o}}=1160\left( \frac{\epsilon }{D_\mathrm{h}}\right) ^{-0.11}. \end{aligned}$$
(11)

By combining Eqs. (9–11), Obot (1988) showed that an estimation of \({\text {Re}}_\mathrm{c}\) can be made by varying \(\alpha\) and \(\epsilon /D_\mathrm{h}\) ratio for conventional rectangular ducts. Morini (2004b) proposed to use the same approach in order to predict the laminar-to-turbulent flow transition in microchannels and called this as Obot–Jones model. This approach has been followed here in order to compare the experimental and numerical results presented in this paper and those available in the open literature with the predictions of the Obot–Jones model (see Fig. 11). It can be seen that the experimental results obtained in this work follow the trendline which represents the predictions of the Obot–Jones model if a value of \(\epsilon /D_\mathrm{h}\) equal to 0.1% is considered except for the two MCs with \(\alpha\) close to 0.7 (MC0 and MC5) which follow a trendline associated with a higher \(\epsilon /D_\mathrm{h}\). Recently, Kim (2016) tested metallic chips each containing 10 parallel MCs of various \(\alpha\). Most of his results with one exception also follow the suggested trend of the Obot–Jones model with \(\epsilon /D_\mathrm{h}\) of almost 1%. Similarly, \({\text {Re}}_\mathrm{c}\) obtained recently by Dirker et al. (2014) also follows the trend exhibited by the results of Kim (2016). The only experimental results obtained with rectangular MCs having \(\alpha \le 0.2\) were obtained by Hong et al. (2012), Hong et al. (2016). These MCs were fabricated in silicon using etching with very small roughness values \(\epsilon\) (2–60 nm). If all the experimental results are compared, current ones seem to suggest values of \({\text {Re}}_\mathrm{c}\) larger than those obtained by the other researcher for a given value of aspect ratio. This might be peculiar to the different test assembly or the material employed to mill the MCs in this work. Except for the works of Hong et al. (2012), Hong et al. (2016), all the other results were obtained with metallic microchannels (mechanically milled) with values of the relative roughness \(\epsilon /D_\mathrm{h}\) of 0.97% for Wibel and Ehrhard (2006), 0.22% for Dirker et al. (2014) and 0.63–2.3% for Kim (2016). On the contrary, in this work milling was performed in PMMA plastic obtaining values of \(\epsilon /D_\mathrm{h}\) between 0.1 and 0.4% which, in general, is lower than the others. This might justify the higher \({\text {Re}}_\mathrm{c}\) observed in the present work. Furthermore, it can also be seen in Fig. 11 that all the presented data follow the Obot–Jones model if \(\epsilon /D_\mathrm{h}\) is between 0.1 and 1%. However, considering the actual values of \(\epsilon /D_\mathrm{h}\) declared by the researchers as well as in the current work, it can be said that Obot–Jones model is able to predict the decrease in \({\text {Re}}_\mathrm{c}\) with increasing \(\alpha\) only qualitatively. In addition, for sake of completeness, it is worth mentioning that Eq. 11 proposed by Idelchick (1986) was suggested for water and not for gases. Although numerical results follow similar trends, they predict a higher \({\text {Re}}_\mathrm{c}\) than experiments for the given \(\alpha\). This can be explained by the fact that in the numerical simulations only the MC is modeled without taking into account manifolds, reducers and other parts of the system that can influence the laminar-to-turbulent transition. In Fig. 11, the analytical correlation presented by Wang et al. (2012) is also plotted for comparison. The correlation shows that \({\text {Re}}_\mathrm{c}\) can go as low as 1220 for \(\alpha =1\). Figure 11 highlights how current findings (both experimental and numerical) and the recent literature are in disagreement with Wang’s correlation.

Fig. 11
figure 11

Comparison of current results of \({\text {Re}}_\mathrm{c}\) and recent literature with Obot–Jones model (Obot 1988; Morini 2004a) for different values of \(\epsilon /D_\mathrm{h}\)

Full size image

Combining all the results presented in Fig. 11, it can be concluded that a clear trend there exists between \(\alpha\) and \({\text {Re}}_\mathrm{c}\) which suggests an increase in \({\text {Re}}_\mathrm{c}\) with a decrease in \(\alpha\). Such an increase becomes significant for lower values of \(\alpha\) (27% in terms of \({\text {Re}}_\mathrm{c}\) when \(\alpha\) goes from 0.5 to 0.1) and is insignificant for higher values of \(\alpha\) (8% when \(\alpha\) goes from 1 to 0.5). Such a small difference is of the same order of magnitude than the typical uncertainty of measured Re in an experimental campaign and for this reason is difficult to observe the difference in \({\text {Re}}_\mathrm{c}\) at higher values of \(\alpha\) (see Table 9). From the numerical results presented in Fig. 11, it can be observed that the difference in terms of \({\text {Re}}_\mathrm{c}\), between MCs having the same \(\alpha\) and different \(D_\mathrm{h}\), diminishes for smaller values of \(\alpha\). This signifies that the effect of \(D_\mathrm{h}\) on the values of \({\text {Re}}_\mathrm{c}\) is only significant for \(\alpha =0.5{-}1\), as discussed in Fig. 10. Nonetheless, \({\text {Re}}_\mathrm{c}\) for all the cases presented in this paper as well as in the open literature is higher than 1700 even for the highest value of \(\alpha\) reported. The decrease in \({\text {Re}}_\mathrm{c}\) is thwarted by compressibility for \(D_\mathrm{h}\le 100~\upmu {\text {m}}\). In fact, compressibility tends to anticipate the laminar-to-turbulent transition significantly for higher values of \(\alpha\). The results presented in this paper confirm that no anomalous earlier transition is generally observed, even for gas flows in MCs, in disagreement with the predictions of the correlation proposed by Wang et al. (2012).

5 Conclusions

Laminar-to-turbulent transition of gas flow in rectangular microchannels is discussed by analyzing the effect of two geometrical parameters, the aspect ratio and the shape of the inlet manifold. At first, the effects of inlet manifold shape on the critical Reynolds number are investigated. Experiments were performed for four different inlet manifold shapes. Results show that the highest value of critical Reynolds number is obtained with an inlet region in which the cross section is gradually and smoothly reduced before the entrance of the channel (Bellmouth geometry). The experimental observations have been complemented by numerical simulations which confirmed the results observed during the experimental runs. A very good match between experimental and numerical results in terms of critical Reynolds number has been observed for all the cases considered.

The effect of the aspect ratio of the rectangular cross section of the microchannel has been studied both experimentally and numerically by means of a series of cases in which the effect due to aspect ratio is separated from the effect of hydraulic diameter of the channel. It has been found that the effect of the aspect ratio on the critical Reynolds number becomes evident for aspect ratios lower than 0.5. A reduction in the value of aspect ratio causes an increase in the critical Reynolds which means a delayed laminar-to-turbulent transition, especially for aspect ratios lower than 0.25.

Finally, the effect of the hydraulic diameter on the laminar-to-turbulent transition was checked for a fixed value of the aspect ratio of the microchannel. It has been observed that the effect of the hydraulic diameter on the flow transition is less significant but in the presence of a gas flow a reduction of the hydraulic diameter determines an early transition due to the compressibility effects.

By comparing the existing experimental and numerical results as well as the recent literature on this subject, it has been shown that critical Reynolds numbers, obtained by testing liquid and gas flows in rectangular microchannels, qualitatively follow the Obot–Jones model. Under the conditions considered in this work (\(D_\mathrm{h}\ge 100~\upmu {\text {m}}\)), the laminar-to-turbulent transition takes place in correspondence of critical Reynolds numbers in the range 1700–3500. In agreement with the conventional theory, no anomalous earlier transition with critical Reynolds numbers less than 1700 is observed.