An experimental study and modelling of roughness effects on laminar flow in microchannels

Three different approaches were used in the present study to predict the influence of roughness on laminar flow in microchannels. Experimental investigations were conducted with rough microchannels 100 to 300μm in height (H). The pressure drop was measured in test-sections prepared with well-controlled wall roughness (periodically distributed blocks, relative roughness k* =k/0.5H≈0.15) and in test-sections with randomly distributed particles anchored on the channel walls (k* ≈0.04–0.13). Three-dimensional numerical simulations were conducted with the same geometry as in the test-section with periodical roughness (wavelength L). A one-dimensional model (RLM model) was also developed on the basis of a discrete-element approach and the volume-averaging technique. The numerical simulations, the rough layer model and the experiments agree to show that the Poiseuille number Po increases with the relative roughness and is independent of Re in the laminar regime (Re<2000). The increase in Po observed during the experiments is predicted well both by the three-dimensional simulations and the rough layer model. The RLM model shows that the roughness effect may be interpreted by using an effective roughness height keff. keff/k depends on two dimensionless local parameters: the porosity at the bottom wall; and the roughness height normalized with the distance between the rough elements. The RLM model shows that keff/k is independent of the relative roughness k* at given k/L and may be simply approximated by the law: keff/k = 1 − (c(ϵ)/2π)(L/k) for keff/k>0.2, where c decreases with the porosity ϵ.


Introduction
The present study is devoted to the problem of laminar flow in rough-wall microchannels. According to the classical point of view, the surface roughness does not influence the laminar regime of flows in ducts of conventional size. However, significant departures from classical theory have been observed in many works on laminar flows in microchannels during the last two decades. Reviews can be found in Papautsky et al. (1999), Sobhan & Garimella (2001), Morini (2004) and Sharp & Adrian (2004). Among the possible reasons for these deviations, the surface finish of the channel walls has been suspected of playing an important role in the hydrodynamics of microchannels. Recently, precise experimental works on smooth microchannels partly clarified this issue. Judy Bavière et al. (2005) and Kohl et al. (2005), among others, have shown that for smooth-wall microchannels, the friction factor does not depart from conventional values for laminar liquid flows up to Reynolds number of about 2000 where the transition to turbulence occurs. These studies have demonstrated that the classical continuum model (conventional mass, Navier-Stokes equations) is relevant for liquid flows in smooth microchannels with height as low as about 5 µm.
Some authors have experimentally attempted to separate the influence of roughness on the friction factor from other effects such as those due to entrance, cross-section shape and fluid properties (table 1). The relative roughness is normalized with the hydraulic diameter D h . Figure 1 shows the experimental data of the three experimental studies that correspond to a basic duct cross-section (round microtube: Mala & Li 1999;Li, Du & Guo 2003;rectangular microchannel: Pfund et al. 2000). The Poiseuille number (Po = f Re, where f is the Fanning friction factor and Re is defined with the bulk velocity and D h ) is normalized with the value corresponding to a smooth 1.8 Mala & Li (1999) 50 µm FST Pfund et al. (2000) 257 µm Li et al. (2003) 128.8 µm Li et al. (2003) 136.5 µm Mala & Li (1999) 152 µm SST Mala & Li (1999) 101 µm FST Mala & Li (1999)  channel. The results show that the surface roughness increases the resistance to the flow. Shen et al. (2006) also obtained a very large increase in Po for rectangular microchannels, but this was most probably due to the entrance effect so their results are not plotted in figure 1. In Mala & Li (1999), the inlet and outlet pressure losses were eliminated by the experimental procedure. The pressure drop was measured at constant flow rate for two tubes which were identical except of different lengths.
Combining the two sets of data then eliminated the inlet and outlet extra pressure losses. No results were given for smooth microtubes so that a comparison with the reference flow in such tubes and the same experimental loop was not achieved in this study. Li et al. (2003) also performed experiments with deionized water in glass, silicon and stainless steel microtubes. They contrasted the behaviour of smooth glass and silicon microtubes, where the Poiseuille number was consistent with the conventional value for macrotubes, with that of rough stainless steel microtubes, where a significant increase was found for Po. Pfund et al. (2000) measured the pressure drop in deionized water flow across high-aspect-ratio rectangular microchannels made in a sandwich structure. Since the test section could be easily assembled and disassembled, they were able to vary the microchannel height by changing a spacer between the same channel walls and also to perform profilometry of the wall surface. When one of the walls of the 257 µm high test section was a rough polyamide plate, Po exhibited a significant increase (∼25 %) above Po smooth . The trends shown in figure 1 were also found in microchannels with trapezoidal cross-section etched in a silicon substrate (table 1). Wu & Ping Cheng (2003) conducted experiments with a large series of trapezoidal silicon microchannels of different roughness. The reported results show an unexpected increase in Po even for very small values of roughness. It is possible that the data were influenced by uncertainties in the channel dimensions or contain other effects, which obscure the influence of surface roughness. This brief review shows that the influence of relative roughness and Reynolds number is not clearly elucidated by experiments. The accuracy in measured microchannel dimensions obviously plays a crucial role in the experimental uncertainty. It is worth underlining that Po depends to the power four on the tube diameter and to the power three on the channel height so that uncertainties in these dimensions strongly influence the measured Poiseuille number. Another experimental difficulty concerns the characterization of rough walls. In particular, the question arises of where the reference surface is located. This issue was discussed by Kandlikar et al. (2005). Contrary to the well-controlled experiments performed in conventional tubes, the surface of microchannels in certain conditions of etching is more like a bumpy and ridged irregular surface than a smooth surface covered with sand grains as in the experiments of Nikuradse (1933). This obviously results in uncertainty in the hydraulic diameter of the test section for microchannels etched with this method.
Various models have been proposed for accounting for roughness effects in laminar flows. For the extreme situation of a very compact arrangement of rough elements, the roughness region may be considered as impermeable to the flow so that the flow may be calculated with a reduced flow area obtained by subtracting the roughness height from the total cross-section. This simple model may also be relevant in some cases for two-dimensional rough elements. For example, Kandlikar et al. (2005) conducted experiments in a 10.03 mm wide rectangular channel with variable gap. The rough elements were parallel sawtooth ridges 73 µm in height, placed 500 µm apart normal to the flow direction in aligned and offset configurations. The authors showed that the conventional law for Po is recovered in the laminar regime when the hydraulic diameter and the theoretical friction factors are calculated with the constricted flow area.
Perturbation methods have been used for computing the flow over two-dimensional wavy walls (Hocking 1976;Tuck & Kouzoubov 1995;Stroock et al. 2002). In this approach, a modified slip boundary condition can be used on a reference plane instead of the physical no-slip condition on the actual surface. Stroock et al. (2002) equivalently introduced an extrapolation length that defines the plane where the no-slip condition must be satisfied. Sarkar & Prosperetti (1996) presented an approximate analysis for a sparse distribution of arbitrarily shaped protrusions. Another approach was proposed by Mala & Li (1999) who introduced a roughness viscosity model in order to interpret their own experimental data. The increase of momentum transfer due to roughness was modelled by the addition of an empirical roughness viscosity to the fluid molecular viscosity: µ eff = µ + µ R . It was assumed that the roughness viscosity µ R has maximum value near the wall and vanishes gradually towards the channel centre. In this model, the roughness viscosity depends on the Reynolds number to account for the inertial effects observed in the experiments. Koo & Kleinstreuer (2003) also proposed another approach by modelling the near-wall region as an equivalent porous medium layer. They were able to reproduce the measurements of Mala & Li (1999) and  by adjusting the permeability of this layer. Bavière et al. (2006) used an analytical approach derived from the porous medium layer model.
The above review reveals the difficulties encountered in the investigation of roughness effects in microchannel flows. This is why numerical simulation may be a useful way of analysing these effects. The well-defined geometry of the numerical model eliminates the uncertainty inherent in measurements. Moreover, secondary phenomena such as the entrance effect or viscous heating can be easily eliminated. Several authors carried out numerical computations in microchannels with rough elements periodically distributed on a smooth surface (table 2). All these studies found that Po increases with the roughness height. Contrary to the experiments, Po was found to be independent of Re, except in Croce & D'Agaro (2004)  height. Hu, Werner & Li (2003) and Bavière et al. (2006) considered two-dimensional microchannels with three-dimensional rectangular prism elements placed on the walls. They expressed the roughness effect on the pressure drop across the channel as a relative channel height reduction depending on the roughness element geometry. Bavière et al. (2006) used the same geometrical arrangement as that of Hu et al. (2003) and found good agreement with these authors for the effective channel height. In their study, the Re range was extended compared to that used by Hu et al. (2003). They also analysed the contribution of drag forces and viscous stresses on rough elements to the pressure drop. Croce et al. (2005) carried out similar computations for conical rough elements placed on the wall. They observed the effect of the cone steepness at a given geometrical obstruction parameter. Croce & D'Agaro (2004) and Wang, Yap & Mujumdar (2005) considered two-dimensional roughness elements of different shapes.
In the first study, the elements were placed regularly or randomly on the wall. The results indicate that the flow field is characterized by a strong recirculation region developing behind the roughness elements. Croce & D'Agaro (2004) obtained a very large increase in Po ( ≈ 100 %) for the highest value of the relative roughness (0.053). The above presentation shows that roughness effects on microchannel flows are at present not well characterized. There are few well-documented experimental studies on this subject. Roughness effects may be obscured by large uncertainties in the measurements and other effects present in microchannel flows. Some issues remain open, such as whether the Poiseuille number changes with the Reynolds number or about the most significant factors which influence the resistance to the flow.
The current work is focused on fully developed laminar flow through rough microchannels. It combines numerical and experimental approaches to clarify the roughness effects in these flows. A test section with a three-dimensional periodic pattern of rough elements was built and tested. The same geometry was considered in a numerical computation of the flow. A one-dimensional model was also developed for accounting wall roughness in microchannel flows. To the best of our knowledge, it is the first time that such a well-defined rough surface was used in an experimental approach of microchannel flows and compared to a numerical model. 2. Experiments 2.1. Periodically rough surfaces Micro mechanical technologies were used to build two test sections with the same design (figure 2). The microchannel walls were composed of a rough side with periodically distributed elements and a parallel smooth side. The rough wall was obtained by deep reactive ion etching. A silicon wafer was etched at a depth k using a mask reproducing the design of the roughness arrangement. The resulting blocks on the rough wall were parallelepipeds of height k and side length d ( = 8 µm) placed in a staggered arrangement with a wavelength L of 16 µm in both the streamwise and spanwise directions (figure 3). The microchannel height H was obtained by chemically etching  a Pyrex plate at the depth Hk. Two cylindrical holes were made in the Pyrex plate for hydraulic connections. The wafer and the Pyrex cover were then anodically bounded to ensure the watertightness of the microchannel. The two microchannels investigated were 3.5 mm in width and 50 mm in length. The roughness height k was equal to 8.6 µm and 10.6 µm for the like microchannel height H equal to 106.8 µm and 153.6 µm giving the relative roughness (k * = k/0.5H ) equal to 0.16 and 0.14, respectively. Note that k * is normalized with H /2 and not with the hydraulic diameter as in other studies. The two microchannels were tested with demineralized water in a closed-loop circuit, which included a pump (Movichrom N CN 3/12, 10 bars, 20 l min −1 ), a 1 µm filter, three flowmeters (Kobold PEL L45, L01, Bronkhorst LFM L2, range 6 l min −1 , 0.2 l min −1 and 0.017 l min −1 ), two piezo resistive strain gauge transmitters (MBS 3000, 16 bars) and two type K thermocouples for the determination of the inlet and outlet temperature. Details of the set-up can be found in Gao, Le Person & Favre-Marinet (2002). The pressure drop was measured using two pressure taps placed on the Pyrex side of the semi-rough part of the microchannel. The upstream pressure tap was located 15 mm from the inlet. The experiments were also conducted with a smooth microchannel, which was made by the same method as the semi-rough ones, except that the silicon wafer was not etched over the microchannel surface. The flow in the microchannels is characterized by the bulk velocity u b =ṁ/ρH , whereṁ is the mass flow rate per unit length in the spanwise direction and ρ is the density. The Reynolds number is defined by Re = 2ṁ/µ and the Poiseuille number represents the dimensionless pressure drop where µ is the dynamic fluid viscosity. Results are given in figure 4, which shows the following.
(i) Po is independent of Re in the laminar regime (Re < 2000) both for the smooth and the semi-rough microchannels. (ii) For the smooth microchannel, the experimental result for Po is in excellent agreement with the theoretical law of Poiseuille flow between two parallel plates (Po = 24). This confirms previous findings of our group (Gao et al. 2002;Bavière et al. 2005) and others.
(iii) The results obtained for the two semi-rough microchannels collapse onto a single curve. The pressure drop is significantly increased (about 20 %) for a relative roughness of about 0.15.
(iv) The pressure drop suddenly increases when Re is higher than about 2000, indicating transition to the turbulent regime. The increase in Po is faster than that given by the Blasius law.

Randomly rough surfaces
Further experiments were conducted in rectangular microchannels obtained by classical machining technology. The same test sections were used by Gao et al. (2002) for their experiments on scale effects in hydrodynamic and heat transfer in microchannels. The active channel walls were two smooth plane bronze blocks, which were separated by a foil of thickness e f , with a hollowed out central part of width w equal to 25 mm (figure 5). The two blocks were rounded off in the upstream part so as to form a convergent channel entrance. They were hand-polished (arithmetical roughness Ra < 0.1 µm). The thickness of the foil fixed the channel height H , which could be varied in the range 0.1-1 mm by steps of 0.1 mm. The other dimensions of the channel were the width w and the length L c ( = 82 mm). Two sumps were machined in the working section at the channel inlet/outlet. Two pressure transducers were flush mounted at the upstream/downstream sump walls. T-shaped connectors were used to take pressure measurements with a differential inductive pressure gauge (HBW PD1/0.1 bars). The circuit described in the previous section was used to test the microchannels. After a first series of hydrodynamic measurements performed with the smooth walls, the test section was roughened.
The surface roughening consisted of an electrochemical deposition of a nickel (Ni) layer of thickness 2 µm (±0.5 µm), together with small silicon carbide (SiC) particles (5 to 7 µm in diameter), on the two channel walls. A mask was used during the electrodeposition process to prevent the roughening of the surfaces supporting the foils. The two resulting Ni + SiC films were dried with a 10 bars nitrogen gas stream.
Consequently, weakly anchored particles were removed before the series of hydrodynamic measurements, leaving some holes in the nickel layers. A schematic view of the transverse cross-section of a randomly rough microchannel is shown in figure 5.
The sandwich structured arrangement allowed the channel height to be varied by changing the thickness of the foil. The thickness of the different foils was measured using a high-precision comparator (digital MITUTOYO 0.001 mm Micrometer) at respective values of 100 µm (±1 µm), 200 µm (±1 µm) and 300 µm (±1 µm). For the smooth case, the thickness of the foil gave the channel height H directly. For the rough case, H was deduced from the foil thickness by subtracting the total thickness of the electro-deposited nickel layers.
The rough surfaces were carefully analysed. In a first step, top-view numerical photographs of the rough surfaces were taken through a Leica ×200 microscope lens and analysed. The photos revealed black dots representing the SiC particles on a light background corresponding to the nickel layer. The planar extension of the dots typically ranged from 5 to 10 µm. The ratio of the surface occupied by the black pixels to the total surface of the picture (90 µm × 120 µm) was estimated to be about 37 %. This work was completed with a two-dimensional optical profilometer (Hommel Somicronic probe with a 30 nm vertical resolution, a 2 µm planar resolution and 100 mm axial course). The probe used for the measurements was chosen for its low level of sensitivity regarding the physical nature of the material analysed. The topography measurements were done at the end of the series of hydrodynamic measurements in several different locations of the test section. The Ni + SiC rough films were found to be transversally uniform throughout the whole width of the channels (w = 25 mm). In each case, these measurements revealed well-marked steplike profiles between the smooth bronze substrate and the rough Ni + SiC surfaces. The nickel layer thickness was found to be 2 µm (±0.5 µm) for both sides. This value is in good agreement with that deduced from the characteristics (time and current intensity) of the electrochemical deposition. Another important result is that for all the measurements, the maximum height between the nickel layer and the top of a particle (parameter k on figure 5) was measured to be slightly over 5 µm. Finally, the planar concentration of particles deduced from these measurements was found to be in rough agreement with that obtained by the photo analysis. Summing up, three foils of thickness 300, 200 and 100 µm were used for the measurements giving the microchannel heights of 296 µm (±2 µm), 196 µm (±2 µm), and 96 µm (±2 µm) with the relative roughness k * ranging from 0.04 to 0.12. The experiments reported in this paper were repeated to verify their reproducibility. This is particularly important for the channels with Ni + SiC surfaces for which the flow could have removed particles. Fortunately, the tests have shown that the friction characteristics were reproducible, which indirectly demonstrated the robustness of the particles anchorage.
The measurements (figure 6) confirm the results found with periodic rough microchannels. Again it is found that Po is independent of Re up to a transition value of about 2000. A significant increase of Po with the relative roughness is also found with this surface finish. These results will be interpreted in § 5.4.

Numerical model
Three-dimensional computations were conducted using a geometrical model similar to that of the experiments with periodic rough walls . The surface roughness consisted of blocks distributed on the smooth walls of a plane microchannel of very large span w and height H (H w). The roughness elements were parallelepipeds of square cross-section of side length d as in the experiments. They were periodically distributed either in aligned or staggered arrangements ( figure 7). The present numerical model follows the details of that developed successively by Hu et al. (2003) and Bavière et al. (2006). For this reason, only a rapid description of the model is presented in the current paper. Owing to periodicity, the computation domain extended over one wavelength in the streamwise (x) and spanwise (z) directions and owing to symmetry, over the half-channel height in the direction normal to the walls (y) (figure 8). A periodic boundary condition was written for the velocity field at the inlet and outlet sides of the computation domain. Symmetry conditions were chosen at the lateral sides and the no-slip velocity condition at all the solid boundaries.
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. A SIMPLEC (Semi Implicit Pressure Linked Equations Consistent) algorithm was used for the computations. The accuracy and the grid independence of the solution were carefully verified. Further details can be found in Bavière et al. (2006). The results are presented later with those of the rough-layer model presented in next section. Koo & Kleinstreuer (2003, 2004 proposed that the roughness region could be modelled by an equivalent porous medium layer adjacent to a clear fluid layer, namely the central part of the channel. They modelled the additional viscous forces due to rough elements in terms of the medium permeability and they used a nonlinear term for accounting inertia forces. Bavière et al. (2006) used a discrete-element approach initially proposed by Taylor, Coleman & Hodge (1985) to compute the rough-wall skin friction in turbulent flows. They developed an analytical one-dimensional model to compute the pressure drop in a microchannel. This part of the current study is an extension of their work. The present rough-layer model (RLM) combines ideas from the porous medium layer model of Koo & Kleinstreuer (2003) and from the approach of Taylor et al. (1985). The roughness layer is modelled by a periodical distribution of discrete elements placed on a smooth bottom wall, as in the numerical model ( § 3). The total resistance to the flow F T due to the pressure gradient acting on the inlet/outlet microchannel cross-sections is balanced by the force F vw due to friction at the smooth bottom wall, the force F vt due to friction at the top of the rough elements and the drag force F d on the rough elements. Bavière et al. (2006) analysed the contribution of the three components to F T and showed that the contribution of F vw decreases and that of F d increases rapidly when the roughness height is increased. The RLM model directly accounts forF vw since it assumes the no-slip velocity condition at the smooth bottom wall. The drag force is determined by the RLM model by modelling a drag coefficient, while modelling the force F vt constitutes the problem of the boundary condition at the roughness layer/clear region interface. In the current study, we have attempted to improve some aspects of Bavière et al. (2006), especially the interfacial condition between the rough layer and the clear fluid layer and the modelling of the drag coefficient of the roughness elements.

Momentum equation
The method of volume averaging was applied to derive the macroscopic momentum equation in the rough layer (Taylor et al. 1985). The approach considers a control volume (CV) of infinitesimal thickness in the direction normal to the wall (figure 8).
Since the structures are periodic, the extension of the CV is limited to one wavelength λ in the x and z directions (λ = L or 2L for the aligned or staggered arrangements, respectively). A plane parallel to the wall determines the cross-section s(y) on each rough element (s(y) = d(y) 2 for parallelepipedic or pyramidal elements, s(y) = πd(y) 2 /4 for conical elements, where d(y) is the local cross-section side length or diameter). The local porosity is defined as the ratio of area open for flow to the total area ε(y) = 1 − s(y) L 2 . (4.1) The flow is supposed to be fully developed and is modelled as one-dimensional. In the volume-averaging technique, two types of average velocity are commonly used as in porous media (Whitaker 1986), namely the Darcy velocity u D (y) and the effective velocity u(y) which are related by However, the Darcy velocity u D is more relevant than the effective velocity because it is representative of the mass flow rate. Then, the local Reynolds number is defined by An integral formulation of the momentum equation is derived for the CV. The equilibrium of the CV results from the competition between the pressure forces acting on its upstream/downstream sides, the viscous shear stresses acting on its upper/lower sides and the drag force due to the rough element (for details, see Bavière et al. 2006). Periodicity allows us to decompose the pressure into a linear part and a periodic component as in Croce & D'Agaro (2004) p(x, y, z) = dp dx x +p(x, y, z). (4.4) Applying this pressure decomposition, we can separate the pressure forces into a contribution from the overall pressure gradient and from the periodic partp. Considering this latter one, the drag force due to the portion of rough element included in the CV is modelled by using a drag coefficient C d δF d = 1 2 µRe d C d u D δy.
(4.5) C d represents the sum of a pressure coefficient associated to the periodical pressure forces acting on the front and rear sides of the rough element: C dp = p/0.5ρu 2 D where p is the averaged pressure difference between these two sides and a friction factor associated to the friction forces: C df = 2τ /0.5ρu 2 D where τ is the averaged wall shear stress on the lateral sides of the rough element. The drag force due to the pressure gradient component is regrouped with the pressure term in the momentum equation This equation has to merge with the momentum equation for the clear region For the case of parallelepipedic roughness elements, ε is independent of y so that the second viscous term of (4.6) vanishes. This case will be discussed in more detail in § 4.3.

4.2.
Drag coefficient modelling In their model, Koo & Kleinstreuer (2004) introduced an ad hoc value of the permeability in order to fit the experimental results of . In the discreteelement model, Taylor et al. (1985) introduced empirical correlations for the drag force coefficient. However, these expressions relate the drag coefficient to the local Reynolds number only and ignore other geometrical parameters such as the distance between rough elements. This is a weakness of their method, as pointed out by Webb (1994). After a first attempt to relate the drag coefficient to the roughness geometry , the current study was undertaken in order to improve the estimation of this coefficient.
The three-dimensional numerical simulations described in the previous section were used to compute the distribution of the drag coefficient C d (y) along a parallelepipedic roughness element as a function of the local Reynolds number for several values of the geometrical parameters. Figure 9 shows that the drag coefficient C d is inversely proportional to Re d . The proportionality factor C r (C r = C d Re d , resistance coefficient hereinafter) varies neither with the global Reynolds number Re nor with the relative roughness height k * and mainly depends on the porosity ε. It is worth emphasizing that this advantage results from applying the pressure decomposition as given by (4.4). The overall pressure gradient is strongly dependent on the global structure of the microchannel represented by the relative roughness, whereas the drag coefficient C d is directly dependent on the local velocity field and the local geometrical structure of roughness.
For each series of computations, a slight increase of the resistance coefficient C r is observed for decreasing Re d , namely towards the bottom wall. This increase may be attributed to three-dimensional effects due to the rough element/bottom wall interaction. Near the bottom wall, the velocity field between two rough elements is strongly influenced by the bottom wall shear stress resulting in a higher friction factor on the rough element sides. The thickness of the boundary layer on the bottom wall increases with the porosity, thus for the case of large ε(=0.94) we can expect intensification of this three-dimensional interaction. Figure 9 confirms that the increase in C r at low Re d is more pronounced in this case. This increase is not, however, of great importance for the global results of the model since it occurs near the bottom wall where the contribution of drag forces is very small, as for the velocity. We are then justified assuming that the resistance coefficient is constant along the roughness element height. This suggests that we should model the drag coefficient by means of two-dimensional numerical simulations of the flow across a bank of rods instead of by computationally expensive three-dimensional simulations. Such computations were then carried out for the aligned and staggered arrangements of rods with square cross-section (Gamrat, Favre-Marinet & Le Person 2007). The coefficient C r was determined as a function of ε and Re d . The following correlations were obtained from these results by using a least-squares method.
The pressure coefficient C dp and the friction factor C df were distinguished. It was found that C df contributes by 50 % (±5 %) to the total drag coefficient. Further computations were carried out for banks of round tubes in the creeping-flow regime and the results were found to be in excellent agreement with Martin, Saltiel & Shyy (1998). 4.3. Boundary conditions After volume averaging, the problem of fluid flow in rough-wall microchannels reduces to a second-order ordinary differential equation. Equation (4.6) can be solved numerically under appropriate boundary conditions. It was assumed that the velocity is equal to zero at the bottom wall (y = 0) whereas symmetry boundary conditions were used at the channel symmetry plane (y = 0.5H ). The formulation of the boundary conditions at a homogeneous fluid/porous layer interface has been analysed by several authors. The survey study of Alazmi & Vafai (2001) can be used as a reference on this subject. For the discrete-element approach, the continuity of Darcy velocity u D is satisfied in the whole computational domain. For conically shaped elements, the continuity of the effective velocity u(y) and the shear stress can be prescribed without restriction. A particular problem arises when cylindrical or prismatic elements are considered. Bavière et al. (2006) assumed the continuity of the effective velocity gradient at the clear fluid/porous layer interface. Because of the abrupt change of ε at the interface, this assumption leads to discontinuity of the Darcy velocity gradient: However, a discontinuity of the effective velocity gradient is expected at the interface owing to the development of velocity boundary layers on the rough-element top surfaces. In the porous medium approach, Ochoa-Tapia & Whitaker (1995) also introduced a jump condition in the stress at the interface after a careful examination of the averaged continuity and momentum equations in a control volume located at the interface. In order to account for this discontinuity, the present model considers a control volume CV N adjacent to the top of the rough elements. CV N is defined with the same extension as CV (figure 8) in the x and z directions and by k 6 y 6 k + δy. The force induced by the rough element on CV N is restricted to the viscous force at its top horizontal surface. For the case of parallelepipedic roughness elements, the model assumes that the boundary layers developing on the top surface and the lateral sides in the top part of the rough element are similar. As a consequence, the shear stress at the rough element top surface can be deduced from the average friction coefficient C df on the lateral sides of the rough element. As explained before, this friction factor is deduced from the two-dimensional simulations, like the average drag coefficient C d .
Finally, the dimensionless system of equations is: where the lengths and velocities are normalized by H /2 and the bulk velocity u b , respectively, and are denoted by the superscript ()*. C r and C df are deduced from the correlations (4.8). The presence of δy * in the second term of (4.10b) is due to the resulting pressure force and viscous force in the momentum equation being proportional to δy* whereas the viscous force on the top surface of the rough element does not depend explicitly on δy*. Contrary to the case of (4.10a), the slide height of the control volume is not eliminated in the momentum equation applied to CV N leading to (4.10b). The solution of system (4.10) depends on three dimensionless geometrical parameters: the relative roughness height k * , the porosity ε of the roughness region and the ratio L * 2 , which indicates the roughness structure fractionation (1/L * 2 is the number of rough elements in a square of side H/2). The term of the interaction between the flow and the rough elements is inversely proportional to L * 2 so that a high degree of fractionation of the rough elements for a given porosity is associated to low values of L * 2 , in other words to a high resistance of the rough layer to the flow. C r and ε are variable with y for non-cylindrical rough elements.
For the conically shaped rough elements, the system of equations is (4.11b) The following boundary conditions are applied for both systems of equations For the case of parallelepipedic rough elements, the following boundary condition is added to the system: (4.13b) The systems of equations (4.10) and (4.11) were discretized and solved by means of a first-order-finite difference method using the software Matlab. Because of the strong variations of velocity in the normal direction, the rough layer was typically discretized in about 100 slices and the extra layer (k * 6 y * 6 k * + δy * ) in 5 slices.

Numerical results and comparison with experiments
5.1. Parallelepipedic roughness elements The RLM model was applied to a rough microchannel with periodically distributed parallelepipedic elements. The geometry was the same as in the numerical model ( figure 8). The influence of the thickness δy of the CV N on the numerical results was found to be negligible. In fact, the Poiseuille number varied by ±1 % when the normalized thickness δy* was changed from 0.002 to 0.02. Figure 10 shows velocity profiles obtained with the RLM model at constant mass flow rate for three different values of the relative rough element heights k * while the porosity and L * 2 were kept constant (ε = 0.75, 1/L * 2 = 1.56). The Poiseuille profile is shown for comparison. As expected, the velocity profiles depend strongly on k * . For this choice of parameters and for moderate values of the roughness height (k * 6 0.4), the velocity profiles exhibit a nearly linear part in the rough layer. A very high value of k * (= 0.8) has been considered. It corresponds to the case of a heat exchanger with pin fins and not to a rough wall. In this case, the flow in the rough layer is nearly twodimensional, except very near the wall, where the velocity profile has to match the no-slip velocity condition for y = 0. The velocity profile for k * = 0.2 is compared with the three-dimensional numerical solution, showing excellent agreement between the two approaches. Figure 11 presents the Poiseuille number as a function of the relative roughness height for three different values of the porosity. The ratio 1/L * 2 was kept constant to 1.56. On the same figure, we have drawn Po max , which corresponds to the Poiseuille The results may also be presented by using an effective roughness height k eff , as in the analysis of Stroock et al. (2002) or equivalently, a penetration depth (kk eff ) of the driving shear into the rough layer. For a given flow rate, k eff defines the flat surface where the no-slip condition must be satisfied by the Stokes flow in order to give the same Po as the actual rough-wall flow. The analysis of the flow over a two-dimensional sinusoidally modulated surface by Stroock et al. suggests plotting k eff /k as a function of the local parameter k/L. Figure 12 shows that the results collapse onto a single curve with this normalization when the relative roughness k * is varied, ε being kept constant. Moreover, the results are almost identical for the aligned and staggered arrangements (highest difference of about 5 % for small k/L). For ε = 0.75, the side length of a rough element is equal to the distance between two successive elements, so that it is tempting to compare the actual results to those corresponding to a sinusoidally modulated surface. The asymptotic trends indicated by Stroock et al. for small-amplitude modulation of the wall compared to the channel height are shown in figure 12. With the present notation, they correspond to k/L< 1/π, k eff k = 1 2 1 + π k L , (4.14) The results obtained with parallelepipedic rough elements are below the theoretical curve drawn for large wavelengths (equation (4.14)). This is obviously due to the three-dimensional structure of the rough layer which reduces the interactions with the flow compared to the two-dimensional situation of the theory when the rough elements are sparsely distributed. The curve corresponding to k/L 1/π has been plotted with the value (c = 0.56) originally given by Hocking (1976) quoted by Stroock et al. For a dense pattern of parallelepipedic elements, their shape obviously limits the penetration depth of the driving shear into the rough layer compared to a sinusoidally modulated surface and this effect is opposed to the previous one. This may explain why the present results are above the curve drawn with c = 0.56. It is striking that the present results are in very good agreement with (4.15) used with c = 0.4 not only for moderate values of k/L, but also for values of k/L as low as 0.2. For this geometry of roughness, a direct comparison with our experiments ( § 2.1) is possible. We verified that k eff /k as given by the RLM model is not sensitive to the boundary condition on the opposite channel wall (smooth or rough). A good agreement is found between the experimental results and the predictions of the RLM model ( figure 12) which, however, seems to underestimate the pressure drop slightly. Unfortunately, the model's results were not known when the test section was built. Figure 12 shows that the value of k eff /k is only 5-10 % less than 1 for the experimental values of ε (=0.75) and other parameters, so that the rough layer is then not far from being impermeable to the flow. It follows that the numerical results are probably not very sensitive to the model assumptions for these experimental flow conditions. Other experiments with higher values of ε and lower values of k/L are desirable to confirm the validation of the model and are planned for the near future.

Conical roughness elements
The RLM model was applied to a rough microchannel with periodically distributed conical elements (base diameter d 0 , height k). The results are compared with the three-dimensional numerical simulations of Croce et al. (2005). As mentioned before, we assumed that the drag coefficient could be estimated on the basis of the twodimensional simulations, thus neglecting three-dimensional aspects of the flow. The following comparison aims at testing this assumption and finding the range of RLM model application. For conical elements, the porosity varies from its minimal value ε 0 at the bottom wall to unity at the rough/clear interface. The relative roughness height k * was kept constant at 0.106 in the computations. The conclusions of the previous section, however, suggest that the results are not sensitive to this parameter. The same presentation as for parallelepipedic rough elements was adopted in figure 13. The abscissa is related to the cone slope (γ = d 0 /2k) by k/L = √ 1 − ε 0 /2γ . The RLM predictions are compared with the results of Croce et al. (2005). As expected, the penetration depth of the driving shear into the rough layer is much larger than for parallelepipedic rough elements (smallest values of k eff /k). The two models agree, showing that the roughness effect is stronger for more compact roughness, i.e. for small values of ε 0 . This effect is, however, less pronounced than for parallelepipedic rough elements. The RLM model and the numerical simulations of Croce et al. (2005) are in good general agreement, especially for the steeper cones (largest values of k/L at given ε 0 ). This slender shape is obviously favourable to the estimation of the drag coefficient by the two-dimensional modelling used in the RLM model. For a milder cone slope (smallest values of k/L), the RLM model underestimates the roughness impact.

Pin fins
The RLM model used in the current study can be easily adapted to predict the pressure drop in heat exchangers where extended surfaces such as pin fins are used. In the model, the periodically arranged roughness may be considered as the distribution of micro pin fins in a heat exchanger. Kosar, Mishra & Peles (2005) measured the pressure drop in a micro heat exchanger with pin fins and zero tip clearance (k = 0.5H ). In other words, the pin fins occupied the whole channel height in their test-section. The authors built microchannels 100 µm in depth etched in a silicon wafer. The pin fins were circular cylinders 50 or 100 µm in diameter d and 100 µm in length, placed transversely to the flow. Since the diameter d is close to the channel height, threedimensional interactions between the pins and the channel walls (endwall effects) are expected in these experiments. Three different configurations of circular fins giving the same porosity (ε = 0.65) are presented. The fins are distributed in the aligned (4ICL -according to the nomenclature of Kosar et al.: H/d = 1, 1/L * 2 = 0.11) or staggered (2SCL: H/d = 1, 1/L * 2 = 0.11, 3SCS: H/d = 2, 1/L * 2 = 0.44) arrangement. The results of Kosar et al. have been normalized as in the present work and plotted in figure 14 for comparison with the RLM model. The model predictions generally compare well with the behaviour of the experimental results for the highest values of Re. The decrease of the experimental Po number when Re is increased remains unexplained. This effect could be due to experimental uncertainties. Both experimental results and the RLM model confirm that the resistance to the flow is increased with the number of fins (higher values of 1/L * 2 ) for given ε. For the small fin diameter case (3SCS), 1/L * 2 is four times higher than for the large fin diameter case (2SCL). However, the Poiseuille number calculated by the RLM model is only about 3.35 times higher than for the 2SCL case. This is probably due to enhanced three-dimensional endwall effects for this latter case. The resistance to the flow is slightly higher (about 14 %) for the staggered arrangement (2SCL) than for the aligned one (4ICL).

Randomly rough surfaces
Because of the randomness of the surface roughness presented in the second series of experiments ( § 2.2), it was not possible to model precisely the drag forces in the rough region. The RLM model then constitutes an approximation of the physical situation since the real surface roughness has to be replaced by an ordered array of elements in this model. Topographical measurements supply the information about the geometrical parameters of roughness. The most important parameters are the thickness of the Ni layer, the roughness height above this layer and the distance between rough elements. The thickness of the Ni layer determines the location of the smooth bottom wall (as in the experimental procedure). The height of roughness elements was assumed to be equal to 6 µm since the size of SiC particles was 5-7 µm. The porosity at the smooth bottom walls (ε 0 = 0.63) was set on the basis of the microscope observation as described in § 2.2. It implies that for a roughness element size of 6 µm, the average spacing between two rows or columns is equal to 9.9 µm. The RLM model was run with the above values of the parameters for the case of parallelepipedic rough elements in a staggered arrangement. The numerical results are plotted in figure 15. The three experimental points correspond to the plateau of figure 6 and collapse onto a single point in this representation. The error bars have been determined with the extreme values of k: 5-7 µm. The uncertainty is very high, especially on the ordinate, which combines uncertainties on k eff and k. The results of Stroock et al.'s analysis are drawn as in figure 12. The constant c (equation (4.15)) has been adjusted to the RLM results. Figure 15 shows that the RLM model is consistent with the experimental results. It seems therefore that the prismatic shape of the rough elements is relevant to the actual surface finish of the test section. It may be remarked that the randomly distributed roughness is characterized by a rather low value of porosity. The very high values of 1/L * 2 (23.7 -225 for H = 96 -296 µm) correspond to k/L > 0.5 in the present case. As a result of this set of parameters, the penetration depth recorded experimentally and predicted by the RLM model (when run with parallelepipeds) is very small. In other words, the pressure drop may be computed with a good approximation by using the reduced flow passage H -2k in these conditions. It may also be remarked that the real roughness is not uniformly distributed as was observed on the photographs of the surface, which is another source of uncertainty. The arrangement of SiC particles could have formed dead-end paths, resulting in a substantial decrease of the permeability of the rough layer and subsequent increase of the pressure drop.

Conclusions
Three different approaches were used in the present study to predict the influence of roughness on laminar flow in microchannels. Experimental investigations were conducted with rough microchannels 100-300 µm in height. The roughness effect was characterized in test sections prepared with well-controlled wall roughness (periodically distributed blocks, relative roughness ≈ 0.15) and in test sections with randomly distributed particles 5-7 µm in size anchored on the channel walls (relative roughness ≈ 0.04-0.13). Three-dimensional numerical simulations were conducted with the same geometry as in the test section with periodical roughness. A onedimensional model (RLM model) was also developed on the basis of a discreteelement approach and the volume-averaging technique. The closure problem of this semi-empirical model consisted in the determination of the drag coefficient of the rough elements as a function of the geometrical parameters of roughness. The threedimensional numerical simulations revealed that this drag coefficient is constant along the roughness element height. As a result, it was obtained by means of a two-dimensional numerical modelling of the flow across banks of rods or tubes. The RLM model was applied with a jump condition for the shear stress at the clear/rough region interface.
The three approaches were applied to microchannels with periodically rough surfaces. The numerical simulations, the RLM model and the experiments agree, showing that the Poiseuille number Po increases with the relative roughness and is independent of Re in the laminar regime (Re < 2000). The increase in Po observed during the experiments is predicted well both by the three-dimensional simulations and the RLM model.
The analysis revealed that the roughness effect may be interpreted by using an effective roughness height k eff for this simplified situation. This parameter defines the flat surface where the no-slip condition must be satisfied by the channel Stokes flow in order to give the same Po as the actual rough-wall flow. When normalized with the actual roughness height, it depends on two dimensionless local parameters: the porosity at the bottom wall; and the roughness height normalized with the distance between the rough elements. The RLM model shows that k eff /k is independent of the relative roughness k * at given k/L and may be approximated simply by the law: k eff /k = 1− (c(ε)/2π)(L/k) for k eff /k > 0.2. The constant c decreases with the porosity. The results of the RLM model are in good agreement with the three-dimensional numerical simulations of Croce et al. (2005) for walls with conical rough elements. When the model is used for a microchannel with pin fins, a fair agreement is found with Kosar et al. (2005).
For the microchannels with randomly distributed roughness investigated in the present work, the experimental Poiseuille number was found to be close to Po max corresponding to the model of a reduced flow area mentioned above. Although the topography of the rough wall was well-characterized, the roughness effect was difficult to model because uncertainties remain such as the shape of the rough microstructures and the inhomogeneities of their distribution. The present analysis suggests that modelling the roughness structure by periodically distributed parallelepipedic elements is relevant to the actual surface finish.
A final remark concerns the scale effect in microchannels with rough walls. The models presented in this paper did not introduce any physical micro-effect, except for that due to roughness. It means that the results of the numerical simulations and of the RLM model are independent of the channel height, which is used as the length scale to normalize the various dimensions of the rough elements. The good agreement with the experimental results indicates that this assumption is valid for microchannels with height larger than about 100 µm. It is worth emphasizing that the roughness effects studied in the present work are expected especially in microchannels where a significant level of relative roughness can be obtained in some devices owing to the small size of the microsystems.