Simulations of Flows via CFD in Microchannels for Characterizing Entrance Region and Developing New Correlations for Hydrodynamic Entrance Length

Devices with microchannels are relatively new, and many correlations are not yet developed to design them efficiently. In microchannels, the flow regime is primarily laminar, where entrance length may occupy a significant section of the flow channel. Therefore, several computational fluid dynamic simulations were performed in this research to characterize the developing flow regime. The new correlations of entrance length were developed from a vast number of numerical results obtained from these simulations. A three-dimensional laminar flow for 37 Reynolds numbers (0.1, 0.2, …, 1, 2, …, 10, 20, …, 100, 200, …, 1000), primarily in low regime with water flow through six rectangular microchannels (aspect ratio: 1, 0.75, 0.5, 0.25, 0.2, 0.125), has been modeled, conducting 222 simulations to characterize flow developments and ascertain progressive velocity profile shapes. Examination of the fully developed flow condition was considered using traditional criteria such as velocity and incremental pressure drop number. Additionally, a new criterion was presented based on fRe. Numerical results from the present simulations were validated by comparing the fully developed velocity profile, friction factor, and hydrodynamic entrance length for Re > 100 in rectangular channels, for which accurate data are available in the literature. There is a need for hydrodynamic entrance length correlations in a low Reynolds number regime (Re < 100). So, the model was run numerous times to generate a vast amount of numerical data that yielded two new correlations based on the velocity and fRe criteria.


Introduction
Microchannel devices, such as heat exchangers [1,2], offer compactness and high heat transfer efficiency, making them promising for the aerospace, computer, and electronics industries. Microchannels also see applications in microreactors [3,4], micromixing, bioanalytical instruments, and lab-on-a-chip. Microchannel hydraulic diameters typically range from 10 µm to 200 µm.
The demand for microchannel devices has significantly increased with the continuous advancements in military and civilian airborne and space-based electronics. To tackle the associated challenges, researchers are actively developing compact microchannel heat exchangers [1,2]. Recent research, exemplified by Bar-Cohen's encyclopedia [5], focuses on fluid dynamics, thermal analyses, and optimization techniques for designing efficient microchannel flow passages. Additionally, researchers are exploring more efficient microchannel designs to handle high-density heat fluxes, such as two-phase flows (particle-liquid or liquid-vapor) and critical heat flux phenomena, using nanofluids in cooling systems [6][7][8]. Thome [9] provides insights into these cutting-edge research topics in the electronics sector.   No deviation from Stokes flow theory was observed within the experimental accuracy. Fully developed velocity profiles were in good agreement with the theory. Developed hydrodynamic entrance length correlations.

K(∞) f app Re
Wiginton and Dalton [24] Analytical X X Fleming and Sparrow [25] Analytical X X X Han [26] Analytical X X X X McComas [27] Analytical X X X Beavers et al. [28] Experimental X X Goldstein and Kried [29] Experimental X Lundgren et al. [30] Analytical X X Shah and London [31] Analytical X Sparrow et al. [32] Experimental X Miller and Han [33] Analytical X Curr et al. [34] Numerical X

Research Objective
Most of the research conducted on rectangular channels occurred during 1960-1970. The researchers usually applied the boundary-layer assumption, covering the applicable laminar Reynolds number range Re ≥ 100 prevalent at their time. With microchannels, the applicable Reynolds number could dip as low as 0.1 for many practical cases. Since there is a lack of information, understanding microchannel fluid mechanics is promising for the future. Therefore, the present study aims to develop correlations for hydrodynamic entrance length for rectangular microchannels for aspect ratios: 1, 0.75, 0.5, 0.25, 0.2, 0.125, and 0.1 ≤ Re ≤ 1000. This broad Reynolds number range will cover the three types of laminar flows described by White [35]: Highly viscous laminar "creeping" motion 1 ≤ Re ≤ 100 Laminar, strong Reynolds number dependence 100 ≤ Re ≤ 1000 Laminar, boundary-layer theory useful The correlations were developed from a significant amount of numerical simulations using a well-known computational fluid dynamic (CFD) software, ANSYS Fluent [36] Version 18.1. The correlations are not only valid for rectangular microchannels, but also for macrochannels under a laminar flow.

Geometry of the Microchannels
A three-dimensional model is shown in Figure 1a for a typical rectangular microchannel. With the assumption of uniform fluid distribution, the system can be analyzed by examining a quarter of a channel, as shown in Figure 1b, by invoking symmetry. The cross-sectional coordinates are x and y, while the axial coordinate is z. The origin is at the center of the channel (bottom left corner of Figure 1b).

Geometry of the Microchannels
A three-dimensional model is shown in Figure 1a for a typical rectangular microchannel. With the assumption of uniform fluid distribution, the system can be analyzed by examining a quarter of a channel, as shown in Figure 1b

 
). Reynolds numbers will vary from 0.1 to 1000 in a logarithmic pattern, as shown in Table 3, for a total of 37 different Reynolds numbers. Six different aspect ratios are considered. The channel's width (2a) is varied to create the various aspect ratios. The variable parameters for the present modeling are the Reynolds number and aspect ratio (α = a/b). Reynolds numbers will vary from 0.1 to 1000 in a logarithmic pattern, as shown in Table 3, for a total of 37 different Reynolds numbers. Six different aspect ratios are considered. The channel's width (2a) is varied to create the various aspect ratios. Note that the high Reynolds numbers can produce impractically high velocities for the prescribed microchannel dimensions. Still, the authors wanted to cover a broad range of Reynolds numbers, and use traditional theory to validate the results due to limited works in the lower or more applicable ranges. Following the criteria of Mehendale et al. [37] (1 µm < D h < 100 µm) and Kandlikar and Grande [38] (1 µm < D h < 200 µm) for microchannels, the largest hydraulic diameter is 100 µm. The hydraulic diameters for each aspect ratio are shown in Table 3. The hydraulic diameter is calculated using traditional theory, Equation (1) [23]. The channel lengths varied based on the required length to determine a fully developed flow, and is further explained in more detail in Section 3.1.

Conservation Equations
For the present study, the governing equations [23,39] are conservation of mass Equation (2) and conservation of momentum Equation (3), assuming that gravitational force is negligible, fluid is incompressible, and the flow is under steady-state condition.
Additionally, the fluid flow is assumed to be laminar and the viscous heating dissipation is considered negligible, which is subsequently proven to be valid in Section 3.3.4.

Numerical Process
The system governing Equations (2) and (3) was solved by the control volume approach using ANSYS Fluent [36]. For the solver setting, the standard pressure-based and steady-state condition were chosen. Using the laminar model, the solution method utilized the SIMPLE scheme with spatial discretization for gradient, pressure, and momentum being least squares cell-based, second order, and second order upwind, respectively.
The models were initialized using the hybrid method, after which the residuals for continuity and velocities were closely monitored. Convergence was achieved for the simulation when the residuals were less than 10 −6 . The model solved for the three velocity components and the pressure throughout the interior computational domain. The data were then exported to MATLAB 2018a [40] for further postprocessing.

Applied Boundary Conditions
The model has a prescribed uniform axial velocity (w) applied at the inlet (z = 0 µm). The specified Reynolds number (Equation (4)) determines the velocity. The pressure outlet boundary condition was applied to the outlet. The no-slip condition (U = 0 m/s) is applied to the channel wall surfaces (x = a, 0 ≤ y ≤ b and y = b, 0 ≤ x ≤ a). Symmetry boundary conditions were applied at x = 0 µm, 0 ≤ y ≤ b, and y = 0 µm, 0 ≤ x ≤ a. The boundary conditions are shown in Figure 2.

Water Properties
The thermophysical properties of water at 308.15 K (35 °C) were used to simulate the fluid. The properties are summarized in Table 4. Several simulations were performed with Al2O3-based ethylene glycol/water nanofluids of 1% and 4% volumetric concentrations. It turned out that no additional information was derived, as our modeling was a single-phase liquid model. A two-phase flow simulation of particles and liquid will show the difference, but that simulation would be much more complex and was not considered in this study.

MATLAB Postprocessing Analysis
The postprocessing analysis was conducted using MATLAB 2018a [40] and determined the following parameters.
The mean velocity [22] is calculated as the average of the cross-sectional plane, as shown in Equation (5), which equals the prescribed inlet axial velocity.
The mean wall shear [22] is determined by averaging the wall shear stress with respect to the perimeter of the channel (Equation (6)) at a given axial location (z location). In Equation (6), x represents the perimeter distance along the wall.

Water Properties
The thermophysical properties of water at 308.15 K (35 • C) were used to simulate the fluid. The properties are summarized in Table 4. Several simulations were performed with Al 2 O 3 -based ethylene glycol/water nanofluids of 1% and 4% volumetric concentrations. It turned out that no additional information was derived, as our modeling was a single-phase liquid model. A two-phase flow simulation of particles and liquid will show the difference, but that simulation would be much more complex and was not considered in this study.

MATLAB Postprocessing Analysis
The postprocessing analysis was conducted using MATLAB 2018a [40] and determined the following parameters.
The mean velocity [22] is calculated as the average of the cross-sectional plane, as shown in Equation (5), which equals the prescribed inlet axial velocity. The mean wall shear [22] is determined by averaging the wall shear stress with respect to the perimeter of the channel (Equation (6)) at a given axial location (z location). In Equation (6), x represents the perimeter distance along the wall.
The local Poiseuille number [22] or Fanning friction factor times Reynolds number can be determined using Equation (7).
The apparent Poiseuille number [22] is calculated using the pressure drop up to the axial location, as presented in Equation (8).
The incremental pressure drop number [22], Equation (9), is the additional pressure loss due to momentum change and the accumulated increment in wall shear between the developing and developed flow [22].
where fRe is the fully developed value and Z* is the dimensionless axial position, Equation (10).
In fully developed flow, the incremental pressure drop number is referred to as K(∞).
Dimensionless hydrodynamic entrance length is defined as where L h is the hydrodynamic entrance length. The velocity ratio is the maximum velocity, U max , divided by the mean velocity, as presented by Equation (12).
The maximum velocity is the centerline velocity computed at the outlet.

Study on the Independence of Channel Length
A length independence study was performed to verify the effects of channel length on hydrodynamic entrance length for a square channel, considering a Reynolds number of 1000 as an example. The study applied the standard fully developed flow criterion of centerline velocity at 99% of U max . The velocity ratio (U*) and the dimensionless hydrodynamic entrance length were determined for four channel lengths (10,20,30, and 40 mm). Comparisons were made with each preceding length to determine when the parameters are independent of the channel length, as shown in Table 5.
From the study, the velocity ratio remains mostly unaffected by the channel length, with the highest difference of 0.26%. The hydrodynamic entrance length dramatically differs between the 10 mm and 20 mm channels by 6.69%. The significant change in entrance length is due to the maximum centerline velocity being an increasing asymptotic function. The slight change to the maximum centerline velocity (0.26%) would affect the hydrodynamic entrance length. At the same time, the difference in length from 20 mm to 30 mm showed a marginal effect on hydrodynamic entrance length. Thus, the initial study had a constant channel length of 20 mm for all aspect ratios and Reynolds numbers. During postprocessing, it was found that most elements were in the fully developed region for low Reynolds numbers and aspect ratios. Thus, results for developing regions were limited. The hydrodynamic entrance length was approximated using the previous data for each subsequent model; the channel lengths were adjusted to 3 times the hydrodynamic entrance length.

Study on the Mesh Independence
A mesh independence study was conducted with channel lengths determined for each model. The modeling adopted orthogonal cartesian grid lines for the xy-plane of the channel cross-section, with smaller elements near the wall as shown in Figure 3. A bias factor increased element size along the axial flow (z-axis) direction. Utilizing the bias factor, 66% of the elements were contained within the suspected developing flow regime.
A mesh study was conducted for each aspect ratio (0.125, 0.20, 0.25, 0.50, 0.75, and 1.0) with a Reynolds number of 0.1 and 1000. Using the lower and upper bounds of Reynolds numbers covers the mesh validation for the Reynolds numbers in between. Due to space limitations, only the square channel mesh (α = 1) is discussed here. The mesh independence study consisted of four configurations, as shown in Table 6. X, Y, and Z elements represent the number of elements in that direction.
Multiple parameters (∆P, L h + , U*, fRe) were computed and analyzed to determine the optimal mesh. Table 7 compares pressure drop (∆P) and dimensionless entrance length (L h + ). The difference is calculated assuming that the tightest mesh, IV, is the more accurate or the so-called "true value".
Micromachines 2023, 14, x FOR PEER REVIEW 10 hydrodynamic entrance length. At the same time, the difference in length from 20 m 30 mm showed a marginal effect on hydrodynamic entrance length. Thus, the initial study had a constant channel length of 20 mm for all aspect r and Reynolds numbers. During postprocessing, it was found that most elements we the fully developed region for low Reynolds numbers and aspect ratios. Thus, result developing regions were limited. The hydrodynamic entrance length was approxim using the previous data for each subsequent model; the channel lengths were adjuste 3 times the hydrodynamic entrance length.

Study on the Mesh Independence
A mesh independence study was conducted with channel lengths determined each model. The modeling adopted orthogonal cartesian grid lines for the xy-plane o channel cross-section, with smaller elements near the wall as shown in Figure 3. A factor increased element size along the axial flow (z-axis) direction. Utilizing the bias tor, 66% of the elements were contained within the suspected developing flow regim    I  15  15  200  45,000  II  20  20  300  80,000  III  25  25  400  250,000  IV  30  30  500 450,000 Mesh I for a Reynolds number of 0.1 displays the highest difference for pressure drop and dimensionless entrance length of 2.08% and 1.25%, respectively. In contrast, the Reynolds number of 1000 for all four meshes agrees with less than 0.70%. These results demonstrate that other than Mesh I, other meshes are reasonably acceptable.
For evaluating the accuracies of velocity ratio (U* = U max /U m ) and fully developed fRe, a comparison is made with known values presented by McComas [27], 2.0962 and 14.227, respectively, as summarized in Table 8. McComas' results were accepted as the "true value" for the difference calculation. As observed in Table 7, Mesh I showed the highest difference of 1.15% and 1.29% for U* and fRe, respectively, while meshes II, III, and IV show differences of less than 0.52% and 0.55%, respectively. From Tables 7 and 8, Mesh IV produces the most accurate results. However, the computational time required for Mesh IV is nearly 2.5 times greater than Mesh III, without significantly increasing accuracy over Mesh III. Thus, Mesh III was chosen as the mesh configuration for subsequent simulations for an aspect ratio of 1. Table 8. Mesh independence study results for an aspect ratio of 1 compared with the theoretical values of McComas [27].

Mesh
Re U* Difference fRe Difference Additional mesh independence studies were conducted for each aspect ratio following the methodology outlined here to arrive at the optimal mesh for subsequent simulations.

Computational Model Validation
The computational model is validated by comparing the velocity profile, U*, and fRe, with published results from analytical analysis and experimental work. The temperature rise due to viscous heating is calculated to determine if viscous heating is negligible, and validated at the end of this section to be true.

Comparison of Velocity Profiles
The fully developed velocity profile for rectangular channels is presented by McComas [27] through Equation (13). Modifications were made to match this paper's geometric orientation, as shown in Figure 4. where Micromachines 2023, 14, x FOR PEER REVIEW 12 of 30

Comparison of Velocity Profiles
The fully developed velocity profile for rectangular channels is presented by McComas [27] through Equation (13). Modifications were made to match this paper's geometric orientation, as shown in Figure 4.
While Equation (13) is accepted to be accurate, it is computationally complex. Although the series converges quickly, 30 terms are required to achieve five significant digits of accuracy. Therefore, due to the complexity of Equation (13), a simpler model was proposed by Purday [42]. Natarajan and Lakshmana [43] expanded upon Purday's model and presented Equation (14), where m and n are functions of aspect ratios.
The fully developed axial velocity profiles from the present numerical results, Equations (13) and (14), for different aspect ratios, are presented in Figure 5, which are in ex- While Equation (13) is accepted to be accurate, it is computationally complex. Although the series converges quickly, 30 terms are required to achieve five significant digits of accuracy. Therefore, due to the complexity of Equation (13), a simpler model was proposed by Purday [42]. Natarajan and Lakshmana [43] expanded upon Purday's model and presented Equation (14a-c), where m and n are functions of aspect ratios.
The fully developed axial velocity profiles from the present numerical results, Equations (13) and (14a-c), for different aspect ratios, are presented in Figure 5 The maximum difference between the present numerical results and the analytical equations is summarized in Table 9. Equation (13) is in excellent agreement with a difference of less than 2%, while Equation (14a-c) has a significantly higher difference, especially for the y-axis velocity, as large as 12.82%. Upon careful data comparison, disagreement occurs near the wall (0.8 ≤ y/b ≤ 1); otherwise, 4% or less was observed in the core region.  (13)) and Natarajan and Lakshmana (Equation (14)). (a) center plane where x = 0; (b) center plane where y = 0.
The maximum difference between the present numerical results and the analytical equations is summarized in Table 9. Equation (13) is in excellent agreement with a difference of less than 2%, while Equation (14) has a significantly higher difference, especially for the y-axis velocity, as large as 12.82%. Upon careful data comparison, disagreement occurs near the wall (0.8 ≤ y/b ≤ 1); otherwise, 4% or less was observed in the core region. Table 9. The maximum axial velocity difference between the present scheme's numerical results and Equations (13) and (14).  (13) and (14) over the entire cross-sectional plane, we observe that Equation (13) has a maximum deviation of 3.9% for the square duct, and less than 2% for other aspect ratios. While Equation (14) is reasonably accurate at the core region, the extremely high difference (≤39%) occurs at the corners.

Comparison of Velocity Ratio
The ratio of maximum velocity to mean velocity (U*) has been presented in the literature by several researchers [25][26][27]30,32,43]. Equations (13) and (14) can be reduced to determine U* by setting x = y = 0 for the centerline velocity, as shown in Equations (15) and (16).   (13) and (14a-c).  (13) and (14a-c) over the entire cross-sectional plane, we observe that Equation (13) has a maximum deviation of 3.9% for the square duct, and less than 2% for other aspect ratios. While Equation (14a-c) is reasonably accurate at the core region, the extremely high difference (≤39%) occurs at the corners.

Comparison of Velocity Ratio
The ratio of maximum velocity to mean velocity (U*) has been presented in the literature by several researchers [25][26][27]30,32,43]. Equations (13) and (14a-c) can be reduced to determine U* by setting x = y = 0 for the centerline velocity, as shown in Equations (15) and (16).
McComas [27] provided a very comprehensive set of results for ten different aspect ratios (1 ≤ α ≤ 0) using Equation (15), while Sparrow et al. [32] experimentally determined the velocity ratio for aspect ratios of 0.5 and 0.2. Their experimental results agreed with the calculations presented by McComas. A comparison between the present numerical results and Equations (15) and (16) is summarized in Table 10. The current numerical computation results agree with Equations (15) and (16), with differences of less than 0.34% and 1.27%, respectively.

Comparison of Fully Developed Friction Factor with Conventional Theory
Several studies have been devoted to the essential practical parameter, the fully developed friction factor. The most widely cited research on this topic could be by Shah and London [31], who developed a correlation for fRe as a function aspect ratio (0 ≤ α ≤ 1), shown in Equation (17).
A rigorous comparison of fRe has been made in Table 11 between the present numerical results and Equation (17) by Shah and London [31] for six aspect ratios. An excellent agreement is observed with the equation of Shah and London with a maximum difference of 0.32%, validating this numerical model.

Evaluation of the Effects of Viscous Heating
Viscous heating was assumed negligible in the numerical model; therefore, validation is necessary. Conservation of energy dictates that the energy loss associated with pressure drop will be converted into thermal energy, which causes the fluid temperature to rise, assuming an adiabatic condition. The energy conversion equation is Equation (18).
where ∆P is the pressure drop, .
V is the volumetric flow, . m is the mass flow rate, c p is the fluid's specific heat, and ∆T is the temperature rise. Equation (18) can be solved for the temperature rise and reduced to Equation (19).
The maximum temperature rise for each aspect ratio would occur with the highest flow rate (Re = 1000). The pressure loss (∆P) is obtained from numerical computation and knowing the thermophysical properties; the temperature rise due to viscous heating is calculated and presented in Table 12. The smallest aspect ratio (0.125) caused the most significant temperature change of 1.69 K due to the high-pressure loss, while other aspect ratios were under 1 K. Considering a temperature change of 1.69 K, the water's density, specific heat, and viscosity would decrease by 0.06%, 0.01%, and 3.29%, respectively. Viscosity shows some change, but is marginal. Therefore, the assumption of negligible viscous heating is valid.

The Criterion of Fully Developed Flow
Various researchers have presented different criteria in determining fully developed flow condition: ratio of centerline velocity to maximum velocity [24][25][26][27]29], or incremental pressure drop number to fully developed incremental pressure drop number [25,28], with ranges of these ratios between 95% and 99%. The most well-known and accepted method for determining the fully developed flow condition is when the centerline velocity achieves 99% of the maximum centerline velocity. This criterion is commonly used due to the ease of measurement from experimental work using image velocimetry [21,32]. The authors will use Equations (20)- (22) for determining the developed flow condition, which will be referred to as the velocity criterion, incremental pressure drop number criterion, and fRe criterion, respectively.
fRe criterion, Equation (22), is a new criterion proposed by the authors. The criterion states that fully developed flow conditions are achieved when f z Re is 99% of fRe. Using the fRe criterion for determining fully developed flow has several merits:

1.
fRe is already well-established from numerous studies by researchers, while fully developed incremental pressure drop number varies between researchers.

2.
fRe is a vital parameter engineers need to know accurately for developing or fully developed flow to determine pressure drop or pumping power requirements. Centerline velocity is less important.

3.
fRe is a parameter describing the fluid-wall interaction, while centerline velocity and incremental pressure depend on the fluid-wall interaction. Thus, fRe will achieve fully developed flow conditions before velocity or incremental pressure drop number.
One of the challenges with using the fRe criterion is that experimental validation would be difficult for measuring the wall shear along the channel with accurate skin friction gages.
A comparison of three fully developed flow criteria is shown in Figure 6 for a square cross-section channel with a Reynolds number of 1000. The velocity and fRe criteria are in reasonable agreement, while the incremental pressure drop number criteria differ greatly. The dimensionless entrance length for velocity and fRe criteria are 0.0733 and 0.0698 at 99%, respectively. Between the two criteria, a difference of 4.72% is observed.
The velocity and fRe criteria were used in all modeling cases to determine the hydrodynamic entrance length. The incremental pressure drop number did not agree with the other criteria, so it was not investigated further. Comparing the results, the difference between them increased with decreasing Reynolds number and aspect ratio.
One of the challenges with using the fRe criterion is that experimental validation would be difficult for measuring the wall shear along the channel with accurate skin fric tion gages.
A comparison of three fully developed flow criteria is shown in Figure 6 for a square cross-section channel with a Reynolds number of 1000. The velocity and fRe criteria are in reasonable agreement, while the incremental pressure drop number criteria differ greatly The dimensionless entrance length for velocity and fRe criteria are 0.0733 and 0.0698 a 99%, respectively. Between the two criteria, a difference of 4.72% is observed. The velocity and fRe criteria were used in all modeling cases to determine the hydro dynamic entrance length. The incremental pressure drop number did not agree with the other criteria, so it was not investigated further. Comparing the results, the difference be tween them increased with decreasing Reynolds number and aspect ratio.

Reynolds Number Influence on Entrance Length
After observing how the criteria compare for a given model of an aspect ratio of 1, a comparison of the velocity and fRe criteria for the entire Reynolds number range is shown in Figure 7. It is a logarithmic plot with a relative difference plot included. The difference between the two criteria increases as the Reynolds number decreases until a Reynolds number of 7, with a maximum difference of 8.06%, which begins to decrease to 5.92% at a Reynolds number of 0.1.

. Reynolds Number Influence on Entrance Length
After observing how the criteria compare for a given model of an aspect ratio of 1, a comparison of the velocity and fRe criteria for the entire Reynolds number range is shown in Figure 7. It is a logarithmic plot with a relative difference plot included. The difference between the two criteria increases as the Reynolds number decreases until a Reynolds number of 7, with a maximum difference of 8.06%, which begins to decrease to 5.92% at a Reynolds number of 0.1.       Figure 7, but with the smallest aspect ratio of 0.125. In this case, the most significant difference is at Reynolds number 3, producing a difference of 64%.

Aspect Ratio Effect on Entrance Length
The aspect ratio also affects the two criteria, as illustrated in Figure 9. At an aspect ratio of 1, the difference is a mere 4%. Still, as the aspect ratio decreases, the difference between the two criteria increases until an aspect ratio of 0.2, then a slight decrease occurs at 0.125.

Aspect Ratio Effect on Entrance Length
The aspect ratio also affects the two criteria, as illustrated in Figure 9. At an aspect ratio of 1, the difference is a mere 4%. Still, as the aspect ratio decreases, the difference between the two criteria increases until an aspect ratio of 0.2, then a slight decrease occurs at 0.125.  A summary of the maximum difference for the different aspect ratios, and the Reynolds number at which it occurs, is presented in Table 13.  A summary of the maximum difference for the different aspect ratios, and the Reynolds number at which it occurs, is presented in Table 13.
The entrance length calculated by the velocity criterion always comes out longer than the fRe criterion. The main conclusion of this investigation is that low Reynolds number and aspect ratio channels need additional studies to establish the cause of the difference. One possible reason for the differences between the velocity and fRe criteria may be the velocity overshoot, which is discussed below.

Velocity Overshoot Observation
Velocity overshoots where the centerline velocity is a local minimum, while local maxima occur near the walls, creating a concave velocity profile instead of a convex one. Velocity overshoots have been studied in detail for circular and parallel plate channels, but to the authors' knowledge, very little has been done for rectangular channels. Figure 10 shows an example of velocity overshoot created from the present numerical results using an aspect ratio of 1, with a Reynolds number of 1000 ( Figure 10a) and Reynolds of 0.1 (Figure 10b). The "X" mark notes the maximum velocity point for that axial location. numbers, and significant at lower Reynolds numbers. Our computational results also found similar behavior, and decreasing the aspect ratio made velocity overshoot more significant.
In Figure 10a, the velocity overshoot starts as the fluid interacts with the wall at the inlet (Z + = 0). Near this location, the velocity overshoot is 1.73% greater than the centerline velocity. The velocity in the overshoot increases until Z + = 2.136 × 10 −4 , showing a velocity of 10.72% greater than the centerline velocity. However, the maximum overshoot for this cross-section occurs near the corner of the duct where two wall surfaces interact with the fluid, creating a velocity 20.21% greater than the centerline velocity, as shown with a contour plot in Figure 11. From this point (Z + ≥ 2.136 × 10 −4 ), the overshoot diminishes as it converges towards the centerline. This happens rapidly for high Reynolds numbers, whereas in Figure 10a, the overshoot converges on the centerline at Z + = 1.145 × 10 −2 . This is roughly 16.40% of the developing region. The effects of a low Reynolds number (Re = 0.10) on velocity overshoot are displayed in Figure 10b. Here, we notice that the velocity difference is lower than that for Re = 1000 of 5.57% at Z + = 0.285 and 15.29% in the corner. The overshoot converges towards the centerline at Z + = 1.602, about 24.08% of the developing region.  In Figure 10a,b, a uniform velocity profile boundary condition is applied at the inlet. Moving down axially, the fluid on the wall obeys the no-slip condition with zero velocity. Thus, the fluid near it accelerates to preserve continuity, creating the overshoot. Moving axially, the velocity overshoot diminishes while converging toward the centerline. From the studies of previous researchers, velocity overshoots are negligible at high Reynolds numbers, and significant at lower Reynolds numbers. Our computational results also found similar behavior, and decreasing the aspect ratio made velocity overshoot more significant.
In Figure 10a, the velocity overshoot starts as the fluid interacts with the wall at the inlet (Z + = 0). Near this location, the velocity overshoot is 1.73% greater than the centerline velocity. The velocity in the overshoot increases until Z + = 2.136 × 10 −4 , showing a velocity of 10.72% greater than the centerline velocity. However, the maximum overshoot for this cross-section occurs near the corner of the duct where two wall surfaces interact with the fluid, creating a velocity 20.21% greater than the centerline velocity, as shown with a contour plot in Figure 11. From this point (Z + ≥ 2.136 × 10 −4 ), the overshoot diminishes as it converges towards the centerline. This happens rapidly for high Reynolds numbers, whereas in Figure 10a, the overshoot converges on the centerline at Z + = 1.145 × 10 −2 . This is roughly 16.40% of the developing region. The effects of a low Reynolds number (Re = 0.10) on velocity overshoot are displayed in Figure 10b. Here, we notice that the velocity difference is lower than that for Re = 1000 of 5.57% at Z + = 0.285 and 15.29% in the corner. The overshoot converges towards the centerline at Z + = 1.602, about 24.08% of the developing region.
The aspect ratio has a significant effect on the velocity overshoot. This has been examined carefully from our numerical results, and displayed in Figures 12 and 13 as a function of the Reynolds number.  The aspect ratio has a significant effect on the velocity overshoot. This has been examined carefully from our numerical results, and displayed in Figures 12 and 13 as a function of the Reynolds number. Figure 12 shows that as the aspect ratio decreases, the difference between the overshoot and centerline velocity decreases. From Figure 12, three zones can be roughly observed: (i) from Re = 10 −1 to 10, the velocity difference is the smallest and is essentially constant; (ii) from Re = 10 to 10 2 , the velocity difference increases and reaches a peak; and Figure 11. Contour plot with velocity overshoot in the corner at Z + = 2.136 × 10 −4 for an aspect ratio of 1 with Reynolds of 1000. "X" marks the location where the maximum velocity is greater than the centerline velocity by 20.21%. Figure 12 shows that as the aspect ratio decreases, the difference between the overshoot and centerline velocity decreases. From Figure 12, three zones can be roughly observed: (i) from Re = 10 −1 to 10, the velocity difference is the smallest and is essentially constant; (ii) from Re = 10 to 10 2 , the velocity difference increases and reaches a peak; and (iii) from Re = 10 2 to 10 3 , the velocity difference decreases and attains a similar magnitude as the low Reynolds number flows (10 −1 to 10). (iii) from Re = 10 2 to 10 3 , the velocity difference decreases and attains a similar magnitude as the low Reynolds number flows (10 −1 to 10). Figure 12. The maximum relative difference between the overshoot and centerline velocity for various aspect ratios as a function of the Reynolds number. Figure 13 shows that the aspect ratio and Reynolds number are inversely related to the distances (Z + ) for velocity overshoot to converge towards the centerline. Thus, for a low aspect ratio, the velocity overshoot distance is higher, while the overshoot is lower. On close analysis, it was found that for an aspect ratio of 0.125, while the velocity overshoot starts in the corner, it quickly moves to the y-axis of symmetry. In this case, the velocity overshoot takes a significant distance to converge toward the centerline. This delay of centerline velocity achieving maximum velocity is why the velocity criterion predicts a longer entrance length than the fRe criterion. Thus, the authors believe that the fRe criterion is more suitable for determining fully developed flow for engineering applications that can accurately predict pressure drop and convective heat transfer coefficient.

Comparison with Previous Correlations
Several prominent authors have researched entrance length for rectangular channels, but most results differ, and there is no universal agreement, as shown in Table 14. All of Figure 12. The maximum relative difference between the overshoot and centerline velocity for various aspect ratios as a function of the Reynolds number. Figure 13 shows that the aspect ratio and Reynolds number are inversely related to the distances (Z + ) for velocity overshoot to converge towards the centerline. Thus, for a low aspect ratio, the velocity overshoot distance is higher, while the overshoot is lower. On close analysis, it was found that for an aspect ratio of 0.125, while the velocity overshoot starts in the corner, it quickly moves to the y-axis of symmetry. In this case, the velocity overshoot takes a significant distance to converge toward the centerline. This delay of centerline velocity achieving maximum velocity is why the velocity criterion predicts a longer entrance length than the fRe criterion. Thus, the authors believe that the fRe criterion is more suitable for determining fully developed flow for engineering applications that can accurately predict pressure drop and convective heat transfer coefficient. (iii) from Re = 10 2 to 10 3 , the velocity difference decreases and attains a similar magnitude as the low Reynolds number flows (10 −1 to 10). Figure 12. The maximum relative difference between the overshoot and centerline velocity for various aspect ratios as a function of the Reynolds number. Figure 13 shows that the aspect ratio and Reynolds number are inversely related to the distances (Z + ) for velocity overshoot to converge towards the centerline. Thus, for a low aspect ratio, the velocity overshoot distance is higher, while the overshoot is lower. On close analysis, it was found that for an aspect ratio of 0.125, while the velocity overshoot starts in the corner, it quickly moves to the y-axis of symmetry. In this case, the velocity overshoot takes a significant distance to converge toward the centerline. This delay of centerline velocity achieving maximum velocity is why the velocity criterion predicts a longer entrance length than the fRe criterion. Thus, the authors believe that the fRe criterion is more suitable for determining fully developed flow for engineering applications that can accurately predict pressure drop and convective heat transfer coefficient.

Comparison with Previous Correlations
Several prominent authors have researched entrance length for rectangular channels, but most results differ, and there is no universal agreement, as shown in Table 14. All of Figure 13. Distance for velocity overshoot to end for various aspect ratios as a function of Reynolds number.

Comparison with Previous Correlations
Several prominent authors have researched entrance length for rectangular channels, but most results differ, and there is no universal agreement, as shown in Table 14. All of the analytical studies [24][25][26][27] applied the boundary layer theory, valid for Reynolds number 100 or greater. Beavers et al. [28] performed an experimental study on sixteen aspect ratios (1/51 to 1) with Reynolds numbers varying from 400 to 3000. They determined the dimensionless entrance length for a square channel was 0.03, and for channels with an aspect ratio less than 0.5, the dimensionless entrance length was 0.015. Goldstein and Kried [29] measured a square duct's developing flow velocity profile using a laser-doppler flowmeter. Most of this work was performed in the 1960s and 1970s for laminar flow in channels, which would be greater than 100 for most "macro" channels. However, the present-day microchannels operate in the 0.1 ≤ Re ≤ 1000 range.  [22] presented the value.
[e] Used the fully developed flow condition: Recently, Ahmad and Hassan [21] conducted microparticle image velocimetry experiments for three square microchannels (500 µm, 200 µm, and 100 µm) over a Reynolds number range of 0.5 to 200. They developed an entrance length correlation, Equation (23), as a function of the Reynolds number. The correlation was within 15% of their data for square channels with a hydraulic diameter of less than 500 µm.
A comparison of hydrodynamic entrance length is made in Figure 14 between the numerical results of the present computation based on velocity and fRe criteria of flow development, and the correlations of Ahmad and Hassan [21], Han [26], and Wiginton and Dalton [24] for a square channel (α = 1). For better visualization, Figure 14 contains three plots expanding across a Reynolds number: (a) 0.1 ≤ Re ≤ 1, (b) 1 ≤ Re ≤ 100, and (c) 100 ≤ Re ≤ 1000.
Our numerical results based on the fRe criterion are in good agreement with Ahmad and Hassan, with a difference of less than 8% (Figure 14a,c), except for Reynolds numbers of 3-30, which increases to 10-21% observed in Figure 14b. The numerical results agreed with Han's results for Reynolds numbers greater than 30, with a difference of less than 7.5%. Note that Ahmad and Hassan and Han overlap for large Reynolds numbers over 100. This is due to Ahmad and Hassan using Han's results in developing their correlation for large Reynolds numbers. The results of Wiginton and Dalton differed more, with the present numerical prediction showing a minimum difference of 9.78% at a Reynolds of 20.
Similarly, our numerical results of entrance length based on velocity criteria for flow development agree with Ahmad and Hassan's correlation, with a maximum difference of 12%. The numerical results agree with Han's results for Reynolds numbers greater than 60, with a difference of less than 3%. The results of Wiginton and Dalton comparatively disagreed more, with numerical results showing a minimum difference of 3.4% at Reynolds number 20, while the rest were above 10%. Figure 15 compares the dimensionless entrance length with varying aspect ratios between the present numerical results (Re = 1000) and Han [26], Wiginton and Dalton [24], McComas [27], and Fleming and Sparrow [25]. As evident from Figure 15, many previous researchers significantly disagree on the entrance length for rectangular channels. From the present research, it can be emphasized that the fRe criterion matches closely (maximum difference: 10%, mean difference: 4.17%) with Han's entrance length prediction. The fRe criterion also seems to show a smooth trend with varying aspect ratios, while the velocity criterion shows a drastic trend change for aspect ratios less than 0.2. Similarly, our numerical results of entrance length based on velocity criteria for flow development agree with Ahmad and Hassan's correlation, with a maximum difference of 12%. The numerical results agree with Han's results for Reynolds numbers greater than 60, with a difference of less than 3%. The results of Wiginton and Dalton comparatively disagreed more, with numerical results showing a minimum difference of 3.4% at Reynolds number 20, while the rest were above 10%. Figure 15 compares the dimensionless entrance length with varying aspect ratios between the present numerical results (Re = 1000) and Han [26], Wiginton and Dalton [24], McComas [27], and Fleming and Sparrow [25]. As evident from Figure 15, many previous researchers significantly disagree on the entrance length for rectangular channels. From the present research, it can be emphasized that the fRe criterion matches closely (maximum difference: 10%, mean difference: 4.17%) with Han's entrance length prediction. The fRe criterion also seems to show a smooth trend with varying aspect ratios, while the velocity criterion shows a drastic trend change for aspect ratios less than 0.2.

An Improved Correlation Was Derived
From the wealth of numerical data generated by numerous simulations for Reynolds numbers varying from 0.1 to 1000 and aspect ratios ranging from 0.125 to 1.0, two correlations for entrance length were developed: one for the velocity criterion, and another for the fRe criterion. A form presented by Chen [44] for circular and parallel plates, Equation (24), which Ahmad and Hassan [21] also recently adopted for rectangular channels, was

An Improved Correlation Was Derived
From the wealth of numerical data generated by numerous simulations for Reynolds numbers varying from 0.1 to 1000 and aspect ratios ranging from 0.125 to 1.0, two correlations for entrance length were developed: one for the velocity criterion, and another for the fRe criterion. A form presented by Chen [44] for circular and parallel plates, Equation (24), which Ahmad and Hassan [21] also recently adopted for rectangular channels, was selected.
In Equation (24), only coefficient B was determined from the curve fit. The lower and upper bounds of the data determined coefficients A and C. Coefficient A is L h /D h at the Reynolds number of 0.1, and coefficient C is L h + at the Reynolds number of 1000. After substantial data analysis using several MATLAB [40] scripts, the correlation coefficients A, B, and C were determined, listed in Table 15. Comparisons of the present numerical data with the proposed correlation are shown in Figure 16. The correlation has a good fit, with a maximum error of 10.28% for an aspect ratio of 0.2, and a mean error of less than 5%. The error means the difference between numerical data and Equation (24). The coefficients A, B, and C yield excellent R 2 values. This table is based on data from all Reynolds numbers studied. The highest error for each aspect ratio usually occurred in the Reynolds number range of 3 to 200, which can be observed as a transition regime from a low Reynolds number (≤1) to a high Reynolds number (≥100) in Figure 16. In Table 16, the coefficients A, B, and C were determined for the fRe criterion, which exhibited an excellent agreement, with a maximum error of 4.39% for an aspect ratio of 0.5. Figure 17 compares the numerical and correlation results using the coefficients of Table 16. In Table 16, the coefficients A, B, and C were determined for the fRe criterion, which exhibited an excellent agreement, with a maximum error of 4.39% for an aspect ratio of 0.5. Figure 17 compares the numerical and correlation results using the coefficients of Table 16. In Table 16, the coefficients A, B, and C were determined for the fRe criterion, which exhibited an excellent agreement, with a maximum error of 4.39% for an aspect ratio of 0.5. Figure 17 compares the numerical and correlation results using the coefficients of Table 16. Using the coefficients A, B, and C determined in Table 16, correlations were developed as a function of aspect ratio. This allows the individual to determine the entrance length for any aspect ratio without interpolating between the data points. For fluid mechanics, the aspect ratio orientation (long side versus short side) for rectangular channels is irrelevant; a channel with an aspect ratio of 0.5 is the same as 2. Thus, the authors wanted to include the "symmetry" nature of the aspect ratio. The "symmetry" nature of the aspect ratio is easily noticeable once the aspect ratio's natural log is considered, as shown in Figure 18. Figure 18a shows the dimensionless entrance length for the Reynolds number of 1000 across a broad range of aspect ratios, while Figure 18b is a similarly drawn plot with the x-axis as ln(α). In Figure 18b, the axis of symmetry can be observed at an aspect ratio of 1 (ln(1) = 0). Additionally, using the aspect ratio's natural log, the aspect ratio's lower and upper bounds define as parallel plates, 0 and ∞ approach −∞ and ∞, respectively. Thus, mathematically, it is easier to include parallel plates. Chen's [44] coefficients for Equation (24) presented by Shah [22] are A = 0.315, B = 0.0175, and C = 0.011. Using the coefficients A, B, and C determined in Table 16, correlations were developed as a function of aspect ratio. This allows the individual to determine the entrance length for any aspect ratio without interpolating between the data points. For fluid mechanics, the aspect ratio orientation (long side versus short side) for rectangular channels is irrelevant; a channel with an aspect ratio of 0.5 is the same as 2. Thus, the authors wanted to include the "symmetry" nature of the aspect ratio. The "symmetry" nature of the aspect ratio is easily noticeable once the aspect ratio's natural log is considered, as shown in Figure 18. Figure 18a shows the dimensionless entrance length for the Reynolds number of 1000 across a broad range of aspect ratios, while Figure 18b is a similarly drawn plot with the x-axis as ln(α). In Figure 18b, the axis of symmetry can be observed at an aspect ratio of 1 (ln(1) = 0). Additionally, using the aspect ratio's natural log, the aspect ratio's lower and upper bounds define as parallel plates, 0 and ∞ approach −∞ and ∞, respectively. Thus, mathematically, it is easier to include parallel plates. Chen's [44] coefficients for Equation (24) presented by Shah [22] are A = 0.315, B = 0.0175, and C = 0.011. The correlations developed by including ln(α) are presented as Equation (25) to Equation (27). Hyperbolic secant functions can represent such symmetric curves. However, Equations (25) and (26) can be reduced to rational polynomials for ease of calculation. For Equation (27), an easy simplification was not possible due to the natural log of the aspect ratio being squared. The correlations developed by including ln(α) are presented as Equation (25) to Equation (27). Hyperbolic secant functions can represent such symmetric curves. However, Equations (25) and (26) can be reduced to rational polynomials for ease of calculation. For Equation (27), an easy simplification was not possible due to the natural log of the aspect ratio being squared.
Using Equation (24) to Equation (27), one can determine the entrance length for channels with an aspect ratio 0 ≤ α ≤ ∞. These correlations have a maximum error of 5.75% and a mean error of 1.87% from the numerically computed data. The relative errors are calculated for several Reynolds number and aspect ratios, and are plotted in Figure 19.

Conclusions
The compilation of earlier published research shows that most studies were on the macrochannel with boundary layer assumption, which may hold for Re ≥ 100. However, most microchannels will operate under Re ≤ 100, and possibly as low as the Reynolds number of 0.1. Therefore, new correlations for determining the hydrodynamic entrance length were necessary. Numerous numerical runs covered Reynolds numbers from 0.1 to 1000, and aspect ratios varying from 0.125 to 1. From the extensive numerical results, an investigation into the fully developed criterion was explored, examining three criteria: velocity, incremental pressure drop number, and fRe. The velocity and fRe criteria were in reasonable agreement for an aspect ratio greater than 0.5 and a Reynolds number greater than 100. The difference between the two criteria for low aspect ratio and Reynolds number may be due to the velocity overshoot. The incremental pressure drop number criterion disagreed with the other two. Thus, this criterion was not pursued any further. Hydrodynamic entrance length correlations were developed for each aspect ratio's velocity and fRe criteria. The velocity criterion correlation (Equation (24) and Table 15) has a mean relative error of less than 5%, and a maximum relative error of 10.28%. The fRe criterion yielded a much better correlation (Equation (24) and Table 16), with a mean relative error of less than 3% and a maximum relative error of 4.39%. Additionally, the fRe criterion correlation was expanded to be a function of aspect ratio using Equation (24) through Equation (27), which is valid for 0.1 ≤ Re ≤ 1000 and 0 ≤ α ≤ ∞, providing a mean relative error less than 2% and a maximum relative error of 5.75%.
Future research must be focused on developing new correlations for friction factors Figure 19. The relative error in hydraulic entrance length prediction via Equation (24) to Equation (27) from numerically computed data.

Conclusions
The compilation of earlier published research shows that most studies were on the macrochannel with boundary layer assumption, which may hold for Re ≥ 100. However, most microchannels will operate under Re ≤ 100, and possibly as low as the Reynolds number of 0.1. Therefore, new correlations for determining the hydrodynamic entrance length were necessary. Numerous numerical runs covered Reynolds numbers from 0.1 to 1000, and aspect ratios varying from 0.125 to 1. From the extensive numerical results, an investigation into the fully developed criterion was explored, examining three criteria: velocity, incremental pressure drop number, and fRe. The velocity and fRe criteria were in reasonable agreement for an aspect ratio greater than 0.5 and a Reynolds number greater than 100. The difference between the two criteria for low aspect ratio and Reynolds number may be due to the velocity overshoot. The incremental pressure drop number criterion disagreed with the other two. Thus, this criterion was not pursued any further. Hydrodynamic entrance length correlations were developed for each aspect ratio's velocity and fRe criteria. The velocity criterion correlation (Equation (24) and Table 15) has a mean relative error of less than 5%, and a maximum relative error of 10.28%. The fRe criterion yielded a much better correlation (Equation (24) and Table 16), with a mean relative error of less than 3% and a maximum relative error of 4.39%. Additionally, the fRe criterion correlation was expanded to be a function of aspect ratio using Equation (24) through Equation (27), which is valid for 0.1 ≤ Re ≤ 1000 and 0 ≤ α ≤ ∞, providing a mean relative error less than 2% and a maximum relative error of 5.75%.
Future research must be focused on developing new correlations for friction factors in the entrance region of microchannels.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.

Conflicts of Interest:
The authors declare no conflict of interest.