Modeling of Laminar Flows in Rough-Wall Microchannels

Numerical modeling and analytical approach were used to compute laminar flows in rough-wall microchannels. Both models considered the same arrangements of rectangular prism rough elements in periodical arrays. The numerical results confirmed that the flow is independent of the Reynolds number in the range 1-200. The analytical model needs only one constant for most geometrical arrangements. It compares well with the numerical results. Moreover, both models are consistent with experimental data. They show that the rough elements drag is mainly responsible for the pressure drop across the channel in the upper part of the relative roughness range.


Introduction
Much work in recent years has gone into the development of microchannel technology. This has been especially useful in electronics, where ever faster and smaller microchips are built, or in biology where lab on chips are currently used. Investigations of liquid flows in rough microducts present a great interest since rough surfaces are commonly found in many practical situations. The size of roughness elements is usually very small compared to macroscale in ducts of conventional size so that roughness effects may be neglected in laminar flows. However, the relative roughness becomes significant when the duct size is reduced to several hundred micrometers or less, and may play a role in the hydrodynamics, even in laminar flows. Surface roughness is due to the process used for making microchannels, and may be high when micromachining is used. In addition, rough walls are potentially interesting for enhancing heat transfer in industrial applications. It is therefore highly desirable to understand the roughness effect on flows in microducts and to quantify its influence on the pressure losses. This has motivated several studies using experimental, analytical and numerical approaches.
Mala and Li ͓1͔ experimentally found an increase in the Poiseuille number for water flows in rough microtubes ranging from 50 m to 254 m in diameter. For Re below 800, they measured an increase in pressure losses ranging from 7% to 15% for relative roughness ͑roughness element height divided by tube diameter͒ between 1.2% and 3.6%. For Re higher than 1000, the Poiseuille number was found to increase significantly with Re. This was interpreted by the authors as the result of an early transition to turbulence in rough microtubes, when compared to conventional size tubes. Guo and coworkers ͓2,3͔ also performed measurements in rough stainless steel microtubes. They concluded that a relative roughness of 3% to 4% increased the pressure drop from 15% to 37% when compared to the smooth case. Judy et al. ͓4͔ conducted experiments in rough stainless steel microtubes from 75 m to 125 m in diameter. There is no precise indication of the roughness height in their paper. They did not conclude that an increase in the pressure losses might be due to roughness elements. In their study, the Poiseuille number was found to be constant in the laminar regime. Measurements were also carried out in rough trapezoidal microchannels ͓5͔, obtained by wet-etching a silicon substrate and covering it with a Pyrex glass plate. The results revealed an increase in the Poiseuille number compared to the conventional theory in the Reynolds range investigated in the present work ͑Ͻ1500͒. For low Re, the increase in pressure losses varied from 8% to 30% for relative roughness ranging from 1.7% to 3%. These investigations supported the idea of an early transition to turbulence due to roughness elements. Pfund et al. ͓6͔ performed experiments in a rough rectangular microchannel 257 m in height. They found a 26% increase in the pressure losses for a mean relative roughness ͑roughness element height divided by channel height͒ of 1.5%. Recently, Baviere et al. ͓7,8͔ investigated high aspect ratio rectangular microchannels ranging from 100 m to 300 m in height. In their experiments, the relative roughness ranged from 1.7% to 5%. For Re between 50 and 1000, the corresponding increase in the Poiseuille number was found to be 5% to 37%.
To summarize, most published experimental results concerning flows in rough microducts show an increase in pressure losses compared to conventional smooth wall ducts. However, some aspects of these results are conflicting. In some experiments, the Poiseuille number increases significantly with Re. On the contrary, other results show Reynolds number independency of the Poiseuille number in a wide range of Re. It is therefore important to establish the relationship between the pressure drop increase and the rough surface characteristics.
An analytical approach of roughness effects was developed by Mala and Li ͓1͔ and Qu et al. ͓5͔. These authors proposed a roughness-viscosity model in a manner similar to the eddyviscosity concept in turbulent flows. The model brought into play a constant which needed to be evaluated from experimental re-sults. More recently, Koo and Kleinstreuer ͓9͔ modeled the nearwall region by using an equivalent porous medium layer. They were able to reproduce the measurements of ͓1,3͔ by adjusting the permeability ͑or equivalently the Darcy factor͒ of this layer. They also applied the model to microjournal bearings ͓10͔.
Numerical computations were performed by Hu et al. ͓11͔͑ denoted HWL hereafter͒, for water flows in two-dimensional microchannels with rectangular prism rough elements. Their calculations were done in channels ranging from 5 mt o5 0mi n height at Reynolds number ranging from 0.002 to 20. The authors expressed the roughness effect as a relative channel height reduction depending on the roughness elements geometry.
The aim of the present paper is to show that a very simple analytical model is able to predict the roughness effect on the flow in rough-wall microchannels. The analytical model is based on the method developed by Taylor et al. ͓12͔ for predicting the roughwall skin friction in turbulent flows. In this paper, we combined analytical and numerical approaches of the flow and compared the results to experimental data ͓8͔ and to published numerical results ͓11͔.

Numerical Approach
The numerical approach pertains to the roughness model developed by Hu et al. ͓11͔. The fluid physical properties, the basic equations, and boundary conditions follow the details of the model proposed by these authors, except that the explored range of Re number was extended to 200.

Geometrical Model of Roughness.
The geometrical model of roughness consists of blocks periodically distributed on the walls of a plane channel of height H ͑Fig. 1͒. The wavelengths are L r and l r , in the longitudinal x and transverse z directions, respectively. The rough elements are parallelepipeds of square cross section ͑side length d͒ and of height k placed on a smooth wall, called thereafter the bottom wall. They are distributed either in symmetrical or asymmetrical arrangements. Series of calculations were performed by varying the dimensions as follows: Roughness height ͑k͓m͔͒: 0.1, 0.2, 0.3, 0.5, 0.7, 1, 2 Roughness side ͑d͓m͔͒: 0.5, 1, 1.5 Roughness spacing ͑L r = l r ͓m͔͒:2 ,3 ,4 Channel height ͑H͓m͔͒: 5, 7.5, 10 2.2 Computation Domain and Boundary Conditions. The computational domain extends over one wavelength in the x and z directions and over the half-channel height ͑H /2͒ in the y direction normal to the wall ͑Fig. 1͒. The longitudinal and transverse dimensions are thus equal to the rough element spacing. The computational domain was treated as part of an extremely long channel so that the flow was considered as fully developed in the longitudinal direction. As a result, the flow properties are periodic in the x direction, except for the pressure which is composed of a linearly decreasing term and a periodic one. This assumption allows for applying periodic boundary conditions on opposite sides of the domain in the longitudinal direction. The periodic boundary condition may be written for the velocity field U͑x,y,z͒ = U͑x + L r ,y,z͒͑ 1͒ where U is the velocity vector.
The lateral sides of the computational domain are chosen as the planes of symmetry of the rough elements parallel to the xy plane. As a result, the velocity gradient of U and V normal to these side planes and the transverse velocity component W are assumed to be zero.
The condition of symmetry at the top plane of the computational domain is written The system of boundary conditions is completed by the no-slip condition on all the solid boundaries.

Equations.
Modeling the flow as laminar, steady and incompressible, the governing equations consist of the following system: Continuity equation

١U =0 ͑4͒
Momentum equation 2.4 Numerical Scheme. Numerical computations of the flow were carried out by using the commercial code FLUENT 6.1.22. The equations were discretized by means of a second order upwind finite volume method. As these equations are nonlinear, a SIM-PLEC ͑semi-implicit pressure linked equations consistent͒ algorithm was used. This algorithm is based on a prediction-correction method, which allows the equations to be linearized and solved iteratively. The pressure under-relaxation factor was set to 0.5. The calculations were performed by means of a double precision solver until the level of residuals decreased below 10 −12 . Computations were performed on orthogonal grids generated by GAMBIT 2.1.2.

Numerical Accuracy.
Grid-convergence tests were performed to verify the mesh accuracy. The tests were conducted with three grids having a number of mesh nodes equal to 32 ϫ 40ϫ 32, 48ϫ 60ϫ 48, 64ϫ 80ϫ 64 in the x , y , z directions, respectively. The pressure gradient was chosen as the control parameter during the computations. The difference in the results as given by the coarse and the intermediate grids was 1.5%, and was 0.7% between the intermediate and the fine grids. The accuracy was further estimated by means of Richardson extrapolation ͓13͔. When two approximations of the solution T 1 and T 2 are computed with two mesh sizes h 1 ͑fine grid͒ and h 2 ͑coarse grid͒, one can estimate a third approximation T 3 whose term of leading order in the Taylor series expansion around the exact value is of higher order than that for T 1 or T 2 . T 3 is given by the Richardson extrapolation

Fig. 1 Computational domain and arrangement of the rough elements
The difference in the pressure gradient as obtained with the intermediate mesh and that obtained from the Richardson extrapolation applied with the intermediate and fine grids was 1.6%. This observation allowed us to assume that the numerical calculations were sufficiently accurate with the grid 48ϫ 60ϫ 48, which was adopted for all the cases.

Analytical Approach
The analytical approach is based on the method developed by Taylor et al. ͓12͔ for predicting the rough-wall skin friction in turbulent flows. It considers the same geometrical model of roughness as the numerical analysis. The flow in a half channel consists of two adjacent layers: A roughness layer with distributed rough elements in the near-wall region and a clear layer in the central region. Like in ͓12͔, the physical effects of roughness are taken into account by using a blockage coefficient and by considering drag forces exerted by the rough elements on the fluid flow. 3.1.1 Blockage Coefficient. The flow velocity is averaged over a surface opened to the flow at the distance y from the wall. Thus, the resulting spatially averaged velocity U͑x , y , z͒ takes into account the blockage effect due to the presence of the rough elements. For periodically distributed blocks, U is dependent on y only and it is sufficient to integrate the velocity over one wavelength in the x and z directions only. U͑y͒ is then defined by where u͑ , y , ͒ stands for the actual flow velocity in the roughness region ͑y ഛ k͒.
and is therefore This is the same expression as in ͓12͔, except that the cross section of the rough elements is square instead of being circular in ͓12͔. In the present model of cylindrical elements, ␤ is independent of y. The maximum blockage effect corresponds to ␤ =0, when the rough elements are adjacent.
We also use a filtration or Darcy velocity u D ͑y͒, in the terminology of flow in porous media. This velocity is defined with the total surface L r l r u D ͑y͒ = 1 L r l r ͵͵ dS͑x͒ u͑,y,͒dd ͑10͒ so that The effective velocity U is higher than the Darcy velocity because it refers to a smaller surface than u D . In the most general case, the average velocity U would be a function of x and z. However, the flow is assumed to be two dimensional and fully developed so that U depends on y only. As a result, the transverse velocity is zero and the inertia terms are also zero in the momentum equation.
In the present approach, the pressure is spatially averaged in the roughness layer like the velocity. However, the result depends on the position x of the integration surface because of the pressure gradient in the channel.
We suppose that p͑x , y͒ is independent of y in the present mean one-dimensional flow. From here, we call P͑x͒ the averaged pressure defined by Eq. ͑12͒.

Momentum Equation.
The equilibrium of a control volume results from the competition between Pressure forces on the upstream and downstream sides of CV Viscous shear stresses on the upper and lower sides of CV Drag forces due to the rough elements The drag force due to the portion of a rough element included in CV is modeled by using a drag coefficient C d Finally, the momentum equation is The expression of the pressure term results from the definition ͑Eq. ͑12͒͒. It may be interpreted by considering that the pressure forces act on a reduced area owing to the blockage effect. In the same way, a reduced area of the CV is exposed to the viscous shear stresses. Equation ͑14͒ is identical to the momentum equation of ͓12͔, except that there are no inertia terms and no Reynolds shear stresses in the present case.

Boundary Conditions.
A difficulty arises for the rough/ clear layer interfacial conditions for the present shape of the roughness elements. In fact, the blockage factor is discontinuous at the interface, going from 1 to a constant value ␤ Ͻ 1 across the interface. This difficulty would not exist for other shapes of the roughness elements, like conical-shaped blocks, as in the Taylor et al. paper ͓12͔. The continuity of the local velocity across the rough/clear layer interface leads to the continuity of the Darcy velocity. As a result, the spatially averaged velocity U͑y͒ is discontinuous like the blockage coefficient at the interface.
As mentioned before, this discontinuity would not exist for roughness elements of gradually decreasing cross section. In fact, it is likely that the blockage effect extends over a small distance within the clear medium away from the interface because the fluid is strongly slowed down in the vicinity of the top of the blocks. The present model does not account for this effect and this is probably a weakness of the model. It could probably be improved by introducing a smoothing of the blockage effect near the interface. This was not attempted in the present work to preserve the simplicity of the model.
In order to express the boundary condition on the shear stress at the interface, we consider a control volume of height ␦y and lengths ␦x, ␦z in the other directions, centered at y = k. The force balance on this CV requires When ␦y tends to zero, the pressure force and the drag force become negligibly small so that Eq. ͑16͒ reduces to or using the Darcy velocity Again, the discontinuity of the shear stress expressed with the Darcy velocity would not exist for an interface with continuously varying blockage factor. The system of boundary conditions is completed by the no-slip condition at the wall bottom and the condition of symmetry on the channel axis in the case of identical rough walls on both sides of the channel.
For a semirough channel formed by a rough wall on one side and a smooth wall on the opposite side, the latter condition is replaced by the no-slip condition on the smooth side ␦F d is also related to the drag coefficient C d by Eq. ͑13͒ so that Eliminating C d by using Eqs. ͑13͒ and ͑22͒, the last term of Eq. ͑14͒ becomes 1 / 2KU / L r l r . It is well known ͓14͔ that inertia effects cannot be totally discarded for creeping flows near extremely long cylinders. However, the numerical simulations of the present work and those of HWL along with some experimental works ͓4,7͔ indicate that roughness effects are Re independent up to Re D h = 200. Consequently we assume that K is a function of y / k and k / d only. The presence of these parameters in the coefficient K reflects the threedimensional nature of the low-Reynolds number flow near the blocks. For creeping flows normal to the major axis of very slender prolate spheroids ͑E ӷ 1͒, the drag force is given by ͓15͔ where E = a / b,2 a = major axis, 2b = minor axis of the spheroid. Following the form of Eq. ͑23͒, we propose for K the following expression: The drag force is, however, also dependent of the geometrical arrangement of the rough elements like in a bank of tubes in a cross-flow exchanger. The coefficient K 1 was therefore deduced from the numerical simulations conducted in typical arrangements of the rough elements. For a given flow field, the elementary drag force ␦F d was obtained by integrating the pressure forces on the front and rear sides of a block slice of height ␦y and the viscous forces on its lateral sides. The numerical results have shown that the coefficient K 1 is nearly independent of y in the upper part of the blocks, as presented later. For k / d ഛ 1, K 1 was found to be 116 for the following geometrical arrangement of the rough elements ͑ d / H = 0.2, 0.1Ͻ k / H Ͻ 0.2, L r / d = l r / d =2͒. For k / d Ͼ 1, the best agreement of the model with the HWL results was obtained for K 2 = 0.5. The dependence of K on the geometrical parameters was not systematically explored in the present work. The same constants were used for all the geometrical arrangements considered in this study. The model could probably be improved with a better estimation of K as a function of the actual arrangement of the rough elements.
3.1.5 Analytical Solution. The system of equations is reduced to dimensionless form by normalizing the lengths with the channel height H and the effective velocity with the bulk velocity of the Poiseuille flow in a smooth channel Modeling K as independent of y, straightforward computations lead to the following equations for the velocity profile in a channel with two rough sides: where A * =1/2K / L r * l r * ␤, C 1 = cosh͑ ͱ A * k * ͒, S 1 = sinh͑ ͱ A * k * ͒ For a semirough channel, the expressions are Integration of the velocity profile gives the flow rate Q as a function of the pressure gradient, the geometrical parameters k * , ␤, and of A * .

Results
The numerical computations were conducted for Re D h varying between 1 and 200. All the results were found to be Re independent. This was also observed by HWL in the range of Re D h = 0.002-20. The symmetrical and asymmetrical arrangements of rough elements gave identical results.

Drag Coefficient of the Roughness Elements.
As can be seen in Fig. 1, the blocks were divided into separate slices of height ␦y for calculating the local force density exerted by the fluid as a function of the local effective velocity U͑y͒. This approach allowed for computing the coefficient K͑y͒ defined by Eq. ͑22͒ from the numerical results. Figure 4 displays the variation of K with the position along a rough element ͑y = 0 corresponds to the bottom wall͒ for a typical value of k * ͑=0.2͒. It shows that K slightly varies in the upper part of the block ͑y / k Ͼ 0.3͒. The strong increase in K observed for the lowest part of the block is obviously due to the influence of the bottom wall and associated three-dimensional effects. The same computation used with a condition of slip on the bottom wall gave a constant value of K all along the block. However, this deviation of K from a constant value is of little consequence on the solution given by the model because the contribution of the drag force to the flow dynamics is of minor importance in the near-wall region as shown in the next section. It follows that the assumption of constant K used in the analytical approach is quite reasonable. Figure 5 compares the velocity profiles given by the analytical and numerical solutions for a rough and a semirough channel. The analytical model slightly overestimates the velocity across the channel. In fact, the deviation is only 3% to 3.5% in the central region of the channel for the two cases studied. It can be concluded that there is a very good agreement between the two solutions, given the crude assumptions of the model and the numerical accuracy. The dimensionless Poiseuille profile is also drawn in Fig. 5 for comparison. The three curves correspond to the same pressure gradient. The presence of roughness significantly reduces the mass flow rate, as expected. The reduction in the maximum velocity is about 28% for the fully rough channel and approximately the half ͑13%͒ for the semirough channel. For this latter case, the velocity profile merges with the smooth channel curve when y * tends to 1 and with the rough channel one when y * tends to 0.

Velocity and Shear Stress Profiles.
Both the analytical and numerical models enable the separation of the role of viscous and drag forces in the flow dynamics. Considering a control volume of height H / 2 based on an elementary cell of sides L r ϫ l r ͑Fig. 1͒, the force F T due to the pressure gradient is balanced by the tangential force F w due to friction on the bottom wall ͑surface ␤L r l r ͒, the tangential force F t due to friction at the top of the rough elements ͑surface ͑1−␤͒L r l r ͒ and the drag force F D on the rough elements. Figure 6 shows the contribution of the three components to the total force F T as a function of the dimensionless roughness height for a fully rough channel. For very small rough elements, the pressure gradient is mainly due to the bottom wall viscous force. This contribution decreases rapidly when the height of the rough elements is increased. In the case of high rough elements, the fluid is strongly slowed down in the rough layer ͑Fig. 5͒ and the velocity gradient decreases at the wall. At the same time, the rough elements merge in a region of higher velocity and the drag force increases accordingly. For k * Ͼ 0.03-0.06, the contribution of the drag force to F T is dominant. The contribution of F t tends to 1 − ␤ when k * tends to zero. It remains approximately constant up to k * Ϸ 0.1. This may be understood as follows: The top of the blocks is located in a region of increasing velocity and is subjected to increasing velocity gradient when k * is increased. In these conditions, F t still contributes to the pressure gradient unlike F w which corresponds to a weakening velocity gradient when k * is increased. For high rough elements ͑k * Ͼ 0.1͒, the contribution of F t decreases and the pressure gradient is mainly due to the drag force. The numerical and analytical models give the same trend for the various phenomena. However, significant differences are present from a quantitative point of view. In particular, the contribution of F t is much weaker in the analytical model. This may be due to the discontinuity of ␤, as discussed before. It is likely that the blockage effect is not completely accounted for by the analytical model so that the velocity gradient is underestimated near the top of the rough elements by this model. In fact, the model ignores the development of boundary layers on the top of the rough elements so that the shear stress is underestimated for y = k. This may also explain why the velocities as given by the analytical model are higher than the numerical ones.

Influence of the Roughness Elements Geometry on the Poiseuille Number.
For given channel height and roughness geometry, the flow is characterized by the Poiseuille number Roughness effects are conveniently interpreted by introducing an apparent channel height H r which corresponds to a smooth channel giving the same flow rate as the rough channel when it is submitted to the same pressure gradient. H r satisfies the following relation for fully developed Poiseuille flow in a plane smooth channel As a result, H r is deduced from the Poiseuille number as given either by the analytical model or the numerical simulations for given channel geometry.
The following figures compare the relative channel height reduction 1 − H r / H as given by the analytical model and the numerical results of the present work. HWL interpolated their numerical results and proposed equations for the channel height reduction as a function of the geometrical parameters. Their results are also displayed on the following graphs.

Influence of the Roughness Element
Height. The apparent channel height reduction is moderate up to k * Ϸ 0.05 ͑Fig. 7͒, then increases rapidly with k * . The graph has also been drawn for very high values of k * ͑0.4͒ which are obviously very far from usual physical situations. There is an excellent agreement between the three approaches. The slight differences of the present numerical results with those of HWL may be explained by the fact that HWL interpolated their results and that inaccuracies may result from this approximation. The expected maximum value for 1 − H r / H is obviously 2k * and is displayed in Fig. 7. The figure reveals that the actual blockage effect is about 25% smaller for the present arrangement of rough elements when k * is high ͑Ͼ0.1͒. The blockage effect is less pronounced for small rough elements.

Influence of the Roughness Element Side Length.
The blockage coefficient directly depends on the roughness element side length d for given wavelengths in x and z directions ͑Eq. ͑9͒͒ and is drawn in Fig. 8. The apparent channel height reduction increases with d * for given k * , as expected ͑Fig. 8͒. The maximum value as given by the analytical model is reached for d *2 = L r * l r * ͑d * = 0.4 in the present case͒ and is close from 2k * . The three models show a continuous decrease of 1 − H r / H when ␤ is increased. Small values of d * correspond to slender blocks placed at relative long distances L r / d * = l r / d * apart. In this case, the HWL results indicate that the roughness effect is still significant although the blockage coefficient becomes very close from 1. This suggests that in this situation the drag force F D exerted by the rough elements still plays an important role in the flow dynamics. The analytical model is consistent with this statement since the drag force ͑or A * in Eq. ͑26͒͒ is independent of d when the drag coefficient K is taken as a constant. However, the model used with K 2 = 0 significantly overestimates the channel height reduction as given by HWL for d * Ͻ 0.1. The agreement is much better when the drag coefficient is used with K 2 = 0.5 ͑Eq. ͑24͒͒. When d * tends to zero, the aspect ratio of the blocks becomes very high and the correction factor K 2 Ln͑k / d͒ in Eq. ͑24͒ plays a more important role in the drag. This suggests that three-dimensional effects must be taken into account in F D in this case.

Influence of the Roughness Element
Spacing. The other way to vary the blockage coefficient is to change the rough elements spacing while keeping their shape constant ͑k * = d * = 0.2 in Fig. 9͒. Unlike the preceding case, the channel height reduction continuously decreases from the maximum value 2k * to zero when the rough element spacing is increased. This result suggests that the role of the drag force becomes less important when the rough elements spacing is increased. This is consistent with the analytical model where the term corresponding to the drag force in Eq. ͑14͒ is inversely proportional to L r l r . The roughness effect becomes vanishingly small when the density of the roughness elements tends to zero while keeping their shape constant. Again, Fig. 9 reveals an excellent agreement of the analytical model with the numerical simulations.

Comparison With Experiments.
In parallel with the analytical and numerical approaches, experiments were conducted in water flows through rough rectangular microchannels ͓7,8͔.I n these experiments, the channel walls consisted of parallel bronze blocks separated by a foil of thickness e f ͑0.1 mmϽ e f Ͻ 1m m ͒ with a hollowed out central part 25 mm in width. The length of the microchannel was 82 mm. The blocks were locally treated by electrochemical deposition of a thin Ni layer ͑thickness e Ni Ϸ 2 m͒, together with small SiC particles ͑5t o7m in height͒. The height of the roughness elements was estimated to be the size of the deposited particles ͑5 m͒ and their spacing was found by visualizations to be around 15 m. The channel height was given by e f −2e Ni . The pressure drop was measured by two pressure transducers flush mounted at the wall of sumps located upstream/ downstream of the channel. The experimental surface finish was obviously far from the regular arrangement of the proposed models. However, the measured pressure drop may tentatively be compared to the predictions of the analytical model. Figure 10 displays the experimental Poiseuille number Po exp deduced from the measured pressure drop and flow rate and the model's results obtained for three values of the rough elements height ͑6.2, 7.2, and 8.2 m͒. The experimental uncertainties are also shown in Fig.  10. The analytical model was used with the following relations: k * = d * = L r * /2=l r * / 2 for three values of k. Numerical computations were performed with k = d =5 m, L r = l r =15 m and the results are also displayed in Fig. 10. The trend exhibited by Po exp is well reproduced by the analytical and the numerical models. Considering the analytical model, the best agreement with the experimental results is obtained for k = 7.2 m, which is slightly higher than the estimated height of the roughness elements. The same tendency may be deduced from the numerical model: A better agreement with the experimental data would have been obtained for higher roughness elements. Apart from the experimental and numerical errors and the uncertainties inherent in the analytical model, this discrepancy may also be due to the experimental surface finish. In fact, it is possible that the layer of nickel was not perfectly smooth but may have included small rough elements besides the largest particles of SiC. As a consequence, the assumption of a smooth bottom wall both in the analytical and numerical models would underestimate the friction on this part of the wall.

Conclusions
Numerical modeling and analytical approach were used to compute laminar flows in rough-wall microchannels. Both models considered the same arrangements of roughness elements in periodical arrays. The numerical results have confirmed that the flow is independent of the Reynolds number in the range 1-200.
The analytical approach was based on a previous work ͓12͔ devoted to the prediction of the rough-wall skin friction in turbulent flows. The analysis has pointed out the difficulties of setting the boundary conditions at the clear/rough layer interface in the case of rectangular prism rough elements. For given geometrical characteristics of the roughness elements, the analytical model only needs two constants K 1 and K 2 to compute the flow. They correspond to the drag coefficient of the rough elements. In most cases, a good approximation of the flow characteristics is obtained by using only one constant, namely, K 1 . Despite the simplicity of the model, a very good agreement is found with the numerical results of the present work and those of Hu et al. ͓11͔. The results show that the rough elements drag is mainly responsible for the pressure drop across the channel in the upper part of the relative roughness range ͑k * Ͼ 0.03-0.06͒. For very high rough elements ͑k * Ͼ 0.1͒, the wall friction has a negligible contribution to the pressure drop. For a concentrated arrangement of the rough elements ͑d * = L r * /2=l r * / 2 = 0.2, Fig. 7͒, a very simple approximation of the roughness effect consists in considering a relative channel height reduction equal to 2k * . The accuracy is about 25% in 1 − H r / H for k * Ͼ 0.1. This blockage effect reduces when the rough elements are small or more dispersed ͑increasing blockage coefficient ␤͒. However, it remains significant for slender blocks, even if the side length becomes very small.
This study was restricted to patterns of periodic rough elements. The good agreement with experimental results obtained with randomly distributed roughness on microchannel walls is therefore striking. It seems that the models catch the main features of the physical phenomena involved in these microflows. The two main ingredients of the analytical model are the blockage coefficient ␤ and the relative roughness height k * . The knowledge of these parameters seems therefore to be of primary importance to predict the roughness effects in rough microchannels. This study will be soon extended to roughness effects on heat transfer in microchannels.