# Modelling Of Roughness Effect On Laminar Flow And Heat Transfer In Rectangular Microchannels

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

► To cite this version:

Gabriel Gamrat, Michel Favre-Marinet, Stéphane Le Person, Roland Bavière, Frédéric Ayela. Mod-
elling Of Roughness Effect On Laminar Flow And Heat Transfer In Rectangular Microchannels. 13th
International Heat Transfer Conference, Aug 2006, sydney, United States. 10.1615/IHTC13.p14.80.
hal-00204895

HAL Id: hal-00204895

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

Submitted on 27 Mar 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 ROUGHNESS EFFECT ON LAMINAR FLOW AND HEAT TRANSFER IN RECTANGULAR MICROCHANNELS

G. Gamrat<fnref n="1" />, M. Favre-Marinet<fnref n="1" />, S. Le Person<fnref n="1" />, R. Baviere<fnref n="1" />,<fnref n="2" />, F. Ayela<fnref n="2" />

## Abstract

The present paper is devoted to the numerical modeling of roughness effects on laminar flow in microchannels. 3D numerical simulations were initially performed and, based on their solutions, a roughness model is proposed. This one-dimensional model is inspired by a discrete-element approach due to Taylor et al. (1985, 1989). In this model, the channel consists in a clear medium adjacent to a porous medium layer where the interactions between the stream and the roughness elements are directly computed by the discrete-element method. Empirical correlations based on 2D numerical simulations of cross-flow through a bank of rods are used for modeling the local drag coefficient and the local Nusselt number of the roughness elements. The appropriate determination of these important parameters allows solving the momentum and energy equations in order to obtain the velocity and temperature profiles. The 1D model solutions are then compared with the 3D numerical solutions showing a very good agreement. The present study shows that roughness significantly influences both the Poiseuille number and the global Nusselt number. The relative increase of the pressure drop is found to be much faster than that of the heat transfer coefficient when the roughness height is increased. Two microchannels were produced with the same three-dimensional roughness arrangement as in the 3D numerical model and experiments were performed with deionized water as the working fluid. The height of microchannels was 106 µm and 153 µm with relative roughness of 16% and 14% respectively. Comparison between experiments conducted in isothermal conditions and model solutions shows the good ability of the numerical models to predict the pressure drop in the rough microchannels, which were tested.

## Nomenclature

<table><tr><td>Cd</td><td>drag coefficient</td></tr><tr><td>C<sub>r</sub></td><td>resistance factor</td></tr><tr><td>Da</td><td>Darcy number</td></tr><tr><td>H</td><td>channel height</td><td>m</td></tr><tr><td>d</td><td>roughness element size</td><td>m</td></tr><tr><td>K</td><td>permeability</td><td>m<sup>2</sup></td></tr><tr><td>k</td><td>roughness element height</td><td>m</td></tr><tr><td>k<sub>f</sub></td><td>thermal conductivity for fluid</td><td>W/mK</td></tr><tr><td>k<sub>s</sub></td><td>thermal conductivity for solid</td><td>W/mK</td></tr><tr><td>L<sub>r</sub></td><td>roughness element spacing in x direction</td><td>m</td></tr><tr><td>l<sub>r</sub></td><td>roughness element spacing in y direction</td><td>m</td></tr><tr><td>Nu</td><td>global Nusselt number (Eq. 21)</td></tr><tr><td>Nu<sub>D</sub></td><td>porous Nusselt number (=hK<sup>0.5</sup>/k<sub>f</sub>)</td></tr><tr><td>Nu<sub>d</sub></td><td>roughness Nusselt number (=hd/k<sub>f</sub>)</td></tr><tr><td>Po</td><td>Poiseuille number</td></tr></table>

<table><tr><td>Re<sub>D</sub></td><td>porous Reynolds number<br>(=u<sub>D</sub>K<sub>max</sub><sup>0.5</sup>/v)</td></tr><tr><td>Re<sub>d</sub></td><td>roughness Reynolds number (=u<sub>D</sub>d/v)</td></tr><tr><td>β</td><td>porosity</td></tr><tr><td>θ</td><td>dimensionless temperature</td></tr></table>

### subscripts

<table><tr><td>a</td><td>axial</td></tr><tr><td>b</td><td>bulk</td></tr><tr><td>D</td><td>Darcy</td></tr><tr><td>f</td><td>fluid, friction</td></tr><tr><td>max</td><td>maximal</td></tr><tr><td>p</td><td>pressure</td></tr><tr><td>s</td><td>solid</td></tr><tr><td>top</td><td>top surface</td></tr><tr><td>w</td><td>wall</td></tr></table>

<footnote n="1">¹ Laboratoire des Ecoulements Géophysiques et Industriels, CNRS-UJF-INPG, 1025 rue de la Piscine, BP 53 X, 38041 Grenoble Cedex, France</footnote>

<footnote n="2">² Centre de Recherches sur les Très Basses Températures, CNRS B.P. 166, 34042 Grenoble Cedex 09, France</footnote>
# 1. Introduction

Recent advances in ever faster and smaller electronic equipments have entailed the development of efficient cooling systems. This has stimulated a strong current of research on flows in microchannels. However many published papers present a strong dispersion of results and no definitive tendency is observed concerning the friction factor and the heat transfer coefficient in microflows. Measurements in microchannels are very difficult and require careful accounting for phenomena that are often negligible in standard flows. Recent experimental works accurately evaluated the uncertainty of measurements and their results are now consistent with the conventional theory. For example, Kohl et al. (2004) experimentally investigated laminar water flow through smooth microchannels with hydraulic diameter ranging from 25 to 100 µm. They used internal pressure measurements to avoid entrance effects. The average value of their results exhibits an increase by only 4% above the theoretical value, whereas the uncertainty is 11%. The authors did not observe an early transition to turbulence. They concluded that the large inconsistencies in previously published data are probably due to instrumentation errors. The same conclusion was drawn by Baviere et al. (2005) from experiments in smooth microchannels ranging from 4.5 to 21 µm.

Some researchers proposed that a possible explanation of the friction factor increase above conventional values could be an effect of Electrical Double Layer. However, a recent work of Phares and Smedley (2004) shows that electroviscous effects are not observed in 150 µm microtubes. They examined also the effect of surface roughness and they concluded that it plays an important role in laminar flow through microtubes. They reported an increase in Poiseuille number values by up to 17% for a relative surface roughness equal to 2.5%. No measurable effect of viscosity was detected for three mixtures of different viscosities.

Judy et al. (2002) performed extensive measurements of pressure losses for laminar flow through different microtubes (two different materials and two different shapes) in the range 15 –150 µm. After a careful examination of possible errors they reported that the uncertainties in the Poiseuille number are as high as 10 –15% for smooth fused silica and 20 – 21.5% for rough steel microtubes. They did not observe significant deviation from Stokes theory, however they concluded that the detection of additional phenomena arising in microflows is very difficult since their influence is below experimental errors at microscales.

Recent experimental research suggests that surface roughness is mainly responsible for deviation from conventional theory in microchannels. Guo and Li (2003) stated that for rough wall microchannels, the form drag of rough elements is predominant and could be regarded as one reason leading to increase in the friction factor. Thus there is a strong interest to explore this effect precisely and to develop simple models which could accurately predict the influence of roughness on fluid flow and heat transfer. Following the conclusion of Judy et al. (2002) about the relative importance of microscale effects and experimental uncertainties, numerical simulation seems to be an adequate tool for exploring the influence of a particular phenomenon on microflows. This was the idea behind a recent study of Croce et al. (2005), who investigated numerically the effects of 3D roughness on microchannel heat transfer and pressure drop. They found a 16% increase in the Poiseuille number and a much smaller increase of heat transfer coefficient for a relative roughness height equal to 2.65%. In their study, conically shaped roughness elements were periodically distributed on the walls of a smooth channel. The numerical study of Hu et al. (2003) reports extensive results about the surface roughness impact on pressure drop in microchannels for cubic shaped roughness elements. They proposed an empirical model which accounts for additional pressure drop in terms of a relative channel height reduction. In other words, the apparent pressure drop is accounted for by the conventional value of the Poiseuille number used with a reduced
channel height. However, the proposed formulas suffer from limitations and cannot be used with other shapes of roughness elements.

Koo and Kleinstreuer (2003) proposed that the roughness region could be modeled by an equivalent porous medium layer (PML hereafter). They modeled the additional viscous forces due to rough elements in terms of the Darcy number and they used a nonlinear term to account for inertia forces. By means of these additional source terms introduced into the Navier–Stokes equation, they were able to reproduce the experimental results of Guo and Li (2003). However, in spite of good agreement between PML model and experimental results they did not relate the Darcy number to the geometrical parameters of roughness. In a next paper, Koo and Kleinstreuer (2005) extended their PML model to heat transfer phenomena in rough wall microconduits. They proposed a one-equation model to account for heat transfer increase due to rough elements. However, such a one-equation model seems to be inadequate for cases where convection heat transfer is predominant so that the temperature distribution in rough elements strongly differs from that in the fluid. In such cases the use of two-equations models is advisable as illustrated by Kim and Kim (1999).

Taylor et al. (1985 and 1989) presented a discrete-element method for predicting skin friction and heat transfer in turbulent flows. This approach is based on a volume averaging technique and considers separately the form drag due to roughness elements and the shear drag at the flat part of the surface. The information about the roughness element shape and size is directly introduced into the momentum equation by the drag coefficient when it is correctly estimated. This approach was used in the work of Baviere et al. (2005), where the resistance coefficient due to rough elements was based on the formula for the drag force on very slender prolate spheroids in creeping flow. The additional constants were adjusted to fit the results of the analytical model to those of numerical simulations.

The current study is devoted to the investigation of the surface roughness effect on fluid flow and heat transfer in microchannels by means of numerical modelling. The geometrical model of roughness follows the details of the model developed by Hu et al. (2003). The paper presents successively 3D numerical simulations and the PML model. The latter one is based on an extended Darcy equation for fluid flow and two equations for heat transfer. The control factors, i.e. the permeability and the local Nusselt number, are determined by means of the numerical simulations and are related to the geometrical parameters. It is shown that 2D simulations, relative to cross-flow through an array of rods, are sufficient to determine the drag force and the local Nusselt number on the rough elements.

## 2. Numerical model

### 2.1. Geometrical model of roughness, computational domain and boundary conditions

The surface roughness is represented by blocks distributed on the smooth walls of a plane microchannel. The simplified present model assumes that the roughness elements are periodically distributed either in aligned or staggered arrangements as presented in Figure 1. The roughness elements are parallelepipeds of square cross-section of side length *d* and height *k*. The wall region is then characterized by a constant porosity $\beta = 1 - d^2/L_r l_r$. The choice of a regular pattern of roughness elements allows using periodical boundary conditions and reduces the computational domain to one period as depicted by dashed lines in Figure 1. The other advantage of periodicity is the possibility to use volume averaging in order to establish the governing equations for the PML model.

Figure 1. Arrangements of rough elements; left: aligned, right: staggered.

<figure><img src="image_2.png" /></figure>
<figure><img src="image_3.png" /><figcaption>Figure 2. Computational domain</figcaption></figure>

Since the PML model principles are the same for staggered and aligned arrangements, the presentation is restricted to the latter one.

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 (Figure 2). The longitudinal and transverse dimensions are thus equal to the rough element spacing. The microchannel wall consists of a block of thickness Hs, heated at uniform heat flux q'' on the external side. The computational domain was treated as the central part of an extremely long channel so that the flow was considered as fully developed in the longitudinal direction and the heat flux distortion due to axial conduction at the channel ends was assumed to be negligible. As a result, the flow properties are periodic in the x direction. This assumption allows applying periodic boundary conditions on the opposite sides of the domain in the longitudinal direction. The periodic boundary condition for any flux may be written:

$$ \vec{\Theta}(x, y, z) = \vec{\Theta}(x + L_r, y, z) \qquad (1) $$

Due to periodicity, state variables as pressure and temperature can be written as the sum of a linear gradient and a periodic component:

$$ \phi(x, y, z) = \frac{d\phi}{dx} x + \tilde{\phi}(x, y, z) \qquad (2) $$

The term $dp/dx$ relates only to the fluid domain and is a priori known in our computations while the term $dT/dx$ is determined by the energy conservation in the fluid and solid domains. It was then assumed that $dT/dx$ in the solid and fluid phases is equal to $q''_0 l_r / \dot{M} c_p$ where $\dot{M}$ is the mass flow rate. Since $dT/dx$ equally contributes to the temperature field in both domains, we can introduce the dimensionless temperature $\theta(x, y, z) = (\tilde{T}(x, y, z) - \tilde{T}_w(x)) / (q''_0 2H / k_f)$ which is periodic:

$$ \theta(x, y, z) = \theta(x + L_r, y, z). \qquad (3) $$

$T_w(x)$ is the mean wall temperature at $y=0$. The boundary condition as given by equation (1) relates to the fluid and solid phases. For the fluid domain, $\vec{\Theta}$ stands for the velocity vector while in the solid block it stands for the heat flux in axial direction $\vec{q}_a$'. The axial heat flux in the solid block results from the constant temperature gradient $dT/dx$ and is then expressed by

$$ \vec{q}_a'' = k_s \frac{q''_0 l_r}{\dot{M} c_p} \vec{i} \qquad (4) $$

where $\vec{i}$ is the unit vector in x direction.

For a large conductivity ratio $k_s / k_f$, the axial heat flux in the fluid was assumed to be negligible.

For the sake of symmetry, the velocity and heat flux gradient normal to the lateral sides and top plane of the computational domain were assumed to be zero. The symmetrical boundary conditions are presented in the form: $\nabla\Theta\cdot\vec{n}=0$, $\vec{\Theta}\cdot\vec{n} = 0$, $\nabla\phi\cdot\vec{n} = 0$ (5)

where $\vec{n}$ is the vector in the direction normal to the symmetry surface.
The system of boundary conditions is completed by the no-slip velocity condition and continuity of the temperature at all fluid/solid interfaces. The flow is characterized by the bulk velocity $u_b$. The Reynolds number is defined by $Re = u_b 2H/v$ and the dimensionless pressure drop is represented by

the Poiseuille number $Po = -\frac{1}{\mu} \frac{dp}{dx} \frac{2H^2}{u_b}$. The Prandtl number is fixed to $Pr = 7$.

## 2.2. Equations

Modelling the flow as laminar, steady and incompressible, the governing equations consist of the following system.

Continuity equation: $$ \nabla \vec{u} = 0 \qquad (6) $$

Momentum equation: $$ \rho(\vec{u} \cdot \nabla \vec{u}) = -\nabla p + \mu \nabla^2 \vec{u} \qquad (7) $$

A further simplification assumes that viscous dissipation can be neglected in the energy equation for the fluid phase which is written as follows:

$$ \rho c_p (\vec{u} \cdot \nabla T) = k_f \nabla^2 T \qquad (8) $$

and for the solid phase:

$$ \nabla^2 T = 0 \qquad (9) $$

## 2.3. 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 non-linear, 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 normalized residuals decreased below $10^{-12}$. Computations were performed on orthogonal grids generated by Gambit 2.1.2.

## 2.4. 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 x 40 x 32, 48 x 60 x 48, 64 x 80 x 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 intermediate and the fine grids was 0.7%. Furthermore the estimation of accuracy by means of Richardson extrapolation was performed. Using the Richardson extrapolation in form:

$$ T_3 = T_1 + (T_1 - T_2) / ((h_2 / h_1)^2 - 1) \qquad (10) $$

and disposing two approximations $T_1$ and $T_2$ for two mesh size $h_1$ (fine grid) and $h_2$ (coarse grid) respectively one can estimate a third approximation $T_3$ whose principal term of error is higher order than that for $T_1$ or $T_2$. The difference between the pressure gradient obtained for the mesh 48 x 60 x 48 and that obtained from the Richardson extrapolation was lower than 2% and this mesh was adopted. The theoretical values for the fully developed laminar flow, $Po_{th}=24$ and $Nu_{th}=8.235$ were recovered with excellent accuracy.

# 3. Porous medium layer model

As mentioned before, Koo and Kleinstreuer (2003) determined the flow resistance (Darcy number) by using the experimental results of Guo and Li (2003). Because of the lack of information about
the distribution and shape of the experimental roughness, they could not relate the flow resistance factor to geometrical parameters. The present PML model considers the same geometrical model of roughness as the numerical analysis. The current section describes the model and its extension to heat transfer phenomena.

## 3.1. Momentum equation

The present model combines ideas from the PML model of Koo and Kleinstreuer (2003) and from a discrete-element approach proposed by Taylor et al. (1985) to compute the rough-wall skin friction in turbulent flows. The approach considers control volumes (CV) of infinitesimal thickness in the direction normal to the wall. It directly accounts for the drag forces exerted on the roughness elements by means of volume averaging inside CV as depicted in Figure 2. The momentum equation in the roughness region results from the balance of forces inside CV

$$ -\frac{dp}{dx} = (p_w - p_l) \frac{d}{L_r l_r} + \tau_s \frac{2d}{L_r l_r} - \mu \frac{d^2 u_D}{dy^2} \quad (11) $$

where $τ_s$, $p_w$ and $p_l$ are respectively the shear stress, the windward and leeward pressures acting on the portion $CV_s$ of the rough element included in CV. The model uses the Darcy velocity $u_D(y)$ defined by

$$ u_D(y) = \frac{1}{L_r l_r} \int_{CV} u(x, y, z) dxdz \quad (12) $$

Applying the state variable decomposition as given by equation (2) allows rewriting the momentum equation (11) in form:

$$ \mu \frac{d^2 u_D}{dy^2} = \frac{dp}{dx} \beta + \frac{1}{2} \rho u_D^2 (Cd_p + Cd_f) \frac{d}{L_r l_r} \quad (13) $$

where $Cd_p = \frac{\tilde{p}_w - \tilde{p}_l}{0.5 \rho u_D^2}$ and $Cd_f = \frac{2 \tau_s}{0.5 \rho u_D^2}$ are the pressure and friction drag coefficients respectively. Their sum represents the total drag coefficient $Cd$.

## 3.2. Energy equations

As for hydrodynamics, the volume averaging technique is applied in order to derive the energy equation. A two-equation model is applied in order to account for the convective heat exchange between the rough element surface and the surrounding fluid. The energy balance inside CV leads to the following energy equation for the fluid phase:

$$ - \rho c_p u_D \frac{dT}{dx} + k_s (1 - \beta) \frac{d^2 T_s}{dy^2} + k_f \beta \frac{d^2 T_f}{dy^2} = 0 \quad (14) $$

and is complemented by the energy equation inside the rough element $CV_s$

$$ k_s (1 - \beta) \frac{d^2 T_s}{dy^2} = q'' \frac{4d}{L_r l_r} \quad (15) $$

where $q''(y)$ is the heat flux exchanged by $CV_s$ with the stream. The temperatures $T_s$ and $T_f$ are the effective solid and fluid temperatures averaged respectively over the area occupied by the solid and the area open for flow.

$$ T_s(y) = \frac{1}{d^2} \int_{CV_s} \tilde{T}_s(x, y, z) dxdz \quad \text{and} \quad T_f(y) = \frac{1}{L_r l_r - d^2} \int_{CV - CV_s} \tilde{T}_f(x, y, z) dxdz \quad (16) $$

The products $(1 - \beta)T_s = T_{sD}$ and $\beta T_f = T_{fD}$ represent the average temperatures over the total area of CV and their sum can be regarded as the Darcy temperature $T_D(y) = T_{sD}(y) + T_{fD}(y)$. $T_D$ has to
be used in order to assure the continuity of average temperature at rough/clear interface. Substituting the left-hand term of equation (15) into equation (14) and eliminating $T_f$ results in:

$$k_f \frac{d^2 T_D}{dy^2} = \rho c_p u_D \frac{q_0'' l_r}{\dot{M} c_p} - Nu_d \frac{4}{L_r l_r} \left(1 - \frac{k_f}{k_s}\right) (T_s - T_D) k_f \quad (17)$$

where $Nu_d(y) = \frac{q_0'' d}{(T_s - T_D)k_f}$ is the local Nusselt number in the sense of the PML model, in fact the global Nusselt number for a slice of rough element $CV_s$.

### 3.3. Drag coefficient modelling

In order to resolve equation (13) it is necessary to model the pressure and friction drag coefficients. The 3D numerical simulations were used in order to yield all the necessary data. Figure 3 presents distributions of $Cd(y)$ as a function of the local Reynolds number $Re_d(y) = u_D(y)d/v$ for several values of the geometrical parameters. It clearly shows that the drag coefficient $Cd$ is inversely proportional to $Re_d$. The proportionality factor $C_r$ (resistance coefficient hereafter) varies neither with the mass flow rate nor with the relative roughness height $k^*=k/0.5H$. This advantage results from applying the pressure decomposition as given by equation (2), after which the pressure drag coefficient $Cd_p$ accounts only for forces due to the periodic component of pressure. A slight increase of the resistance coefficient $C_r$ at low $Re_d$ number is observed whatever the roughness considered. It could be explained by inadequately predicted pressure forces at the vicinity of the smooth bottom wall where velocity and pressure forces tend to zero. In fact, pressure forces are overestimated by the pressure interpolation scheme which leads to variations of $C_r$ when approaching the smooth bottom wall. This uncertainty in pressure forces estimation is however of weak importance for the global results because it occurs in a region where the contribution of drag forces is very small like the velocity. It is then justified to assume that the coefficient $C_r$ is constant along the roughness element height. Figure 3 shows that $C_r$ only depends on the porosity $\beta$ and increases when $\beta$ is decreased. Thus, it suggests to model the coefficient $C_r$ by means of 2D numerical simulations instead of computationally expansive 3D simulations. For 2D cross-flow through an array of rods and after applying the relation $C_r = CdRe_d$ the momentum equation (13) is reduced to

$$0 = \frac{dp}{dx} \beta + \frac{\mu}{K} u_D \quad (18)$$

where $K = 2L_r l_r / C_r$ is the permeability (in the numerical model, $L_r = l_r = 2\mu m$). Equation (18) is similar in form to the conventional Darcy equation, except that the porosity $\beta$ appears after applying the pressure decomposition.

<figure><img src="image_4.png" /><figcaption>Figure 3. Distribution of drag coefficient along rough element</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Figure 4. Variation of permeability with the porous Reynolds number for several values of porosity</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>Figure 5. Porous Nusselt number vs. porous Reynolds number for several values of porosity</figcaption></figure>

Figure 4 shows the variation of *K* as a function of the Reynolds number *Re*<sub>D</sub> for several values of porosity. The numerical results are approximated by $K = a + bRe_D$. The values of $K_{max} = a$ were used to define a length scale in the definition of *Re*<sub>D</sub> number $Re_D = u_D\sqrt{K_{max}}/\nu$. Regardless the porosity, a slight variation of *K* with *Re*<sub>D</sub> appears when *Re*<sub>D</sub> exceeds a value of order O(10). In fact, laminar microchannel flows are characterized by very small values of the permeability, so that *Re*<sub>D</sub> does not reach this high-range for moderate roughness height. Thus it is expected that the drag forces due to surface roughness are proportional to the mean velocity. This is consistent with our 3D

simulations and with the study of Croce et al. (2005) where the variations of *Po* with *Re* were negligible. However, the 2D modelling suppresses the influence of recirculation zones induced by the flow upon the rough elements which could be stronger for other roughness geometry, in particular like in flows over 2D ribs as reported by Croce and D'Agaro (2004). In such cases a stronger dependence of *Po* with *Re* is expected.

## 3.4. Nusselt number modelling

The Nusselt number $Nu_d$ occurring in equation (17) is modelled by means of 2D numerical simulations. Since the model uses the temperature average in the sense of Darcy $T_D$, the heat transfer coefficient is defined by $h = q'' / (T_s - T_D)$. This definition misses the usual meaning of the *h* coefficient where the fluid bulk temperature $T_b$ is used instead of $T_D$. This drawback results from the volume averaging technique which does not allow using $T_b$ in the energy equation. This issue will be again discussed about the results.

The 2D simulations already used for the estimation of permeability, were extended to the heat transfer problem. The boundary condition inside solid elements was supposed to be uniform volume source heating. Figure 5 shows the variations of $Nu_D$ as a function of $Re_D$ for the same set of porosity as presented in Figure 4. The permeability $K$ was used as the length scale for $Nu_D$, i.e. $Nu_D = Nu_d \sqrt{K} / d$. Contrary to the permeability, the $Nu_D$ number does not exhibit a continuous drop when the porosity is decreased but one can observe a minimum value of $Nu_D$ for $\beta \approx 0.75$. Comparison of Figure 4 and Figure 5 shows that *h* is more sensitive than *K* to $Re_D$. It may be probably attributed to the value of the Prandtl number $Pr = 7$ used in the current study, which means that the thermal boundary layers developing on the roughness surface are thinner than the corresponding velocity boundary layers.

## 3.5. Boundary conditions

After volume averaging, the problem of fluid flow and heat transfer in rough wall microchannels reduces to a system of second order ordinary differential equations. It can be solved numerically under appropriate boundary conditions. It was assumed that the fluid and solid temperature are equal ($\theta_s = \theta_D$) and the velocity is equal to zero at the bottom wall ($y=0$). Symmetry boundary conditions were used at the channel symmetry plane ($y=0.5H$). The continuity of temperature $T_D$ and Darcy velocity $u_D$ is satisfied in the whole computational domain. The boundary condition at the PML/clear region interface raises a particular problem. The development of velocity and thermal boundary layers on the rough elements top surface contributes to the discontinuity of shear
stress and heat flux at the PML/clear region interface in the volume averaging approach. In order to account for these discontinuities the present model considers a control volume $CV_N$ defined by $k \le y \le k + \delta y$ where the momentum and energy balances are written. Since $K$ and $Nu_D$ are mostly modelled by means of 2D numerical simulations, it was necessary to assume that the developing velocity and thermal boundary layers on the side and top surfaces of roughness elements are similar. As a consequence, the shear stress at the rough elements top surface is deduced from the value of the average friction drag coefficient $Cd_f$ which was found at the side surfaces of rough elements. In the same way, the second boundary condition for the energy equation in the solid phase assumes identical convection heat transfer coefficient $h$ for the top surface and the side surface of rough element. The dimensionless system of governing equations is:

$$ 
\begin{array}{l}
0 = \frac{\beta Po}{8} - \frac{u^*}{Da} + \frac{d^2 u^*}{dy^{*2}} \\
\\
0 \le y^* \le k^* \qquad \frac{d^2 \theta_D}{dy^{*2}} = \frac{u^*}{4} - Nu_D \frac{d^*}{\sqrt{Da}} (\theta_s - \theta_D) \frac{H^{*2}}{L_r l_r} \left(1 - \frac{k_f}{k_s}\right) \\
\qquad \frac{d^2 \theta_s}{dy^{*2}} = Nu_D \frac{d^*}{\sqrt{Da}} (\theta_s - \theta_D) \frac{H^{*2}}{d^{*2}} \frac{k_f}{k_s} \\
\\
k^* \le y^* \le k^* + \delta y^* \qquad 0 = \frac{Po}{8} - \frac{u^*}{Da_f} \frac{d^*}{2\delta y^*} + \frac{d^2 u^*}{dy^{*2}} \\
\qquad \frac{d^2 \theta_D}{dy^{*2}} = \frac{u^*}{4} - Nu_{top} \frac{d^*}{\sqrt{Da}} (\theta_s - \theta_D) \frac{d^*}{4\delta y^*} \frac{H^{*2}}{L_r l_r} \left(1 - \frac{k_f}{k_s}\right) \\
\\
k^* + \delta y^* \le y^* \le 1 \qquad 0 = \frac{Po}{8} + \frac{d^2 u^*}{dy^{*2}} \\
\qquad \frac{d^2 \theta_D}{dy^{*2}} = \frac{u^*}{4}
\end{array}
\qquad (19) $$

and boundary conditions:

$$ y^* = 0 \qquad u^* = 0, \theta_s = \theta_D $$

$$ y^* = k^* \qquad u^*(k^{*-}) = u^*(k^{*+}), \quad \theta_D(k^{*-}) = \theta_D(k^{*+}) $$

$$ y^* = 1 \qquad \frac{du^*}{dy^*} = \frac{d\theta_D}{dy^*} = 0 $$

where: $Da = \frac{K}{(H/2)^2}$, $Da_f = \frac{2L_r l_r}{Cd_f Re_d (H/2)^2}$ and $Nu_{top} = (2.25 - \beta)Nu_D$ are the Darcy number, the friction component of Darcy number and the Nusselt number at the rough element top surface respectively. The lengths are normalized by $H/2$ and the dimensionless velocity is defined as: $u^* = u_D / u_b$. The relation $Nu_{top} = (2.25 - \beta)Nu_D$ results from the hypothesis that the heat transfer coefficient $h$ at the top and side surfaces of rough elements is equal. This empirical relation takes into account the fact that in cases of dense packed rough elements (small $\beta$) convection heat transfer takes place almost only at the side walls. The relation is valid when the convective heat transfer in porous layer is predominant. The system of equations (19) was discretized and solved by means of a first order finite difference method using Matlab.
<figure><img src="image_7.png" /><figcaption>Figure 6. Dimensionless velocity profiles at constant mass flow rate</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Figure 7. Fluid and solid dimensionless temperature profiles. Continuous lines: PML model; symbols: 3D simulations</figcaption></figure>

## 4. Results

Figure 6 shows velocity profiles at constant mass flow rate for three different thicknesses of PML while the Darcy number was kept constant (constant porosity $\beta$ of PML). The velocity profile for $k*= 0.2$ is compared with the 3D numerical solution, showing excellent agreement of the results. As it was expected, the velocity profiles strongly depend on the PML thickness. The apparent velocity gradient immediately above the PML corresponds to the overall pressure gradient.

Figure 7 shows the profiles of the fluid $\theta_f$ and solid $\theta_s$ temperatures defined by equations (16), for the same PML thickness and Darcy number as in Figure 6. As for the velocity profiles, a very good agreement is found between 3D and PML model solutions. The apparent discontinuity of fluid temperature $\theta_f$ results from the discontinuity of porosity $\beta$ and from the strong heat flux at the top surface of rough elements. The observed differences between fluid and solid temperatures indicate that in the case of strong convection the assumption of thermal equilibrium is not appropriate. The rough elements of high thermal conductivity when compared to the conductivity of water ($k_f k_s = 5.1 * 10^{-5}$ in the current study), conduct heat directly towards the region of higher velocities. Since the highest heat sink occurs in the central part of the channel, the conductive heat transfer rate in solid elements increases with the roughness height. The heat rate is mostly conducted through the rough elements causing a small temperature gradient both in the fluid and solid phases because the heat flux through the fluid is small and the thermal conductivity is high in the solid phase. It is then clear why the temperature profiles showed in Figure 7 become flatter when the PML thickness increases. The observed small temperature gradient is in agreement with the assumption reported by Qu et al. (2000). In the current study it leads however to the opposite conclusion, namely that roughness contributes to increase the convective heat transfer so that the modified $Nu$ relationship proposed by Qu et al. (2000) is refuted by the current work. In fact, heat transfer is then enhanced thanks the efficiency of conduction along the roughness elements as in the study of Sahiti et al. (2005). Figure 8 presents the influence of roughness on $Po$ and the global $Nu$ number at constant flow rate and constitutes the synthesis of the present results. The Poiseuille number has been already introduced in section 2.2 and the $Nu$ number is equal to $Nu = 1/\theta_b$ where $\theta_b$ is the dimensionless bulk temperature defined as follows:

$$ \theta_b = \frac{1}{u_b H} \int \theta_f(y) u_D(y) dy \quad (20) $$

It should be reminded here that the PML model based on volume averaging uses and reports $\theta_f(y)$ which is the mean dimensionless temperature instead of the local bulk temperature so it does not account appropriately for the enthalpy rate convected by the fluid when the velocity varies in x and z directions. The presence of roughness introduces such variations of velocity so that $Nu$ based on
<figure><img src="image_9.png" /><figcaption>Figure 8. Increase of Po and Nu with the relative roughness height. Continuous lines: PML model; symbols: 3D simulations</figcaption></figure>

the definition given by equation (20) could be overestimated. This uncertainty introduced by the model in the Nusselt number increases with the roughness height and could be very important for high values of $k^*$. Figure 8 shows an increase of Po over the theoretical value. Po strongly increases for small values of Da as expected. On the other hand Nu increases with the PML thickness not as much as Po. As for Po the growth of Nu depends on the Darcy number. However, in the lower range of Da this dependency is not much pronounced. The results of the PML model agree very well with the results obtained from the 3D numerical simulations. This again confirms the consistency of the PML model with the 3D numerical

simulations. The results show that the ratio Nu/Po decreases when the PML thickness is increased and this trend is more pronounced for small Da.

The present results were compared with experiments performed on water flows through semi-rough rectangular microchannels. Two test-sections with the same geometrical arrangements of roughness as in the model were investigated (Figure 9). The rough wall was obtained by etching a silicon wafer at a depth $k$ using a mask reproducing the design of the roughness arrangement. The microchannel height was obtained by etching a pyrex plate at the depth $H-k$. The wafer and the pyrex cover were then anodically bounded to make a microchannel of height $H$ with a rough wall on one side and a smooth wall on the other side. The microchannels investigated were 106 µm or 153 µm in height, 3 mm in width and 32 mm in length. The relative roughness ($k/0.5H$) was 16% and 14% respectively. Both microchannels were tested in a closed-loop circuit. The pressure drop was measured by two pressure transducers placed upstream-downstream of the test-section. The flow rate was measured by high-accuracy flowmeters (Gao et al., 2002). The PML model was run with boundary conditions corresponding to a semi-rough channel. Empirical parameters of equations (19) were obtained by means of 2D simulations for the flow and heat transfer through the banks of rods in the staggered arrangement. Figure 10 shows an excellent agreement between the measured Poiseuille number and the predictions of the PML model. The measurements and the model are in consistency to show that

i) Po is independent of Re (no inertia effect) in the laminar regime ($Re < 2000$)

ii) The pressure drop is significantly increased (about 20%) for a relative roughness of about 15%

An extension of the experimental study to heat transfer measurements is planned.

<figure><img src="image_10.png" /><figcaption>Figure 9. Roughness elements etched in silicon wafer ($L_r=l_r=16\mu m$)</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Figure 10. Variation of Po with Re – comparison between experimental results and model predictions</figcaption></figure>
5. Conclusions

The 3D and PML models were adapted and developed in order to predict the influence of roughness on laminar flow in microchannels. In spite of its simplicity, the PML model is able to predict the pressure drop and heat transfer increase due to surface roughness. It was shown that the control parameters i.e. the local *Nu* number and the permeability only depend on the porosity and can then be efficiently modeled by means of simple 2D numerical simulations. This is confirmed by the agreement observed between results of the PML model and 3D numerical simulations. The PML model reports an increase of *Po* with the PML thickness while the increase of *Nu* is much weaker. This result is in agreement with the conventional fins theory.

References

Baviere R., Gamrat G., Favre-Marinet M., Le Person S., 2006, Modelling of laminar Flows in Rough-Wall Microchannels, J. Fluids Engineering (to appear in July 2006)

Baviere R., Ayela F., Le Person S., Favre-Marinet M., 2005, Experimental characterization of water flow through smooth rectangular microchannels, Phys. Fluids, 17, 098105, 1-4

Croce G., D'Agaro P., 2004, Numerical analysis of roughness effect on microtube heat transfer, Superlattices and Microstructures, 35, 601-616

Croce G., D'Agaro P., Nonino C., Zani F., 2005, Three-Dimensional Roughness Effect on Microchannel Heat Transfer and Pressure Drop, ECI Int. Conference on Heat Transfer and Fluid Flow in Microscale, Castelvecchio Pascoli

Gao P., Le Person S., Favre-Marinet M., 2002, Scale effects on hydrodynamics and heat transfer in two-dimensional mini and microchannels, Int. J. Thermal Sciences, 41, 1017-1027

Guo Z., Li Z., 2003, Size Effect on Microscale Single-Phase Flow and Heat Transfer, Int. J. Heat and Mass Transfer, 46, 149-159

Hu Y., Werner C., Li D., 2003, Influence of Three-Dimensional Roughness on Pressure-Driven Flow Through Microchannels, J. Fluids Engineering, 125, 871-879

Judy J., Maynes D., Webb B.W., 2002, Characterization of frictional pressure drop flows through microchannels, Int. J. Heat and Mass Transfer, 45, 3477-3489

Kim S.J., Kim D., 1999, Forced Convection in Microstructures for Electronic Equipment Cooling, J. Heat Transfer, 121, 639-645

Kohl M.J., Abdel-Khalik S.I., Jeter S.M., Sadowski D.L., 2004, An Experimental Investigation of Microchannel Flow with Internal Pressure Measurements, Int. J. Heat and Mass Transfer, 48, 1518-1533

Koo J., Kleinstreuer C., 2003, Liquid Flow in Microchannels: Experimental Observations and Computational Analyses of Microfluidics Effects, J. Micromechanics and Microengineering, 13, 568-579

Koo J., Kleinstreuer C., 2005, Analysis of Surface Roughness Effects on Heat Transfer in Micro-conduits, Int. J. Heat and Mass Transfer, 48, 2625-2634

Phares D.J., Smedley G.T., 2004, A study of laminar flow of polar liquids through circular microtubes, Physics of Fluids, 16, 1267-1272

Qu W., Mala Gh. M., Li D., 2000, Heat transfer for water flow in trapezoidal silicon microchannels, Int. J. Heat and Mass Transfer, 43, 3925-3936

Sahiti N., Durst F., Dewan A., 2005, Heat transfer enhancement by pin elements, Int. J. Heat and Mass Transfer, 48, 4738-4747

Taylor R.P., Coleman H.W., Hodge B.K., 1985, Prediction of Turbulent Rough-Wall Skin Friction Using a Discrete Element Approach, J. Fluids Engineering, 107, 251-257

Taylor R.P., Coleman H.W., Hodge B.K., 1989, Prediction of Heat Transfer in Turbulent Flow Over Rough Surfaces, J. Heat Transfer, 111, 568-572