MODELING OF ROUGHNESS EFFECT ON LAMINAR FLOW AND HEAT TRANSFER IN RECTANGULAR MICROCHANNELS

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published 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.

T he present paper is dev oted to the numerical modeling of roug hness effects on laminar flow in microchannels. 3 D numerical simulations were initially performed and, based on their solutions, a roug hness model is proposed. T his one-dimensional model is inspired by a discrete-element approach due to T aylor et al. ( 1 9 8 5 , 1 9 8 9 ). In this model, the channel consists in a clear medium adj acent to a porous medium layer where the interactions between the stream and the roug hness elements are directly computed by the discrete-element method. Empirical correlations based on 2 D numerical simulations of cross-flow throug h a bank of rods are used for modeling the local drag coefficient and the local Nusselt number of the roug hness elements. T he appropriate determination of these important parameters allows solv ing the momentum and energ y equations in order to obtain the v elocity and temperature profiles. T he 1 D model solutions are then compared with the 3 D numerical solutions showing a v ery g ood ag reement. T he present study shows that roug hness sig nificantly influences both the P oiseuille number and the g lobal Nusselt number. T he relativ e increase of the pressure drop is found to be much faster than that of the heat transfer coefficient when the roug hness heig ht is increased. T wo microchannels were produced with the same threedimensional roug hness arrang ement as in the 3 D numerical model and ex periments were performed with deioniz ed water as the working fluid. T he heig ht of microchannels was 1 0 6 µm and 1 5 3 µm with relativ e roug hness of 1 6 % and 1 4 % respectiv ely. Comparison between ex periments conducted in isothermal conditions and model solutions shows the g ood ability of the numerical models to predict the pressure drop in the roug h microchannels, which were tested.
No m e n c l a t u r e

. In t r o d u c t i o n
Recent adv ances in ev er faster and smaller electronic equipments hav e entailed the dev elopment of efficient cooling systems. T his has stimulated a strong current of research on flows in microchannels. Howev er many published papers present a strong dispersion of results and no definitiv e tendency is observ ed concerning the friction factor and the heat transfer coefficient in microflows. M easurements in microchannels are v ery difficult and require careful accounting for phenomena that are often neg lig ible in standard flows. Recent ex perimental works accurately ev aluated the uncertainty of measurements and their results are now consistent with the conv entional theory. F or ex ample, Kohl et al. ( 2 0 0 4 ) ex perimentally inv estig ated laminar water flow throug h smooth microchannels with hydraulic diameter rang ing from 2 5 to 1 0 0 µm. T hey used internal pressure measurements to av oid entrance effects. T he av erag e v alue of their results ex hibits an increase by only 4 % abov e the theoretical v alue, whereas the uncertainty is 1 1 % . T he authors did not observ e an early transition to turbulence. T hey concluded that the larg e inconsistencies in prev iously published data are probably due to instrumentation errors. T he same conclusion was drawn by B av iere et al. ( 2 0 0 5 ) from ex periments in smooth microchannels rang ing from 4 .5 to 2 1 µm.
Some researchers proposed that a possible ex planation of the friction factor increase abov e conv entional v alues could be an effect of Electrical Double Layer. Howev er, a recent work of P hares and Smedley ( 2 0 0 4 ) shows that electrov iscous effects are not observ ed in 1 5 0 µm microtubes. T hey ex amined also the effect of surface roug hness and they concluded that it plays an important role in laminar flow throug h microtubes. T hey reported an increase in P oiseuille number v alues by up to 1 7 % for a relativ e surface roug hness equal to 2 .5 % . No measurable effect of v iscosity was detected for three mix tures of different v iscosities.
J udy et al. ( 2 0 0 2 ) performed ex tensiv e measurements of pressure losses for laminar flow throug h different microtubes ( two different materials and two different shapes) in the rang e 1 5 -1 5 0 µm. A fter a careful ex amination of possible errors they reported that the uncertainties in the P oiseuille number are as hig h as 1 0 -1 5 % for smooth fused silica and 2 0 -2 1 .5 % for roug h steel microtubes. T hey did not observ e sig nificant dev iation from Stokes theory, howev er they concluded that the detection of additional phenomena arising in microflows is v ery difficult since their influence is below ex perimental errors at microscales.
Recent ex perimental research sug g ests that surface roug hness is mainly responsible for dev iation from conv entional theory in microchannels. Guo and Li ( 2 0 0 3 ) stated that for roug h wall microchannels, the form drag of roug h elements is predominant and could be reg arded as one reason leading to increase in the friction factor. T hus there is a strong interest to ex plore this effect precisely and to dev elop simple models which could accurately predict the influence of roug hness on fluid flow and heat transfer. F ollowing the conclusion of J udy et al. ( 2 0 0 2 ) about the relativ e importance of microscale effects and ex perimental uncertainties, numerical simulation seems to be an adequate tool for ex ploring the influence of a particular phenomenon on microflows. T his was the idea behind a recent study of Croce et al. ( 2 0 0 5 ), who inv estig ated numerically the effects of 3 D roug hness on microchannel heat transfer and pressure drop. T hey found a 1 6 % increase in the P oiseuille number and a much smaller increase of heat transfer coefficient for a relativ e roug hness heig ht equal to 2 .6 5 % . In their study, conically shaped roug hness elements were periodically distributed on the walls of a smooth channel. T he numerical study of Hu et al. ( 2 0 0 3 ) reports ex tensiv e results about the surface roug hness impact on pressure drop in microchannels for cubic shaped roug hness elements. T hey proposed an empirical model which accounts for additional pressure drop in terms of a relativ e channel heig ht reduction. In other words, the apparent pressure drop is accounted for by the conv entional v alue of the P oiseuille number used with a reduced channel heig ht. Howev er, the proposed formulas suffer from limitations and cannot be used with other shapes of roug hness elements. Koo and Kleinstreuer ( 2 0 0 3 ) proposed that the roug hness reg ion could be modeled by an equiv alent porous medium layer ( P M L hereafter). T hey modeled the additional v iscous forces due to roug h elements in terms of the Darcy number and they used a nonlinear term to account for inertia forces. B y means of these additional source terms introduced into the Nav ier -Stokes equation, they were able to reproduce the ex perimental results of Guo and Li ( 2 0 0 3 ). Howev er, in spite of g ood ag reement between P M L model and ex perimental results they did not relate the Darcy number to the g eometrical parameters of roug hness. In a nex t paper, Koo and Kleinstreuer ( 2 0 0 5 ) ex tended their P M L model to heat transfer phenomena in roug h wall microconduits. T hey proposed a oneequation model to account for heat transfer increase due to roug h elements. Howev er, such a oneequation model seems to be inadequate for cases where conv ection heat transfer is predominant so that the temperature distribution in roug h elements strong ly differs from that in the fluid. In such cases the use of two-equations models is adv isable as illustrated by Kim and Kim ( 1 9 9 9 ).
T aylor et al.
( 1 9 8 5 and 1 9 8 9 ) presented a discrete-element method for predicting skin friction and heat transfer in turbulent flows. T his approach is based on a v olume av erag ing technique and considers separately the form drag due to roug hness elements and the shear drag at the flat part of the surface. T he information about the roug hness element shape and siz e is directly introduced into the momentum equation by the drag coefficient when it is correctly estimated. T his approach was used in the work of B av iere et al. . T he choice of a reg ular pattern of roug hness elements allows using periodical boundary conditions and reduces the computational domain to one period as depicted by dashed lines in F ig ure 1 . T he other adv antag e of periodicity is the possibility to use v olume av erag ing in order to establish the g ov erning equations for the P M L model. Since the P M L model principles are the same for stag g ered and alig ned arrang ements, the presentation is restricted to the latter one. T he computational domain ex tends ov er one wav eleng th in the x and z directions and ov er the half channel heig ht ( H/ 2 ) in the y direction normal to the wall ( F ig ure 2 ). T he long itudinal and transv erse dimensions are thus equal to the roug h element spacing . T he microchannel wall consists of a block of thickness H s , heated at uniform heat flux ' ' 0 q on the ex ternal side. T he computational domain was treated as the central part of an ex tremely long channel so that the flow was considered as fully dev eloped in the long itudinal direction and the heat flux distortion due to ax ial conduction at the channel ends was assumed to be neg lig ible. A s a result, the flow properties are periodic in the x direction. T his assumption allows applying periodic boundary conditions on the opposite sides of the domain in the long itudinal direction. T he periodic boundary condition for any flux may be written: ( 1 ) Due to periodicity, state v ariables as pressure and temperature can be written as the sum of a linear g radient and a periodic component: T he term dx dp / relates only to the fluid domain and is a priori known in our computations while the term dx dT / is determined by the energ y conserv ation in the fluid and solid domains. It was then assumed that dx dT / in the solid and fluid phases is equal to where M & is the mass flow rate. Since dx dT / equally contributes to the temperature field in both domains, we can which is periodic: ( is the mean wall temperature at 0 = y . T he boundary condition as g iv en by equation ( 1 ) relates to the fluid and solid phases. F or the fluid domain, Θ r stands for the v elocity v ector while in the solid block it stands for the heat flux in ax ial direction ' ' a q r . T he ax ial heat flux in the solid block results from the constant temperature g radient dx dT / and is then ex pressed by where i r is the unit v ector in x direction. F or a larg e conductiv ity ratio f s k k / , the ax ial heat flux in the fluid was assumed to be neg lig ible. F or the sake of symmetry, the v elocity and heat flux g radient normal to the lateral sides and top plane of the computational domain were assumed to be z ero. T he symmetrical boundary conditions are presented in the form: ( 5 ) where n r is the v ector in the direction normal to the symmetry surface.
( 1 0 ) a n d d i s p o s i n g t w o a p p r o x i m a t i o n s T 1 a n d

Porous medium layer model
( 1 1 ) w her e s τ , p w a n d p l a r e r es p ec t i v el y t he s hea r s t r es s , t he w i n d w a r d a ( 1 5 ) w h e r e ) (  o n t i n u i t y o f a v e r a g e t e m p e r a t u r e a t r o u g h / c l e a r i n t e r f a c e . S u b s t i t u t i n g t h e l e f t -h a n d t e r m o f e q u a t i o n ( 1 5 ) i n t o e q u a t i o n ( 1 4 ) a n d e l i m i n a t i n g T f r e s u l t s i n :    . Comparis on of F ig u re 4 and F ig u re 5 s how s that h is more s ens itiv e than K to R e D . I t may be probably attribu ted to the v alu e of the P rand tl nu mber P r = 7 u s ed in the c u rrent s tu d y, w hic h means that the thermal bou nd ary layers d ev eloping on the rou g hnes s s u rf ac e are thinner than the c orres pond ing v eloc ity bou nd ary layers .

Boundary conditions
A f ter v olu me av erag ing , the problem of f lu id f low and heat trans f er in rou g h w all mic roc hannels red u c es to a s ys tem of s ec ond ord er ord inary d if f erential eq u ations . I t c an be s olv ed nu meric ally u nd er appropriate bou nd ary c ond itions . I t w as as s u med that the f lu id and s olid temperatu re are eq u al (θ s =θ D ) and the v eloc ity is eq u al to z ero at the bottom w all (y =0 ) . S ymmetry bou nd ary c ond itions w ere u s ed at the c hannel s ymmetry plane (y =0 . 5 H ) . T he c ontinu ity of temperatu re T D and D arc y v eloc ity u D is s atis f ied in the w hole c ompu tational d omain. T he bou nd ary c ond ition at the P M L / c lear reg ion interf ac e rais es a partic u lar problem. T he d ev elopment of v eloc ity and thermal bou nd ary layers on the rou g h elements top s u rf ac e c ontribu tes to the d is c ontinu ity of s hear s tres s and heat f lu x at the P M L / c lear reg ion interf ac e in the v olu me av erag ing approac h. I n ord er to ac c ou nt f or thes e d is c ontinu ities the pres ent mod el c ons id ers a c ontrol v olu me CV N d ef ined by y k y k δ + ≤ ≤ w here the momentu m and energ y balanc es are w ritten. S inc e K and Nu D are mos tly mod elled by means of 2 D nu meric al s imu lations , it w as nec es s ary to as s u me that the d ev eloping v eloc ity and thermal bou nd ary layers on the s id e and top s u rf ac es of rou g hnes s elements are s imilar. A s a c ons eq u enc e, the s hear s tres s at the rou g h elements top s u rf ac e is d ed u c ed f rom the v alu e of the av erag e f ric tion d rag c oef f ic ient C d f w hic h w as f ou nd at the s id e s u rf ac es of rou g h elements . I n the s ame w ay, the s ec ond bou nd ary c ond ition f or the energ y eq u ation in the s olid phas e as s u mes id entic al c onv ec tion heat trans f er c oef f ic ient h f or the top s u rf ac e and the s id e s u rf ac e of rou g h element. T he d imens ionles s s ys tem of g ov erning eq u ations is : F ig u re 7 s how s the prof iles of the f lu id θ f and s olid θ s temperatu res d ef ined by eq u ations (1 6 ) , f or the s ame P M L thic k nes s and D arc y nu mber as in F ig u re 6 . A s f or the v eloc ity prof iles , a v ery g ood ag reement is f ou nd betw een 3 D and P M L mod el s olu tions . T he apparent d is c ontinu ity of f lu id temperatu re θ f res u lts f rom the d is c ontinu ity of poros ity β and f rom the s trong heat f lu x at the top s u rf ac e of rou g h elements . T he obs erv ed d if f erenc es betw een f lu id and s olid temperatu res ind ic ate that in the c as e of s trong c onv ec tion the as s u mption of thermal eq u ilibriu m is not appropriate. T he rou g h elements of hig h thermal c ond u c tiv ity w hen c ompared to the c ond u c tiv ity of w ater (k f /k s = 5 . 1 * 1 0 -3 in the c u rrent s tu d y) , c ond u c t heat d irec tly tow ard s the reg ion of hig her v eloc ities . S inc e the hig hes t heat s ink oc c u rs in the c entral part of the c hannel, the c ond u c tiv e heat trans f er rate in s olid elements inc reas es w ith the rou g hnes s heig ht. T he heat rate is mos tly c ond u c ted throu g h the rou g h elements c au s ing a s mall temperatu re g rad ient both in the f lu id and s olid phas es bec au s e the heat f lu x throu g h the f lu id is s mall and the thermal c ond u c tiv ity is hig h in the s olid phas e. I t is then c lear w hy the temperatu re prof iles s how ed in F ig u re 7 bec ome f latter w hen the P M L thic k nes s inc reas es . T he obs erv ed s mall temperatu re g rad ient is in ag reement w ith the as s u mption reported by Q u et al.  H ow ev er, in the low er rang e of D a this d epend enc y is not mu c h pronou nc ed . T he res u lts of the P M L mod el ag ree v ery w ell w ith the res u lts obtained f rom the 3 D nu meric al s imu lations . T his ag ain c onf irms the c ons is tenc y of the P M L mod el w ith the 3 D nu meric al s imu lations . T he res u lts s how that the ratio N u /P o d ec reas es w hen the P M L thic k nes s is inc reas ed and this trend is more pronou nc ed f or s mall D a .
T he pres ent res u lts w ere c ompared w ith ex periments perf ormed on w ater f low s throu g h s emi-rou g h rec tang u lar mic roc hannels . T w o tes t-s ec tions w ith the s ame g eometric al arrang ements of rou g hnes s as in the mod el w ere inv es tig ated (F ig u re 9 ) . T he rou g h w all w as obtained by etc hing a s ilic on w af er at a d epth k u s ing a mas k reprod u c ing the d es ig n of the rou g hnes s arrang ement. T he mic roc hannel heig ht w as obtained by etc hing a pyrex plate at the d epth H-k.