# Modelling of Laminar Flows in Rough-Wall Microchannels

Roland Bavière, Gabriel Gamrat, Michel Favre-Marinet, Stéphane Le Person

► To cite this version:

Roland Bavière, Gabriel Gamrat, Michel Favre-Marinet, Stéphane Le Person. Modelling of Laminar Flows in Rough-Wall Microchannels. Journal of Fluids Engineering, 2006, 128, pp.734-741. 10.1115/1.2201635. hal-00204663

HAL Id: hal-00204663

https://hal.science/hal-00204663v1

Submitted on 13 Feb 2020

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# 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.

**R. Bavière**
Laboratory of Geophysical and
Industrial Flows, Grenoble University,
UJF-INPG-CNRS BP 53X,
38041 Grenoble Cedex 9, France and
Centre de Recherches sur les Très
Basses Températures,
CNRS B.P. 166,
38042 Grenoble Cedex 09, France

**G. Gamrat**
**M. Favre-Marinet<fnref n="1" />**
**S. Le Person**
Laboratory of Geophysical and
Industrial Flows, Grenoble University,
UJF-INPG-CNRS, BP 53X,
38041 Grenoble Cedex 9, France

## 1 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, Bavière 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 eddy-viscosity concept in turbulent flows. The model brought into play a constant which needed to be evaluated from experimental re-

<footnote n="1">¹Corresponding author.</footnote>
<figure><img src="image_2.png" /><figcaption>Fig. 1 Computational domain and arrangement of the rough elements</figcaption></figure>

sults. More recently, Koo and Kleinstreuer [9] modeled the near-wall 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 µm to 50 µm in 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 rough-wall 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].

## 2 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.

### 2.1 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[\mu m]$): 0.1, 0.2, 0.3, 0.5, 0.7, 1, 2

Roughness side ($d[\mu m]$): 0.5, 1, 1.5

Roughness spacing ($L_r=l_r[\mu m]$): 2, 3, 4

Channel height ($H[\mu 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) \quad (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.

$$\text{at } z=0 \text{ and } z=l_r \quad dU/dz=0, \quad dV/dz=0, \quad W=0 \quad (2)$$

The condition of symmetry at the top plane of the computational domain is written

$$\text{at } y=H/2 \quad dU/dy=0, \quad V=0, \quad dW/dy=0 \quad (3)$$

The system of boundary conditions is completed by the no-slip condition on all the solid boundaries.

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

Continuity equation

$$\nabla U = 0 \quad (4)$$

Momentum equation

$$\rho U \nabla U = - \nabla P + \mu \nabla^2 U \quad (5)$$

### 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 SIMPLEC (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.

### 2.5 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
<figure><img src="image_3.png" /><figcaption>Fig. 2 Sketch of a control volume over the rough surface. Side view.</figcaption></figure>

$$T_3 = T_1 + \frac{T_1 - T_2}{(h_2/h_1)^2 - 1} \quad (6)$$

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.

## 3 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 Equations.** The present section is devoted to the derivation of the momentum equation in the roughness layer. Using the method of Taylor et al., conservation equations for mass and momentum are derived for control volumes (CV) of infinitesimal thickness $\delta y$ in the direction normal to the wall, extending over a large distance compared to the rough elements spacing in the other directions (respectively, $\delta x$, $\delta z \gg L_r$, $l_r$, Fig. 2). Flow quantities are averaged over the lengthscales $\delta x$, $\delta z$ as in homogenization methods.

**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

$$U(y) = \frac{1}{L_r l_r - d^2} \int \int_{dS(x)} u(\xi, y, \zeta) d\xi d\zeta \quad \text{for } 0 \le y \le k \quad (7)$$

where $u(\xi, y, \zeta)$ stands for the actual flow velocity in the roughness region ($y \le k$). Integration is performed over the surface $dS(x)$ open for flow in an elementary cell of sides $l_r \times L_r$ (Fig. 3). This surface is independent of $y$ in the case of cylindrical rough elements.

The local blockage coefficient in the $x$ direction is the fraction of the area open for flow through a slice of thickness $\delta y$ in a plane $x=Cst$. This fraction is obviously a function of $x$ in the roughness layer. Like Taylor et al. [12], we consider the averaged value of the local blockage coefficient over a large distance $\delta x$. For the particular case of periodically distributed rough elements, the averaged blockage coefficient $\beta$ is given by

<figure><img src="image_4.png" /><figcaption>Fig. 3 Surface open for flow. Top view.</figcaption></figure>

$$\beta = \frac{1}{L_r l_r} \int \int_{dS(x)} d\xi d\zeta \quad (8)$$

and is therefore

$$\begin{align}
\beta &= 1 - \frac{d^2}{L_r l_r} \quad \text{for } 0 \le y \le k \tag{9} \\
\beta &= 1 \quad \text{for } k \le y \nonumber
\end{align}$$

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, $\beta$ is independent of $y$. The maximum blockage effect corresponds to $\beta=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) = \frac{1}{L_r l_r} \int \int_{dS(x)} u(\xi, y, \zeta) d\xi d\zeta \quad (10)$$

so that

$$U(y) = u_D(y)/\beta \quad (11)$$

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.

$$\bar{p}(x,y) = \frac{1}{L_r l_r - d^2} \int \int_{dS(x)} p(\xi, y, \zeta) d\xi d\zeta \quad \text{for } 0 \le y \le k \quad (12)$$

We suppose that $\bar{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).

**3.1.2 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<sub>d</sub>

$$
\delta F_d = \frac{1}{2} \rho U^2(y) C_d d \delta y \quad (13)
$$

Finally, the momentum equation is

$$
0 = -\frac{d}{dx}(\beta P) + \frac{d}{dy}\left(\mu\beta\frac{dU}{dy}\right) - \frac{1}{2}\rho C_d d\frac{U^2}{L_r l_r} \quad \text{for } 0 \le y \le k \quad (14)
$$

$$
0 = -\frac{dP}{dx} + \mu \frac{d^2 U}{dy^2} \quad \text{for } k \le y \le H/2
$$

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.

**3.1.3 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 $\beta < 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 dis-
continuous like the blockage coefficient at the interface.

$$
\beta U(y = k^{-}) = U(y = k^{+}) \quad (15)
$$

As mentioned before, this discontinuity would not exist for rough-
ness 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 intro-
ducing 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 $\delta y$ and
lengths $\delta x$, $\delta z$ in the other directions, centered at $y=k$. The force
balance on this CV requires

$$
\left[
\begin{aligned}
& \mu \frac{dU}{dy} \Big|_{y=k+\delta y/2} - (1-\beta)\mu \frac{dU}{dy} \Big|_{y=k^{+}} - \beta\mu \frac{dU}{dy} \Big|_{y=k-\delta y/2} \Big] \delta x \delta z \\
& + \left[ dP \frac{(1+\beta)}{2} \delta z - \frac{1}{2} \rho U^2(k) C_d \frac{d\delta x \delta z}{2L_r l_r} \right] \delta y = 0
\end{aligned}
\right]
\quad (16)
$$

When $\delta y$ tends to zero, the pressure force and the drag force
become negligibly small so that Eq. (16) reduces to

$$
\mu \frac{dU}{dy} \Big|_{y=k^{-}} = \mu \frac{dU}{dy} \Big|_{y=k^{+}} \quad (17)
$$

or using the Darcy velocity

$$
\frac{1}{\beta} \mu \frac{du_D}{dy} \Big|_{y=k^{-}} = \mu \frac{du_D}{dy} \Big|_{y=k^{+}} \quad (18)
$$

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.

$$
U(y=0) = 0 \tag{19}
$$

$$
\left. \frac{\partial U}{\partial y} \right|_{y=H/2} = 0 \qquad (20)
$$

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

$$
U(y=H)=0 \qquad (21)
$$

**3.1.4 Drag Coefficient.** The model relates the drag force $\delta F_d(y)$ exerted by the fluid on a slice of rough element of height $\delta y$ to the local velocity $U(y)$. The starting point of the analysis considers a single rough element oriented across a fluid stream of uniform velocity $U_\infty = U(y)$. This flow is mainly characterized by the local Reynolds number $Re_d = \rho U(y) d/\mu$, which is based on a typical dimension of a block cross section, namely, the side length $d$. $Re_d$ is small in most microchannel flows: $Re_d \approx 5$ at the top of the rough elements for a typically high value of $Re_{D_h} = 2000$ in a channel of relative roughness 2.5% and for elements of relative width ($d/H$) about 2.5%. The flow is therefore dominated by viscous effects in the vicinity of the blocks.

Dimensional analysis shows that the drag force may be written

$$
\delta F_d(y) = \frac{1}{2} K \left( \frac{k}{d'}, \frac{y}{k'}, \text{Re}_d \right) \left( \mu \frac{U(y)}{d} \right) d \delta y \quad (22)
$$

$\delta F_d$ is also related to the drag coefficient $C_d$ by Eq. (13) so that $K=C_d\text{Re}_d$. Eliminating $C_d$ by using Eqs. (13) and (22), the last term of Eq. (14) becomes $1/2K\mu U/L_r l_r$.

It is well known [14] that inertia effects cannot be totally dis-
carded for creeping flows near extremely long cylinders. How-
ever, 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<sub>Dh</sub>=200. Conse-
quently 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 three-
dimensional nature of the low-Reynolds number flow near the
blocks. For creeping flows normal to the major axis of very slen-
der prolate spheroids (E>1), the drag force is given by [15]

$$
D = \frac{4\pi\mu U 2a}{\frac{1}{2} + Ln \frac{2a}{b}} \quad (23)
$$

where $E=a/b$, $2a$=major axis, $2b$=minor axis of the spheroid.

Following the form of Eq. (23), we propose for K the following
expression:

$$
K = \frac{K_1}{1 + K_2 L n \frac{k}{d}} \quad \text{for } \frac{k}{d} \ge 1 \qquad (24)
$$

$$
K = K_1 \quad \text{for } \frac{k}{d} \le 1
$$

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₁ 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<sub>d</sub> 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<sub>1</sub> is nearly independent of y in the upper part of
the blocks, as presented later. For $k/d \le 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

$$U^* = \frac{U}{-\frac{1}{\mu} \frac{dP}{dx} \frac{H^2}{12}} \quad (25)$$

Modeling $K$ as independent of $y$, straightforward computations lead to the following equations for the velocity profile in a channel with two rough sides:

for $0 \le y^* \le k^*$

$$U^*(y^*) = \frac{12}{A^*} \left\{ 1 + \frac{\sqrt{A^*} \left( \frac{1}{2} - k^* \right) \sinh(\sqrt{A^*} y^*) - \cosh[\sqrt{A^*} (k^* - y^*)]}{C_1} \right\} \quad (26)$$

for $k^* \le y^* \le \frac{1}{2}$

$$U^*(y^*) = \frac{12\beta}{A^*} \left[ 1 + \frac{\sqrt{A^*} \left( \frac{1}{2} - k^* \right) S_1 - 1}{C_1} \right] + 6[y^* - y^{*2} + k^*^2 - k^*]$$

where $A^*=1/2K/L_r^*l_r^*\beta$, $C_1=\cosh(\sqrt{A^*}k^*)$, $S_1=\sinh(\sqrt{A^*}k^*)$. For a semirough channel, the expressions are

For $0 \le y^* \le k^*$

$$U^*(y^*) = \frac{12}{A^*} (1 - e^{-\sqrt{A^*} y^*}) + 2D_1 \sinh(\sqrt{A^*} y^*) \quad (27)$$

For $k^* \le y^* \le 1$

$$U^*(y^*) = (y^* - 1) \left[ 2D_1 \sqrt{A^*} C_1 + 12 \left( k^* + \frac{E_1}{\sqrt{A^*}} \right) \right] + 6(1 - y^{*2})$$

where $E_1=\exp(-\sqrt{A^*}k^*)$

$$D_1 = \frac{6}{\sqrt{A^*} (k^* - 1) C_1 - \beta S_1} \left[ \frac{k^{*2}}{2} - \frac{1}{2} + (1 - k^*) \left( k^* + \frac{E_1}{\sqrt{A^*}} \right) + \beta \left( \frac{1 - E_1}{A^*} \right) \right]$$

Integration of the velocity profile gives the flow rate $Q$ as a function of the pressure gradient, the geometrical parameters $k^*$, $\beta$, and of $A^*$.

## 4 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.

<figure><img src="image_5.png" /><figcaption>Fig. 4 Distribution of the drag law coefficient $K$ along a rough element. $K^*=d^*=L_r^*/2=l_r^*/2=0.2$, $\beta=0.75$, Re=200.</figcaption></figure>

### 4.1 Drag Coefficient of the Roughness Elements.
As can be seen in Fig. 1, the blocks were divided into separate slices of height $\delta 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.

### 4.2 Velocity and Shear Stress Profiles.
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

<figure><img src="image_6.png" /><figcaption>Fig. 5 Velocity profiles. $k^*=0.1$, $d^*=L_r^*/2=l_r^*/2=0.2$, $\beta=0.75$, Re=1; semirough channel, plus indicates numerical computations, thin lines indicate analytical model; fully rough channel, open circle indicate numerical computations, thick lines indicate analytical model; smooth channel, dotted lines indicate Poiseuille flow.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>Fig. 6 Contribution of viscous and drag forces to the pressure gradient. $d' = L_r^*/2 = l_r^*/2 = 0.2$, $\beta = 0.75$. Solid diamonds and dotted lines indicate friction on the bottom wall: $F_{vw}/F_T$, solid circles and dash-dotted lines indicate friction at the top of the rough elements: $F_{vt}/F_T$, solid triangles and dashed lines indicate drag force: $F_D/F_T$. Solid symbols correspond to numerical computations, continuous lines correspond to the analytical model.</figcaption></figure>

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.

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 \times l_r$ (Fig. 1), the force $F_T$ due to the pressure gradient is balanced by the tangential force $F_{vw}$ due to friction on the bottom wall (surface $\beta L_r l_r$), the tangential force $F_{vt}$ due to friction at the top of the rough elements (surface $(1-\beta)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_{vt}$ tends to $1-\beta$ when $k^*$ tends to zero. It remains approximately constant up to $k^* \approx 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_{vt}$ still contributes to the pressure gradient unlike $F_{vw}$ which corresponds to a weakening velocity gradient when $k^*$ is increased. For high rough elements ($k^* > 0.1$), the contribution of $F_{vt}$ 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_{vt}$ is much weaker in the analytical model. This may be due to the discontinuity of $\beta$, 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

<figure><img src="image_8.png" /><figcaption>Fig. 7 Influence of the roughness element height on the channel height reduction. $d' = l_r^*/2 = l_r^*/2 = 0.2$, $\beta = 0.75$. Open circles indicate numerical computations, thick lines indicate analytical model, dotted lines indicate HWL's results, dashed-dotted lines indicate $1-H_r/H=2k^*$.</figcaption></figure>

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.

### 4.3 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

$$P_0 = - \frac{\frac{dP}{\mu} 2H^3}{Q} \qquad (28)$$

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

$$24 = - \frac{\frac{dP}{\mu} 2H_r^3}{Q} \qquad (29)$$

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.

$$\frac{H_r}{H} = \left( \frac{24}{P_0} \right)^{1/3} \qquad (30)$$

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.

#### 4.3.1 Influence of the Roughness Element Height.
The apparent channel height reduction is moderate up to $k^* \approx 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.
<figure><img src="image_9.png" /><figcaption>Fig. 8 Influence of the roughness element side length on the channel height reduction. $k^*=L_r^*/2=l_r^*/2=0.2$. Dashed lines indicate $\beta$. Open circles indicate Numerical computations, analytical model, dashed-dotted lines indicate $K_2=0$, thick line indicate $K_2=0.5$, dotted lines indicate HWL's results.</figcaption></figure>

### 4.3.2 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 $\beta$ 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.

### 4.3.3 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.

## 4.4 Comparison With Experiments.
In parallel with the analytical and numerical approaches, experiments were conducted in water flows through rough rectangular microchannels [7,8]. In these experiments, the channel walls consisted of parallel bronze blocks separated by a foil of thickness $e_f$ (0.1 mm $<e_f<1$ mm) 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}$

<figure><img src="image_10.png" /><figcaption>Fig. 9 Influence of the roughness element spacing on the channel height reduction. $k^*=d^*=0.2$. Same symbols as in Fig. 8.</figcaption></figure>

$\approx 2$ µm), together with small SiC particles (5 to 7 µm 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

<figure><img src="image_11.png" /><figcaption>Fig. 10 Poiseuille number. Comparison of numerical computations and analytical model with experimental results. ♦ Experimental results; ○ Numerical computations, $k=5$ µm; analytical model, dotted lines indicate $k=6.2$ µm, thick lines indicate 7.2 µm, dashed lines indicate 8.2 µm, dash-dotted lines indicate $Po=24$.</figcaption></figure>
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.

## 5 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 periodic 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=I_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 $\beta$). 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 $\beta$ 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.

## Acknowledgment

This research was supported by the CNRS. Roland Bavière and Gabriel Gamrat were supported by the French Ministry of Education which is gratefully acknowledged. The authors are also very grateful to Eric Chainet for the roughening of the microchannel walls.

## Nomenclature

$C_d$ = drag coefficient, dimensionless

$d$ = width of roughness elements, m

$D_h$ = hydraulic diameter, m

$F_D$ = drag force on roughness elements, N

$H$ = channel height, m

$K$ = drag law coefficient, dimensionless

$k$ = height of roughness elements, m

$L_r$ = longitudinal distance between two roughness elements, m

$l_r$ = transverse distance between two roughness elements, m

$P$ = pressure, pa

$P_o$ = Poiseuille number, dimensionless

$Q$ = volumetric flow rate per unit length, $m^3 s^{-1} m^{-1}$

$Re_d$ = Reynolds number based on local effective velocity and dimension of rough elements, dimensionless

$Re_{D_h}$ = Reynolds number based on bulk velocity and hydraulic diameter, dimensionless

$U$ = spatial average of longitudinal velocity, $m s^{-1}$

$u_D$ = filtration or Darcy velocity, $m s^{-1}$

$V$ = spatial average normal velocity, $m s^{-1}$

$x$ = abscissa along the channel, m

$y$ = coordinate normal to the wall, m

$z$ = transverse coordinate, m

## Greek symbols

$\rho$ = fluid density, $kg m^{-3}$

$\tau_w$ = wall shear stress, $N m^{-2}$

$\mu$ = dynamic viscosity, $kg m^{-1} s^{-1}$

$\nu$ = kinematic viscosity, $m^2 s^{-1}$

## References

[1] Mala, G. M., and Li, D., 1999, “Flow Characteristics of Water in Microtubes,” Int. J. Heat Fluid Flow, 20, pp. 142–148.

[2] Li, Z. X., Du, D. X., and Guo, Z. Y., 2003, “Experimental Study on Flow Characteristics of Liquid in Circular Microtubes,” Microscale Thermophys. Eng., 7, pp. 253–265.

[3] Guo, Z. Y., and Li, Z., 2003, “Size Effect on Microscale Single-Phase Flow and Heat Transfer,” Int. J. Heat Mass Transfer, 46, pp. 149–159.

[4] Judy, J., Maynes, D., and Webb, B. W., 2002, “Characterization of Frictional Pressure Drop for Liquid Flows Through Microchannels,” Int. J. Heat Mass Transfer, 45, pp. 3477–3489

[5] Qu, W., Mala, G. M., and Li, D., 2000, “Pressure-Driven Water Flows in Trapezoidal Silicon Microchannels,” Int. J. Heat Mass Transfer, 43, pp. 353–364.

[6] Pfund, D., Rector, D., Shekarriz, A., Popescu, A., and Welty, J., 2000, “Pressure Drop Measurements in a Microchannel,” AIChE J., 46, pp. 1496–1507.

[7] Bavière, R., Ayela, F., Le Person, S., and Favre-Marinet, M., 2004, “An Experimental Study of Water Flow in Smooth and Rectangular Microchannels,” Second International Conference on Microchannels and Minichannels, Rochester, June 2004, S. G. Kandlikar ed., ASME, New York, pp. 221–228.

[8] Bavière, R., “Etude de l’Hydrodynamique et des Transferts de Chaleur dans des Microcanaux,” Ph.D. thesis, UJF Grenoble University, http://tel.ccsd.cnrs.fr/.

[9] Koo, J., and Kleinstreuer, C., 2003, “Liquid Flow in Microchannels: Experimental Observations and Computational Analyses of Microfluidic Effects,” J. Micromech. Microeng., 13, pp. 578–579.

[10] Kleinstreuer, C., and Koo, J., 2004, “Computational Analysis of Wall Roughness Effects for Liquid Flow in Micro-Conduits,” J. Fluids Eng., 126, pp. 1–9.

[11] Hu, H., Werner, C., and Li, D., 2003, “Influence of Three-Dimensional Roughness on Pressure-Driven Flow Through Microchannels,” J. Fluids Eng., 125, pp. 871–879.

[12] Taylor, R. P., Coleman, H. W., and Hodge, B. K., 1985, “Prediction of Turbulent Rough-Wall Skin Friction Using a Discrete Element Approach,” J. Fluids Eng., 107, pp. 251–257.

[13] Roache, P. J., 1997, “Quantification of Uncertainty in Computational Fluid Dynamics,” Annu. Rev. Fluid Mech., 29, pp. 123–160.

[14] White, F. M., 1991, Viscous Fluid Flow, McGraw-Hill, New York.

[15] Sherman, F. S., 1990, Viscous Flow, McGraw-Hill, New York.