NUMERICAL SIMULATION OF FLOW THROUGH MICROCHANNELS WITH RANDOM ROUGHNESS

The aim of the study is to determine the effect of a randomly generated rough surface on the laminar flow of a fluid in a microchannel. Two-dimensional axially symmetric microchannels with a circular cross-section in the range of Reynolds number Re = 100-1700 were considered. Flow numerical simulations were performed using the Ansys / Fluent software.


Introduction
The rapid development of miniaturization observed in the last decade has provided new scientific problems and technological challenges. Microdevices and microcomponents are becoming increasing widely used in numerous industrial applications from biomedicine to fuel cells. Most of these devices have microfluidic systems with various functionalities, from process to cooling (Kmiotek & Kucaba-Piȩtal, 2018;Kordos & Kucaba-Pietal, 2018;Marzec, 2021;Zaremba et al., 2018). Unfortunately, the direct translation of the flow models and heat transfer problems developed and tested for macrochannels is not justified for microchannels, because they do not take into account the phenomena significant in the microscale and the discrepancies increase with the reduction of the characteristic dimension and the surface treatment method. For example, as the channel size decreases, there is an increase in the effect of roughness on fluid flow.
While the flows in rough macrochannels are very well known, the micro scale has not been fully researched yet. Therefore, in recent decades, a lot of experimental research has been carried out on microchannels. Some of these studies indicate a change in the nature of the flow with a lower critical Reynolds number value than in the macro scale, as well as a greater roughness coefficient occurring in channels with a small hydraulic diameter. It is suggested by the increase of the roughness effect together with the decrease of the channel size (Dai et al., 2014;Zhang et al., 2010).
The article presents how the flow in the micochannel is influenced by the occurrence of a randomly generated surface roughness, i.e., a set of unevenness (peaks and pits) on the real surface with relatively small spacing between vertices (Kandlikar, 2006;Whitehouse, 2004), In the macro scale, the material from which the element was made, the type of processing, and processing parameters undoubtedly have the greatest influence on the surface roughness. The surface can most often be characterized as a combination of two profileswaviness and roughness (some surfaces also show a shape error). Waviness is measured over a much greater distance than roughness. The difference between the profiles is shown in Fig. 1. To describe the surface structure in the macroscale, various surface parameters are defined by means of numerical values. Many of them do not have a significant impact on the macroscale flow, however, they are particularly important in microchannels, because the considered inequalities are then larger. They can be divided into three groups (Fig. 2) describing the following features: amplitude, spacing, and slope, of which amplitude is of key importance in microflows. The most common and commonly used parameter describing the amplitude is the Ra parameter.
Its value describes the arithmetic mean of the profile deviation from the mean. The calculation of this parameter is shown in Fig. 3. The mean profile line (i) is determined by calculating the mean profile height along the sampling length. All negative deviations are converted to positive (ii). The mean line from the modified profile (iii) is then calculated. The Ra parameter alone is not sufficient to fully describe the surface roughness, but it is most commonly used in the case of microns. The Ra parameter is determined by the formula (Whitehouse, 2004): where:  Surface roughness is the result of the simultaneous interaction of many independent factors, both random and determined, and as a result it has a very complex microgeometry. The size of the surface roughness depends on the type of material and, above all, on the type of its processing. Roughness is caused, among others, by the processes of decohesion, plastic deformation in the cutting zone and formation of chip segments, friction of the tool contact surface against the machined surface, and chip friction against the machined surface (Whitehouse, 2004).
When machining surfaces, even on a microscale, surface roughness is inevitable. There are many types of micro-machining methods with machining precision from 10 −2 μm to 5 μm, such as EDM (electrical discharge machining), ECM (electrochemical machining), pickling, micro-milling and so on (Chu et al., 2014). Depending on the production process, the microchannels can therefore even have a surface roughness in the order of the size of the channel. To properly design flow microdevices, it is necessary to determine the physical laws governing the fluid flow and heat transfer in the microchannel, taking into account the method of making the microdevice and the resulting roughness of the microchannel surface. (Kandlikar, 2006).
According to the knowledge of macrosystems, when the relative roughness is less than 5%, its influence on the coefficient of friction is negligible (Kandlikar, 2006). For microscale channels, experimental and numerical results showed that surface roughness has a significant effect on heat transfer and heat transfer (Ansari & Zhou, 2020;Lu et al., 2020), (Mirmanto & Karayiannis, 2012) . For example, the experiment of Kandlikar et al. indicated that for a 0.62 mm pipe with a relative roughness height of 0.355%, the influence of roughness on the coefficient of friction and heat conduction was significant (Kandlikar, 2005).
In publication (Dai et al., 2014), the authors reviewed 33 scientific articles (a total of 5569 data were collected) for flows in micro-and mini-channels with different wall roughness. The aim was to determine the effect of roughness on the friction coefficients. The authors concluded that if the relative roughness height is <1%, it has little effect on the friction coefficient and the critical Reynolds number. A value of 1% has been suggested as the threshold for distinguishing between smooth and rough micro-and mini-channels.
Despite the availability of a large amount of experimental data, the flows in microscopes are still not fully investigated and described. Due to the scale of the phenomenon, numerical simulations seem to be a good tool. The data obtained through computational methods allows the analysis of many aspects difficult to grasp in an experiment, which will enable understanding of the basic physics of the problem.
Most often in the literature, the surface roughness is modeled with simple geometric shapes, e.g., triangles, rectangles, squares, ellipses, trapeziums (Lalegani et al., 2018;Zhang et al., 2010), (Croce, 2016) and the fractal geometry method (Jia & Song, 2019) or Gaussian function (Pelević & van der Meer, 2016) . However, due to the stochastic nature of roughness, a better approach seems to be to randomly generate surface irregularities, which will be described later in the article.
Literature analyzes show that roughness affects the streamline distribution. With a high roughness value, this can lead to the flow being detached near the wall and the formation of recirculation zones. Detachment of flow and recirculation areas are most likely the main causes of increased friction and the formation of a pressure gradient (Zhang et al., 2010).
The aim of the study was to estimate the effect of roughness on the laminar flow of fluid in microchannels, as well as to select an appropriate method of roughness modeling to achieve the best compliance with the experimental data.

Research methodology
Two-dimensional axially symmetrical microchannels with a circular cross-section were considered. On the wall of the microchannel there are elements simulating the roughness of the channel generated on the basis of the normal distribution, with the density function: The area was obtained by generating a random value from 0 to 1 as a random variable.
The average was the radius of the channel (i.e., 25 μm), while the standard deviation was defined as the product of 1.25 and (Zhang et al., 2010). The radius ri is randomly generated by the above-mentioned function. The variable xi changed in the axial direction by a constant value equal to ( +1 = + ). Having these 2 parameters, it was possible to generate a random (one-dimensional) curve. To achieve the three-dimensional effect, further curves with spacing are generated. The microchannel is modeled as axisymmetric, therefore the curves are averaged (i.e., rays having the same x variable are averaged). As a result, we get "i" points with radius ri and location on the xi axis. This is visualized in Fig.   4. Due to the randomness of the process, Ra of each surface is monitored, then the mean Ra of the generated surfaces is calculated. The final surface is accepted if the mean Ra (calculated from all surface components) does not differ by more than 0.3% from the target Ra. The flow was assumed to be two-dimensional, axisymmetric, incompressible, and steady.
The diameter of the channel was a characteristic dimension for the Reynolds number. The test fluid was water at room temperature (fluid density ρ = 998 kg/m 3 , dynamic viscosity μ = 0.001 Pa·s). The microchannel geometry is presented in Fig. 5.
The parameters adopted for the calculations are presented in Table 1, while the parameter values used for the modeling of the random surface are presented in Table 2. The flow in the microchannel results from the principles of fluid conservation (Kandlikar, 2006): 1. Principles of conservation of mass: The system of equations adopted for modeling can be written in the form of equations [1]: • Equation of continuity: • Momentum equations (Navier -Stokes): The commercial ANSYS Fluent software, using the finite volume method, was used to solve the equations.
The boundary conditions adopted for the analysis are: • At the inlet to the duct, a profile was set up at the outlet of the reference pressure • zero tangential velocity on the channel wall (also impenetrability of the walls), • axial symmetry.
The Reynolds number was determined from the formula (Kandlikar, 2006 To determine the effect of roughness on the flow in the microchannel, the formulas used in macrochannels and milichannels were used. The dimensionless coefficient of friction (Fanning, Dracy) was determined (Kandlikar, 2005). The Fanning friction coefficient (fF) is defined as: The Dracy coefficient of friction depends on: the type of flow, wall roughness, pipe geometry (length, diameter) and is most often determined using the Moody diagram. Dracy's coefficient of friction for laminar flow can be determined based on the Hagen-Poiseuille law (Kandlikar, 2005): In macro and microscale in laminar flow, the coefficient of frictional loss depends on the Reynolds number and not on the roughness (Kandlikar, 2005).
According to (Dai et al., 2014) after calculating the coefficient fD, it is possible to calculate the corrected roughness coefficient matching the experimental data obtained in microchannel tests using, for example, the Gaussian distribution model, which assumes that the roughness behaves in accordance with the Gaussian normal distribution. The correlation was developed based on the standard deviation of the roughness profile. A correction factor has been introduced that modifies the friction coefficient fD. The relative roughness of the discussed correlation is: where: ε [m] -the roughness height of the inner surface of the microchannel.
The construction of the microchannel geometry takes into account the entrance length. The entrance length was calculated from the formula (Kandlikar, 2006): The entrance length of 6 mm was adopted for Re = 2100. The value of the exit length was set at 0.5 mm.
Grid Convergence Index (GCI) relative to the average velocity in the section located in the middle of the channel length was determined to assess the correctness of the mesh [−] -współczynnik zagęszczenia siatki, [−] -the approximation order of the calculations (the assumed value is 2).
In order to optimize the mesh, a study of the GCI coefficient was performed. These studies have shown that with an element size of 0.1 µm, the GCI coefficient is less than 0.4%.
The mesh was generated from elements of the prism and tetrahedron type.

Research results
To investigate the effect of roughness on the flow in the microchannel, the following were analyzed: • changing the shape of the roughness element -as a surface composed of points that meet the Gaussian distribution for the flow of fluid in a channel with a circular cross-section, • a change in the Reynolds number from Re = 100 to Re = 1700.
Based on the calculations, the coefficients of linear loss, pressure loss were determined, and the current lines and velocity vector distributions were generated. Linear loss coefficients were compared to those derived from classical theory, Gaussian random distribution, and experimental data (Mohiuddin Mala & Li, 1999).
The effect of the scatter of results resulting from the pseudo-random distribution of points from which the geometry is built was also investigated.
Due to the use of randomly generated geometry, it was required to create a statistical sample.
For this purpose, 30 sample random surfaces of Ra=1,75 μm (Fig. 6). The Darcy coefficient (f_rand) was calculated from the formula (11), on the basis of the mean of these 30 models, for Reynolds numbers 100-1700. Then, the standard deviation was calculated (Table 3).  The results in Table 3 show that the greatest uncertainty concerns the higher values of Reynolds number, for Re = 1700 -99.8% of the results are within +/-5.3% of the mean value, while for Re = 100 only within +/-2.1%.
Moreover, a comparison the Darcy roughness coefficient was made of the results from the random surface simulation (f_rand), the coefficients resulting from the classical theory (fclassic = 64/Re) and the results of the experimentally obtained coefficients (fexp) (Mohiuddin Mala & Li, 1999), as well as those generated by the Gaussian random distribution formula (fgaussian) calculated from the formula (12). In the equation, fD has been replaced by fclassic. All results were compared to the experimental data (Table 4). Based on the results contained in Table 4 and Fig. 5, it can be concluded that the randomly generated surface is very consistent with the experimental results. The maximum relative difference is about 4% for Reynolds number 1700. For smaller Reynolds numbers the error does not exceed 2.2%. Using the data from the classical theory, it is possible to compute the corrected coefficients of a correlation using a random Gaussian distribution. It gives relatively good results not exceeding 8% error for Reynolds numbers from 0 to 1100.
To visualize the flow and present the influence of the Reynolds number on the flow, the distributions of the streamlines and velocity vectors for the Reynolds numbers 100, 1100 and 2100 were generated for an exemplary configuration selected from 30 models (Fig. 6, Fig. 7). To investigate the pressure drop in the microchannel, 26 measurement planes have been defined. The pressures were related to the cross section at the inlet to the channel. The inlet velocity is 2.01 m/s (Re = 100) and 22.11 m/s (Re = 1100). Ra is 1.75 μm. The results are presented for 3 exemplary models (surf_1, surf_2, surf_3) selected from the group of 30 created ( Fig. 8., Fig. 9.). • Using the method of generating random geometry presented in the article, a very good agreement with the experimental data was obtained, especially in the range of lower Reynolds numbers, where the relative error was only 2.1%.
• The main reason of observed pressure drop is formation of recirculation zones in the roughness pits.
• The one of parameters influencing the length of the recirculation zones are the roughness height. As a result of increasing the speed, the vortices increase with the filling of a given cavity.
• Friction losses decreases as the Reynolds number increases.
• Comparing several methods of calculating the Darcy roughness coefficient, it was found that generating random geometry is the method whose results are closest to the data obtained by experiments. The coefficients resulting from the classical theory generate from 11.1 to 35.5% error (depending on the Reynolds number), therefore it is not the recommended method for calculating the linear loss factor in microchannels. On the other hand, the corrected coefficients using the formula obtained from the Gaussian distribution model achieve an amazingly good result from 2.3 to 8% (for Reynolds numbers in the range 100 -1100). The discussed model works worst for Reynolds numbers in the range 1100-1700, reaching 22% of the relative error.
• The model created on the basis of pseudorandom numbers highlights the nonlinearity of the pressure drop along the channel axis (for higher Reynolds numbers -eg Re = 1100). This is expected on a micro scale where there are abrupt changes in parameters due to an increased proportion of roughness relative to the geometry dimension.