Roland Bavière

Gabriel Gamrat

Michel Favre-Marinet

Stéphane Le

S Le Person

Modelling of Laminar Flows in Rough-Wall Microchannels

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.

Introduction

Much work in recent years has gone into the development of microchannel technology. This has been especially useful in electronics, where ever faster and smaller microchips are built, or in biology where lab on chips are currently used. Investigations of liquid flows in rough microducts present a great interest since rough surfaces are commonly found in many practical situations. The size of roughness elements is usually very small compared to macroscale in ducts of conventional size so that roughness effects may be neglected in laminar flows. However, the relative roughness becomes significant when the duct size is reduced to several hundred micrometers or less, and may play a role in the hydrodynamics, even in laminar flows. Surface roughness is due to the process used for making microchannels, and may be high when micromachining is used. In addition, rough walls are potentially interesting for enhancing heat transfer in industrial applications. It is therefore highly desirable to understand the roughness effect on flows in microducts and to quantify its influence on the pressure losses. This has motivated several studies using experimental, analytical and numerical approaches.

Mala and Li ͓1͔ experimentally found an increase in the Poiseuille number for water flows in rough microtubes ranging from 50 m to 254 m in diameter. For Re below 800, they measured an increase in pressure losses ranging from 7% to 15% for relative roughness ͑roughness element height divided by tube diameter͒ between 1.2% and 3.6%. For Re higher than 1000, the Poiseuille number was found to increase significantly with Re. This was interpreted by the authors as the result of an early transition to turbulence in rough microtubes, when compared to conventional size tubes. Guo and coworkers ͓2,3͔ also performed measurements in rough stainless steel microtubes. They concluded that a relative roughness of 3% to 4% increased the pressure drop from 15% to 37% when compared to the smooth case. Judy et al. ͓4͔ conducted experiments in rough stainless steel microtubes from 75 m to 125 m in diameter. There is no precise indication of the roughness height in their paper. They did not conclude that an increase in the pressure losses might be due to roughness elements. In their study, the Poiseuille number was found to be constant in the laminar regime. Measurements were also carried out in rough trapezoidal microchannels ͓5͔, obtained by wet-etching a silicon substrate and covering it with a Pyrex glass plate. The results revealed an increase in the Poiseuille number compared to the conventional theory in the Reynolds range investigated in the present work ͑Ͻ1500͒. For low Re, the increase in pressure losses varied from 8% to 30% for relative roughness ranging from 1.7% to 3%. These investigations supported the idea of an early transition to turbulence due to roughness elements. Pfund et al. ͓6͔ performed experiments in a rough rectangular microchannel 257 m in height. They found a 26% increase in the pressure losses for a mean relative roughness ͑roughness element height divided by channel height͒ of 1.5%. Recently, Baviere et al. ͓7,8͔ investigated high aspect ratio rectangular microchannels ranging from 100 m to 300 m in height. In their experiments, the relative roughness ranged from 1.7% to 5%. For Re between 50 and 1000, the corresponding increase in the Poiseuille number was found to be 5% to 37%.

To summarize, most published experimental results concerning flows in rough microducts show an increase in pressure losses compared to conventional smooth wall ducts. However, some aspects of these results are conflicting. In some experiments, the Poiseuille number increases significantly with Re. On the contrary, other results show Reynolds number independency of the Poiseuille number in a wide range of Re. It is therefore important to establish the relationship between the pressure drop increase and the rough surface characteristics.

An analytical approach of roughness effects was developed by Mala and Li ͓1͔ and Qu et al. ͓5͔. These authors proposed a roughness-viscosity model in a manner similar to the eddyviscosity concept in turbulent flows. The model brought into play a constant which needed to be evaluated from experimental re-sults. More recently, Koo and Kleinstreuer ͓9͔ modeled the nearwall region by using an equivalent porous medium layer. They were able to reproduce the measurements of ͓1,3͔ by adjusting the permeability ͑or equivalently the Darcy factor͒ of this layer. They also applied the model to microjournal bearings ͓10͔.

Numerical computations were performed by Hu et al. ͓11͔͑ denoted HWL hereafter͒, for water flows in two-dimensional microchannels with rectangular prism rough elements. Their calculations were done in channels ranging from 5 mt o5 0mi n height at Reynolds number ranging from 0.002 to 20. The authors expressed the roughness effect as a relative channel height reduction depending on the roughness elements geometry.

The aim of the present paper is to show that a very simple analytical model is able to predict the roughness effect on the flow in rough-wall microchannels. The analytical model is based on the method developed by Taylor et al. ͓12͔ for predicting the roughwall skin friction in turbulent flows. In this paper, we combined analytical and numerical approaches of the flow and compared the results to experimental data ͓8͔ and to published numerical results ͓11͔.

Numerical Approach

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

Geometrical Model of Roughness.

The geometrical model of roughness consists of blocks periodically distributed on the walls of a plane channel of height H ͑Fig. 1͒. The wavelengths are L r and l r , in the longitudinal x and transverse z directions, respectively. The rough elements are parallelepipeds of square cross section ͑side length d͒ and of height k placed on a smooth wall, called thereafter the bottom wall. They are distributed either in symmetrical or asymmetrical arrangements. Series of calculations were performed by varying the dimensions as follows:

Roughness height ͑k͓m͔͒: 0.1, 0.2, 0.3, 0.5, 0.7, 1, 2 Roughness side ͑d͓m͔͒: 0.5, 1, 1.5

Roughness spacing ͑L r = l r ͓m͔͒:2 ,3 ,4 Channel height ͑H͓m͔͒: 5, 7.5, 10

Computation Domain and Boundary

Conditions. The computational domain extends over one wavelength in the x and z directions and over the half-channel height ͑H /2͒ in the y direction normal to the wall ͑Fig. 1͒. The longitudinal and transverse dimensions are thus equal to the rough element spacing. The computational domain was treated as part of an extremely long channel so that the flow was considered as fully developed in the longitudinal direction. As a result, the flow properties are periodic in the x direction, except for the pressure which is composed of a linearly decreasing term and a periodic one. This assumption allows for applying periodic boundary conditions on opposite sides of the domain in the longitudinal direction. The periodic boundary condition may be written for the velocity field U͑x,y,z͒ = U͑x + L r ,y,z͒͑ 1͒

where U is the velocity vector.

The lateral sides of the computational domain are chosen as the planes of symmetry of the rough elements parallel to the xy plane. As a result, the velocity gradient of U and V normal to these side planes and the transverse velocity component W are assumed to be zero. at z = 0 and z = l r dU/dz =0, dV/dz =0, W =0 ͑2͒

The condition of symmetry at the top plane of the computational domain is written at y = H/2 dU/dy =0, V =0, dW/dy =0 ͑3͒

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

Equations.

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

Continuity equation

١U =0 ͑4͒

Momentum equation

U ١ U =-١P + ١ 2 U ͑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 SIM-PLEC ͑semi-implicit pressure linked equations consistent͒ algorithm was used. This algorithm is based on a prediction-correction method, which allows the equations to be linearized and solved iteratively. The pressure under-relaxation factor was set to 0.5. The calculations were performed by means of a double precision solver until the level of residuals decreased below 10 -12 . Computations were performed on orthogonal grids generated by GAMBIT 2.1.2.

Numerical Accuracy.

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

T 3 = T 1 + T 1 -T 2 ͑h 2 /h 1 ͒ 2 -1 ͑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.

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 ␦y in the direction normal to the wall, extending over a large distance compared to the rough elements spacing in the directions ͑respectively, ␦x, ␦z ӷ L r , l r , Fig. 2͒. Flow quantities are averaged over the lengthscales ␦x, ␦z as in homogenization methods.

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͒ = 1 L r l r -d 2 ͵͵ dS͑x͒ u͑,y,͒dd for 0 ഛ y ഛ k ͑7͒
where u͑ , y , ͒ stands for the actual flow velocity in the roughness region ͑y ഛ k͒. Integration is performed over the surface dS͑x͒ open for flow in an elementary cell of sides l r ϫ 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 ␦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 ␦x. For the particular case of periodically distributed rough elements, the averaged blockage coefficient ␤ is given by

␤ = 1 L r l r ͵͵ dS͑x͒ dd ͑8͒
and is therefore

␤ =1- d 2 L r l r for 0 ഛ y ഛ k ͑9͒ ␤ = 1 for k ഛ y
This is the same expression as in ͓12͔, except that the cross section of the rough elements is square instead of being circular in ͓12͔. In the present model of cylindrical elements, ␤ is independent of y. The maximum blockage effect corresponds to ␤ =0, when the rough elements are adjacent.

We also use a filtration or Darcy velocity u D ͑y͒, in the terminology of flow in porous media. This velocity is defined with the total surface L r l r u D ͑y͒ = 1

L r l r ͵͵ dS͑x͒ u͑,y,͒dd ͑10͒ so that U͑y͒ = u D ͑y͒/␤ ͑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.

p ¯͑x,y͒ = 1 L r l r -d 2 ͵͵ dS͑x͒ p͑,y,͒dd for 0 ഛ y ഛ k ͑12͒
We suppose that p ¯͑x , y͒ is independent of y in the present mean one-dimensional flow. From here, we call P͑x͒ the averaged pressure defined by Eq. ͑12͒.

Momentum Equation. The equilibrium of a control volume results from the competition between

Pressure forces on the upstream and downstream sides of CV Viscous shear stresses on the upper and lower sides of CV Drag forces due to the rough elements The drag force due to the portion of a rough element included in CV is modeled by using a drag coefficient C d

␦F d = 1 2 U 2 ͑y͒C d d␦y ͑13͒
Finally, the momentum equation is

0=- d dx ͑␤P͒ + d dy ͩ␤ dU dy ͪ-1 2 C d d U 2 L r l r for 0 ഛ y ഛ k ͑14͒ 0=- dP dx + d 2 U dy 2 for k ഛ y ഛ 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.

Boundary Conditions.

A difficulty arises for the rough/ clear layer interfacial conditions for the present shape of the roughness elements. In fact, the blockage factor is discontinuous at the interface, going from 1 to a constant value ␤ Ͻ 1 across the interface. This difficulty would not exist for other shapes of the roughness elements, like conical-shaped blocks, as in the Taylor et al. paper ͓12͔. The continuity of the local velocity across the rough/clear layer interface leads to the continuity of the Darcy velocity. As a result, the spatially averaged velocity U͑y͒ is discontinuous like the blockage coefficient at the interface.

␤U͑y = k -͒ = U͑y = k + ͒͑ 15͒
As mentioned before, this discontinuity would not exist for roughness elements of gradually decreasing cross section. In fact, it is likely that the blockage effect extends over a small distance within the clear medium away from the interface because the fluid is strongly slowed down in the vicinity of the top of the blocks. The present model does not account for this effect and this is probably a weakness of the model. It could probably be improved by introducing a smoothing of the blockage effect near the interface. This was not attempted in the present work to preserve the simplicity of the model.

In order to express the boundary condition on the shear stress at the interface, we consider a control volume of height ␦y and lengths ␦x, ␦z in the other directions, centered at y = k. The force balance on this CV requires

ͫ ͯ dU dy ͯ y=k+␦y/2 -͑1-␤͒ͯ dU dy ͯ y=k + -␤ͯ dU dy ͯ y=k-␦y/2 ͬ ␦x␦z + ͫ dP ͑1+␤͒ 2 ␦z - 1 2 U 2 ͑k͒C d d␦x␦z 2L r l r ͬ ␦y =0 ͑16͒
When ␦y tends to zero, the pressure force and the drag force become negligibly small so that Eq. ͑16͒ reduces to

ͯ dU dy ͯ y=k - = ͯ dU dy ͯ y=k +

͑17͒

or using the Darcy velocity

1 ␤ ͯ du D dy ͯ y=k - = ͯ du D dy ͯ y=k +

͑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 ͑19͒ ͯ ץU ץy ͯ y=H/2 =0 ͑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 Dimensional analysis shows that the drag force may be written

U͑y = H͒ =0 ͑21͒
␦F d ͑y͒ = 1 2 Kͩ k d , y k ,Re d ͪͩ U͑y͒ d ͪd␦y ͑22͒
␦F d is also related to the drag coefficient C d by Eq. ͑13͒ so that

K = C d Re d .
Eliminating C d by using Eqs. ͑13͒ and ͑22͒, the last term of Eq. ͑14͒ becomes 1 / 2KU / L r l r .

It is well known ͓14͔ that inertia effects cannot be totally discarded for creeping flows near extremely long cylinders. However, the numerical simulations of the present work and those of HWL along with some experimental works ͓4,7͔ indicate that roughness effects are Re independent up to Re D h = 200. Consequently we assume that K is a function of y / k and k / d only. The presence of these parameters in the coefficient K reflects the threedimensional nature of the low-Reynolds number flow near the blocks. For creeping flows normal to the major axis of very slender prolate spheroids ͑E ӷ 1͒, the drag force is given by ͓15͔

D = 4U2a 1 2 + Ln 2a b ͑23͒
where E = a / b,2 a = major axis, 2b = minor axis of the spheroid. Following the form of Eq. ͑23͒, we propose for K the following expression:

K = K 1 1+K 2 Ln k d for k d ജ 1 ͑24͒ K = K 1 for k d ഛ 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 1 was therefore deduced from the numerical simulations conducted in typical arrangements of the rough elements. For a given flow field, the elementary drag force ␦F d was obtained by integrating the pressure forces on the front and rear sides of a block slice of height ␦y and the viscous forces on its lateral sides. The numerical results have shown that the coefficient K 1 is nearly independent of y in the upper part of the blocks, as presented later. For k / d ഛ 1, K 1 was found to be 116 for the following geometrical arrangement of the rough elements ͑ d / H = 0.2, 0.1Ͻ k / H Ͻ 0.2, L r / d = l r / d =2͒. For k / d Ͼ 1, the best agreement of the model with the HWL results was obtained for K 2 = 0.5. The dependence of K on the geometrical parameters was not systematically explored in the present work. The same constants were used for all the geometrical arrangements considered in this study. The model could probably be improved with a better estimation of K as a function of the actual arrangement of the rough elements.

3.1.5 Analytical Solution. The system of equations is reduced to dimensionless form by normalizing the lengths with the channel height H and the effective velocity with the bulk velocity of the Poiseuille flow in a smooth channel

U * = U - 1 dP dx H 2 12 ͑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 ഛ y * ഛ k * U * ͑y * ͒ = 12 A * Ά1+ ͱ A * ͩ 1 2 -k * ͪsi nh͑ ͱ A * y * ͒ -cosh͓ ͱ A * ͑k * -y * ͔͒ C 1 · ͑26͒ for k * ഛ y * ഛ 1 2 U * ͑y * ͒ = 12␤ A * ΄1+ ͱ A * ͩ 1 2 -k * ͪS 1 -1 C 1 ΅+6͓ y * -y * 2 + k * 2 -k * ͔
where

A * =1/2K / L r * l r * ␤, C 1 = cosh͑ ͱ A * k * ͒, S 1 = sinh͑ ͱ A * k * ͒
A * ͑1-e -ͱ A * y * ͒ +2D 1 sinh͑ ͱ A * y * ͒ ͑27͒ For k * ഛ y * ഛ 1 U * ͑y * ͒ = ͑y * -1͒ ͫ 2D 1 ͱ A * C 1 +12 * E 1 ͱ A * ͪͬ +6͑1-y * 2 ͒ where E 1 = exp͑-ͱ A * k * ͒ D 1 = 6 ͱ A * ͑k * -1͒C 1 -␤S 1 ͫ k *2 2 - 1 2 + ͑1-k * ͒ ͩ k * + E 1 ͱ A * ͪ + ␤ͩ 1-E 1 A * ͪ ͬ
Integration of the velocity profile gives the flow rate Q as a function of the pressure gradient, the geometrical parameters k * , ␤, and of A * .

Results

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

Drag Coefficient of the Roughness Elements.

As can be seen in Fig. 1, the blocks were divided into separate slices of height ␦y for calculating the local force density exerted by the fluid as a function of the local effective velocity U͑y͒. This approach allowed for computing the coefficient K͑y͒ defined by Eq. ͑22͒ from the numerical results. Figure 4 displays the variation of K with the position along a rough element ͑y = 0 corresponds to the bottom wall͒ for a typical value of k * ͑=0.2͒. It shows that K slightly varies in the upper part of the block ͑y / k Ͼ 0.3͒. The strong increase in K observed for the lowest part of the block is obviously due to the influence of the bottom wall and associated three-dimensional effects. The same computation used with a condition of slip on the bottom wall gave a constant value of K all along the block. However, this deviation of K from a constant value is of little consequence on the solution given by the model because the contribution of the drag force to the flow dynamics is of minor importance in the near-wall region as shown in the next section. It follows that the assumption of constant K used in the analytical approach is quite reasonable.

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

Influence of the Roughness Elements Geometry on the Poiseuille Number.

For given channel height and roughness geometry, the flow is characterized by the Poiseuille number Po = -1

dP dx 2H 3 Q ͑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 = -1

dP dx 2H r 3 Q ͑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.

H r H

= ͩ 24 Po ͪ 1/3

͑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.

Influence of the Roughness Element

Height. The apparent channel height reduction is moderate up to k * Ϸ 0.05 ͑Fig. 7͒, then increases rapidly with k * . The graph has also been drawn for very high values of k * ͑0.4͒ which are obviously very far from usual physical situations. There is an excellent agreement between the three approaches. The slight differences of the present numerical results with those of HWL may be explained by the fact that HWL interpolated their results and that inaccuracies may result from this approximation. The expected maximum value for 1 -H r / H is obviously 2k * and is displayed in Fig. 7. The figure reveals that the actual blockage effect is about 25% smaller for the present arrangement of rough elements when k * is high ͑Ͼ0.1͒. The blockage effect is less pronounced for small rough elements. The analytical model is consistent with this statement since the drag force ͑or A * in Eq. ͑26͒͒ is independent of d when the drag coefficient K is taken as a constant. However, the model used with K 2 = 0 significantly overestimates the channel height reduction as given by HWL for d * Ͻ 0.1. The agreement is much better when the drag coefficient is used with K 2 = 0.5 ͑Eq. ͑24͒͒. When d * tends to zero, the aspect ratio of the blocks becomes very high and the correction factor K 2 Ln͑k / d͒ in Eq. ͑24͒ plays a more important role in the drag. This suggests that three-dimensional effects must be taken into account in F D in this case.

Influence of the Roughness Element

Spacing. The other way to vary the blockage coefficient is to change the rough elements spacing while keeping their shape constant ͑k * = d * = 0.2 in Fig. 9͒. Unlike the preceding case, the channel height reduction continuously decreases from the maximum value 2k * to zero when the rough element spacing is increased. This result suggests that the role of the drag force becomes less important when the rough elements spacing is increased. This is consistent with the analytical model where the term corresponding to the drag force in Eq. ͑14͒ is inversely proportional to L r l r . The roughness effect becomes vanishingly small when the density of the roughness elements tends to zero while keeping their shape constant. Again, Fig. 9 reveals an excellent agreement of the analytical model with the numerical simulations. 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 10. The trend exhibited by Po exp is well reproduced by the analytical and the numerical models. Considering the analytical model, the best agreement with the experimental results is obtained for k = 7.2 m, which is slightly higher than the estimated height of the roughness elements. The same tendency may be deduced from the numerical model: A better agreement with the experimental data would have been obtained for higher roughness elements. Apart from the experimental and numerical errors and the uncertainties inherent in the analytical model, this discrepancy may also be due to the experimental surface finish. In fact, it is possible that the layer of nickel was not perfectly smooth but may have included small rough elements besides the largest particles of SiC. As a consequence, the assumption of a smooth bottom wall both in the analytical and numerical models would underestimate the friction on this part of the wall.

Conclusions

Numerical modeling and analytical approach were used to compute laminar flows in rough-wall microchannels. Both models considered the same arrangements of roughness elements in periodical arrays. The numerical results have confirmed that the flow is independent of the Reynolds number in the range 1-200.

The analytical approach was based on a previous work ͓12͔ devoted to the prediction of the rough-wall skin friction in turbulent flows. The analysis has pointed out the difficulties of setting the boundary conditions at the clear/rough layer interface in the case of rectangular prism rough elements. For given geometrical characteristics of the roughness elements, the analytical model only needs two constants K 1 and K 2 to compute the flow. They correspond to the drag coefficient of the rough elements. In most cases, a good approximation of the flow characteristics is obtained by using only one constant, namely, K 1 . Despite the simplicity of the model, a very good agreement is found with the numerical results of the present work and those of Hu et al. ͓11͔. The results show that the rough elements drag is mainly responsible for the pressure drop across the channel in the upper part of the relative roughness range ͑k * Ͼ 0.03-0.06͒. For very high rough elements ͑k * Ͼ 0.1͒, the wall friction has a negligible contribution to the pressure drop. For a concentrated arrangement of the rough elements ͑d * = L r * /2=l r * / 2 = 0.2, Fig. 7͒, a very simple approximation of the roughness effect consists in considering a relative channel height reduction equal to 2k * . The accuracy is about 25% in 1 -H r / H for k * Ͼ 0.1. This blockage effect reduces when the rough elements are small or more dispersed ͑increasing blockage coefficient ␤͒. However, it remains significant for slender blocks, even if the side length becomes very small.

This study was restricted to patterns of periodic rough elements. The good agreement with experimental results obtained with randomly distributed roughness on microchannel walls is therefore striking. It seems that the models catch the main features of the physical phenomena involved in these microflows. The two main ingredients of the analytical model are the blockage coefficient ␤ and the relative roughness height k * . The knowledge of these parameters seems therefore to be of primary importance to predict the roughness effects in rough microchannels. This study will be soon extended to roughness effects on heat transfer in microchannels.

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. ϭ dynamic viscosity, kg m -1 s -1 ϭ kinematic viscosity, m 2 s -1

Nomenclature

Fig. 1

 1 Fig. 1 Computational domain and arrangement of the rough elements

Fig. 2 Fig. 3

 23 Fig. 2 Sketch of a control volume over the rough surface. Side view.

  3.1.4 Drag Coefficient. The model relates the drag force ␦F d ͑y͒ exerted by the fluid on a slice of rough element of height ␦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 ϱ = U͑y͒. This flow is mainly characterized by the local Reynolds number Re d = U͑y͒d / , which is based on a typical dimension of a block cross section, namely, the side length d.R e d is small in most microchannel flows: Re d Ϸ 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.

  For a semirough channel, the expressions are For 0 ഛ y * ഛ k * U * ͑y * ͒ = 12

Fig. 4

 4 Fig. 4 Distribution of the drag law coefficient K along a rough element. k * = d * = L r * /2=l r * /2=0.2, ␤ = 0.75, Re= 200.

Fig. 6

 6 Fig. 6 Contribution of viscous and drag forces to the pressure gradient. d * = L r * /2=l r * /2=0.2, ␤ = 0.75. Solid diamonds and dotted lines indicate friction on the bottom wall: F w / F T , solid circles and dash-dotted lines indicate friction at the top of the rough elements: F t / 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.

Fig. 7

 7 Fig. 7 Influence of the roughness element height on the channel height reduction. d * = l r * /2=l r * /2=0.2, ␤ = 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 * .

4. 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͔.I n these experiments, the channel walls consisted of parallel bronze blocks separated by a foil of thickness e f ͑0.1 mmϽ e f Ͻ 1m m ͒ with a hollowed out central part 25 mm in width. The length of the microchannel was 82 mm. The blocks were locally treated by electrochemical deposition of a thin Ni layer ͑thickness e Ni Ϸ 2 m͒, together with small SiC particles ͑5t o7m in height͒.

  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.

Fig. 8

 8 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 ␤. 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.

Fig. 9

 9 Fig. 9 Influence of the roughness element spacing on the channel height reduction. k * = d * = 0.2. Same symbols as in Fig. 8.

  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 Po ϭ 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 ϭ fluid density, kg m -3 w ϭ wall shear stress, N m -2


G Grégoire

M Favre-Marinet

Saint Amand F Julien

Modeling of Turbulent Fluid Flow Over a Rough Wall With or Without Suction

The turbulent flow close to a wall with two-dimensional roughness is computed with a two-layer zonal model. For an impermeable wall, the classical logarithmic law compares well with the numerical results if the location of the fictitious wall modeling the surface is considered at the top of the rough boundary. The model developed by Wilcox for smooth walls is modified to account for the surface roughness and gives satisfactory results, especially for the friction coefficient, for the case of boundary layer suction.

Introduction

Although numerical modeling of turbulent flows is continuously progressing, some questions remain open and need further investigation in many practical situations. Among them, the problem of modeling the flow over rough surface corresponds to a very common situation in industrial or geophysical applications and it has therefore received much attention from many researchers. A review of the different works on this subject can be found in Patel ͓1͔. The simplest way of modeling high Reynolds number wall flows, and historically the oldest one, is to use a wallfunction approach and to match the mean velocity profile to the logarithmic law near the wall. When a rough wall is considered, the velocity profile is well described by introducing a shift, called the roughness function and denoted ⌬B, in the logarithmic law of the wall, ͓2,3͔. Such a wall-function approach is desirable in many situations since it avoids detailed computations of the flow in the viscous sublayer and in the buffer layer and therefore saves large computation time. However, Patel ͓1͔ underlines ''the uncertainty in the dependence of ⌬B on the size and type of roughness and also in the effective location of the fictitious wall, from which the distance is measured.'' Patel showed some inadequacy of the logarithmic law for modeling a turbulent flow over a wavy wall, which can be considered as a kind of rough surface and application of the wall-function approach is questionable for this type of surfaces.

Turbulent flows over porous walls with blowing or suction have been extensively investigated. A simple model was proposed by Stevenson ͓4͔ and in a modified version by Simpson et al. ͓5͔. Wilcox ͓6͔ suggested modifying the von Karman's constant and obtained results consistent with the experimental investigations of Andersen et al. ͓7͔. However, these models did not specify the surface roughness and to our knowledge the problem of modeling a turbulent flow over a very rough wall with suction has not been considered yet.

The present study is concerned with turbulent flows over a wall with two-dimensional periodic roughness. This type of roughness is used in the paper and pulp industry and in industrial heat exchangers, for example. The purpose of this work is twofold, first to clarify the issue of the effective location of the wall for this type of very rough impermeable surface and secondly to propose a modification of the suction or blowing laws when they are applied to such rough walls.

The flows were computed by using the Fluent™ CFD software with a two-layer zonal turbulence model, which is suitable for modeling the flow in the near-wall region. The law of the wall commonly used in the wall-function approach was compared to the results of the numerical computation.

Numerical Approach

2.1 Flow Conditions. The present study is motivated by optimization of the flow in pressure screens, which are used in the recycling paper industry to filter the paper pulp, ͓8͔. In these devices, the pulp is forced to pass through very fine slots ͑width 0.1 mm͒, which are machined in the screen wall in order to retain the contaminants whereas the useful fibers are entrained by the flow across the slots. In order to avoid floculation of fibers, a tangential motion parallel to the wall is imposed to the flow ͑for details, see ͓8͔͒. This tangential velocity component is due to the entrainment of the fluid by two-dimensional foils moving in the pulp close to the wall ͑typical distance: 12 mm͒. In a reference frame fixed to the screen wall, the flow near the surface may be considered as a steady shear layer of uniform thickness perturbed periodically by the passage of the foils. Unsteady phenomena are localized during the passage of the foils and may be disregarded during the rest of the time in a first approach. As a result, the wall region is modeled in the present study as a layer of uniform thickness submitted to steady flow conditions.

The shape of the rough wall in the present work corresponds to that of industrial pressure screens. It consists of a periodic jagged profile ͑Fig. 1͒. The geometry considered and the flow are twodimensional. The flow rate across the slots and the flow rate of the tangential motion are adjusted independently in the computation as well as in the industrial situation. The particular case of null flow rate across the slots corresponds to impermeable wall.

Computation Domain and Boundary Conditions.

The wall profile is characterized by periodic steps of height k s , equal to 1.2 mm and of wavelength L, equal to 3.2 mm ͑Fig. 1͒. Only three waves of the pattern are considered in the computation domain. The domain height H is kept constant for all the computa-tions, according to the assumptions discussed above. A more general case would correspond to the problem of developing boundary layers, but is not considered here.

The computed flow is assumed to be fully developed in the x-direction. It is therefore supposed to repeat periodically with the wavelength of the wall profile. Since the length of the computation domain is a multiple of the wavelength, periodic boundary conditions are assumed at the inlet and outlet sections ͑i͒ and ͑o͒ of the domain. In the case of impermeable wall, the flow rate Q per unit span across these sections is prescribed. In addition, a condition of symmetry is assumed on the upper boundary ͑u͒ of the domain, which is located at a distance H from the bottom of the troughs of the rough surface ͑Fig. 1͒.

When suction is applied through the wall, the velocity distribution in the slots is assumed to be uniform ͑suction velocity V p ). A uniform downward velocity V por through the upper boundary is introduced in order to balance the outflow through the slots:

V por ϭV p d/L.
(1)

It is worth noting that the suction velocity V por is constant along the upper boundary in keeping with the assumption that u is independent of x along the computation domain. Strictly speaking, the velocity field is periodic in x, but the influence of the jagged profile is limited to a narrow band near the surface ͑Section 3.1͒.

In the case of suction, it was no longer possible to set simultaneously a condition of symmetry and the above condition ͑Eq. ͑1͒͒ on ͑u͒. A constant tangential velocity V t was then assigned on ͑u͒͑ details are given in Section 3.3͒.

Turbulence Model.

The numerical computation employed the two-layer zonal model ͑denoted TLZM thereafter͒ first developed by Wolfstein ͓9͔ and later by Chen and Patel ͓10͔. The present version is a combination of the k-RNG ͑renormalization group͒ model used far from the wall and a mixing length model used near the wall. The boundary between the two layers is defined by a turbulent Reynolds number, Re y 1 , based on the distance y 1 of a point ͑cell center in the numerical computation͒ to the nearest wall

Re y 1 ϭk 1/2 y 1 / (2)
where k is the turbulent kinetic energy.

For Re y 1 Ͻ200, the one-equation model of Wolfstein is used. The turbulent viscosity is computed from

t ϭC k 1/2 l l ϭc 1 y 1 ͑ 1Ϫexp͑ϪRe y 1 /A ͒͒.
(3)

The dissipation rate is modeled by introducing a second length scale l ϭk 3/2 /l .

l ϭc 1 y 1 ͑ 1Ϫexp͑ϪRe y 1 /A ͒͒. (4)
The constants proposed by Chen and Patel ͓10͔ are

c 1 ϭC Ϫ3/4 , A ϭ70, A ϭ2c 1 (5)
with ϭ0.41, C ϭ0.09

2.4 Numerical Scheme. We employed a hybrid mesh in order to optimize accuracy while insuring a reasonable time of computation. Rectangular cells were used near the wall in order to control most efficiently the distance to the wall and to have sufficient accuracy in the viscous sublayer. At least ten cells were placed in the near wall region (Re y 1 Ͻ200). Triangular cells were used in the external flow (Re y 1 Ͼ200) in order to minimize the cells number while keeping a low cells skewness when the size of the cells was increased. The whole mesh consisted of about 15,000 cells.

The equations were discretized by means of a second-order accurate finite volume method. As these equations are nonlinear, a SIMPLE ͑semi-implicit pressure linked equations͒ algorithm was used. This algorithm is based on a prediction-correction method, which allows the equations to be linearized. The drawback of this iterative method is the convergence slowness. It was checked in these calculations that the number of iterations is proportional to the number of cells. The computation started with uniform values of the various physical quantities. The software Fluent™ V6.0 was used on a PC HP Vectra Xu.

Numerical Accuracy.

Tests were conducted for a typical case of flow with a moderate suction rate (V p ϭ3.5 m/s, V t ϭ15 m/s). The kinetic turbulent energy k and the rate of dissipation of the turbulent kinetic energy are the physical quantities, which converge with the slowest rate. The level of normalized residuals for k and reached 10 Ϫ4 for both quantities after 2500 iterations, 5 10 Ϫ6 and 10 Ϫ5 , respectively, after 10,000 iterations. Additionally, the flow rate near the wall was computed to follow the convergence of the calculations. It was defined in a cross section over a wall crest (xϭL)b y

q w ϭ ͵ yϭ0 yϭ0.5 mm udy (6)
The variation of q w during 25,000 iterations showed that a minimum ͑Ϫ6% relative to the final value͒ was reached after 2500 iterations. The asymptotic value was approached with a good accuracy ͑ϩ2%͒ after 10,000 iterations. A first computation ͑case 1͒ was conducted with a grid of 15,700 cells, then the grid was refined in the region of high dissipation rate ͑downstream of the separation point at the wall crest͒. This case 2 used 18179 cells. Further refinement ͑case 3, 20,555 cells͒ was considered near the stagnation point, near the vortex center ͑Section 3.1͒ and in the region of high values of k. Finally, the grid of case 2 was otherwise refined ͑case 4, 19,118 cells͒ near the separation point and in the shear layer, which develops in the downstream direction behind this point. The influence of the resolution of the grid on the mean velocity profile over the wall crest was then investigated. Table 1 compares the computed values of the x-velocity with the results of case 4, which was considered as the more accurate solution, for several distances to the wall. According to Table 1, the mean velocity profile is computed with an accuracy better than 1% over most of the test section for cases 2 and 3. Errors can be estimated to 2-3% for case 1. The differences with the most refined case grow to 4 -7% very near to the wall. The skin-friction coefficient is also slightly affected by the grid refinement ͑differences less than 1% between cases 2, 3 and 4, 7% between cases 1 and 4͒.

y/k s V/V t ͑Case 4͒⌬ V/V t ͑Case 1͒⌬ V/V t ͑Case 2͒⌬ V/V t ͑Case

Results

Description of the Flow.

A typical result showing the mean flow field is drawn in Fig. 2 for a moderate suction rate (V p /V t ϭ1/3). The back-facing steps produce flow separation and create recirculation bubbles between two successive crests. Two particular streamlines may be considered for helping description of the flow. A streamline L A starting slightly upon a wall crest ends at a stagnation point A on the opposite side of the wall profile. On the other hand, a streamline L B starting at a wall crest C ends at a stagnation point B on the backside of the wall profile. L B is then the boundary of the recirculation bubble. L A and L B determine a streamtube corresponding to the flow through the slot. When the wall is impermeable, the flow features are essentially the same, except that only one particular streamline L A is to be considered. In this case, L A starts from C and it is the boundary of the recirculation bubble.

It is worth noting that the mean streamlines are only slightly perturbed by the wall even for small distances to the crests ͑Fig. 2͒.

Check of the Logarithmic Law for Impermeable

Rough Wall. The purpose of the present section is to compare the mean velocity profile resulting from the preceding TLZM computations with the logarithmic law which is used in the wallfunction approach for the case of impermeable wall (V p ϭ0). In this approach, the mean velocity profile is described in the fully rough regime, ͓3͔,b y

u͑Y ͒/u*ϭ1/ ln͑Y /k s ͒ϩB (7)
where u* denotes the friction velocity equal to ͱ w / ( w is the wall shear stress͒. Bϭ8.5 for very rough walls, ͓3͔.

Y is the coordinate normal to the wall. In fluid mechanics textbooks, the origin of Y is not always clearly defined. In fact, the formulation of Eq. ͑7͒ implies that the wall is represented by a fictitious plane surface, which corresponds to the origin of Y. For very rough walls, it is crucial to clearly specify the effective location of this fictitious plane surface when using Eq. ͑7͒. The complementary issue of the present discussion concerns therefore the position of the fictitious wall which gives the best agreement of the logarithmic law with the mean velocity profile as given by the previous numerical computations. For sake of clarity, the results will be compared in a reference frame which origin is located at the top of the rough surface ͑point C of Fig. 1͒.yi st h e ordinate in this reference frame. If y w denotes the location of the virtual surface introduced in the logarithmic law ͑Fig. 1͒, Y and y are related by

Y ϭyϪy w . ( 8
)
The mean velocity profile is then given by u͑ y ͒/u*ϭ1/ ln͑͑yϪy w ͒/k s ͒ϩB.

In order to compare the solutions associated to different positions of the virtual wall, u* was adjusted for each value of y w in Eq. ͑9͒ so as to give the prescribed mass flow rate Q through the cross section CM ͑Fig. 1͒ over a crest. The relation between Q and u* was obtained by integrating the logarithmic law ͑Eq. ͑9͒͒ from yϭ0t oHϪk s . The following relations are obtained:

u*ϭ Q D͑HϪk s ,y w ͒ϪD͑ 0,y w ͒ϩB͑ HϪk s ͒ (10)
with D͑ y,y w ͒ϭ͑ yϪy w ͒͑ln͑͑yϪy w ͒/k s ͒Ϫ1 ͒/ (11)

V t /u*ϭ1/ ln͓͑H/k s Ϫ1Ϫy w /k s ͔͒ϩB. (12)
In addition, the skin-friction coefficient was obtained from the above relations by

C f ϭ2͑u*/V t ͒ 2 . ( 13
)
The equivalent coefficient was computed with TLZM by integrating the total stress exerted on the wall profile contour C in the computational domain ͑1: curvilinear abscissa, element of length dl͒. 14) where p is the pressure, n x and x are, respectively, the x-components of the normal unit vector and of the computed wall shear stress. C f includes the pressure and friction effects for this nonhorizontal surface. This definition of an equivalent friction coefficient is the same as in Taylor et al. ͓11͔. It is here nondimensionalized with the length 3L of the equivalent plane virtual wall. The above calculations were performed for k s ϭ1.2 mm, H ϭ12 mm, Qϭ0.15 m 3 /s/m which were of industrial interest. The Reynolds number based on external velocity and half-channel height is V t H/ϭ2.10 5 . It should be noted that the present geometry corresponds to a very rough wall (k s ϩ Ϸ1000). Table 2 and Fig. 3 show the results obtained for two extreme positions of the virtual wall, at the crests (y w ϭ0) and at the bottom of the troughs (y w ϭϪk s ) of the actual wall profile. The mean velocity profile computed over a crest by TLZM is also plotted in Fig. 3. For this type of roughness, Fig. 3 clearly shows that the best agreement with the numerical results is obtained when the virtual wall is placed at the crests of the actual wall.

C f ϭ 1 V t 2 /2 1 3L ͩ ͵ C pn x dlϩ ͵ C x dl ͪ (
This conclusion is confirmed when computed velocity profiles are considered at three evenly distributed sections L/3 apart ͑Fig. 4͒. The velocity profiles almost collapse when the normalized distance to the wall, y/k s is higher than 0.4. Again the logarithmic law used with y w ϭ0 ͑Eq. ͑9͒͒ gives a good approximation of the velocity profiles. The actual mean velocity profile departs, however, significantly from the semi-logarithmic law for y/k s Ͻ0.3. Velocities are higher in this region because a mixing layer develops between the upper part of the flow and the recirculating vortex near the wall.

These results show that the details of the flow between the crests ͑separation, reattachment͒ as shown on Fig. 2, have little effect on the mean flow field at a sufficient distance (Ϸ0.4k s )t o the wall. It is important to remark that the logarithmic plot of all the profiles uses the same mean friction velocity defined by Eq. ͑13͒ and ͑14͒ and not the local value of the shear stress. This latter choice would have led to incoherent results since the shear stress is extremely variable on the wall and may even be equal to zero at the stagnation point A ͑Fig. 2͒.

For the case of an impermeable wall and for the present type of roughness ͑jagged wall profile͒, it may be concluded that the logarithmic law is a good approximation of the mean velocity profile over the whole rough surface. In this case, the location of the virtual wall modeling the rough surface must be taken at the top of the ridged wall.

Suction Effects.

When suction is applied through the regularly spaced slots of the basket, the features of the flow are not very different from the impermeable case, as it was remarked in Section 3.1. In particular, the streamline pattern of the mean flow over the wall crests is not strongly affected by suction. Therefore, modeling suction effects over this rough wall may be tentatively undertaken by considering uniform suction through a fictitious plane porous wall. Considering the results of the impermeable case, it seems appropriate to locate this virtual wall at the crests of the wall.

The problem of suction through a smooth wall was studied by Stevenson ͓4͔ who proposed a mixing-length approach to model the flow. In this model, the mixing length is still supposed proportional to the distance to the wall whereas a source term due to suction is introduced in the momentum budget in the near-wall region. Integration in y of the x-momentum equation results in

2 y 2 ͯ ץu ץy ͯͩ du dy ͪ ϭu* 2 ϪV por u (15)
where V por is related to the suction velocity by Eq. ͑1͒.

Integrating Eq. ͑15͒ in y yields

2/V por ϩ ͑͑1ϩV por ϩ u ϩ ͒ 1/2 Ϫ1 ͒ϭ1/ ln͑ y ϩ ͒ϩC (16)
where V por ϩ ϭϪV por /u*. The constant of integration C is obtained by matching Eq. ͑16͒ to the logarithmic law for smooth impermeable walls when V por ϩ tends to 0. This gives CϭBϭ5.5.

When a rough wall is considered and modeled by a virtual porous plane surface, the same method may be applied, but the right-hand side of Eq. ͑16͒ has to be matched with the logarithmic law corresponding to rough surfaces. In this case, we obtain CϭBϪ⌬Bϭ5.5Ϫ1/ ln͑1ϩ0.3k s ϩ ͒ (17) where ⌬B is the shift of the logarithmic law for a rough surface.

For the very high values of k s ϩ considered in the present study (k s ϩ Ϸ1000), 1 is neglected in the logarithm in Eq. ͑17͒.

So finally, Stevenson's law modified for a rough surface is

u ϩ ͑ y,k s ͒ϭ1/ ln͑ y/k s ͒ϩBϪV por ϩ /4͑1/ ln͑ y/k s ͒ϩB ͒ 2 . ( 18
)
We tried another approach by modifying Wilcox's formulation ͓6͔ for smooth surfaces with suction.

Starting from the experimental work of Andersen et al. ͓7͔, Wilcox obtained a good approximation of the velocity profiles by introducing the following variation of the von Karman's constant:

ϭ/͑1ϩ⌶V por ϩ ͒ with ⌶ϭ3.11ϩ0.61 ln͑ y ϩ ͒ (19)
For very rough impermeable surfaces, roughness effects lead to destruction of the viscous sublayer and the distance to the wall is then correctly scaled with the roughness height k s instead of the viscous sublayer thickness /u*. Assuming the same physical process for rough surfaces with suction, we propose to modify the law of Wilcox by replacing y ϩ by y/k s in Eq. ͑19͒. The mean velocity profile is then described by

u ϩ ͑ y,k s ͒ϭ1/ ln͑ y/k s ͒ϩB ϪV por ϩ /͑3.11ϩ0.61 ln͑ y/k s ͒͒ln͑ y/k s ͒. ( 20
)
The similarity between the modified Stevenson's ͑Eq. ͑18͒͒ and Wilcox's laws ͑Eq. ͑20͒͒ is apparent. Using one of these two laws, it is possible to calculate the friction velocity and the flow rate across the computation domain defined in Section 2.2 if the tangential velocity V t is given.

When the modified Stevenson's law ͑Eq. ͑18͒͒ is used, the following expression is obtained for the friction velocity:

u*/V t ϭ1/ZϩZ/4V por /V t (21)
with Zϭ/ϩB ϭln͑H/k s Ϫ1 ͒.

When the modified Wilcox's law ͑Eq. ͑20͒͒ is considered, u* is given by

u*/V t ϭ͑1ϩ1/V por /V t ͑ 3.11ϩ0.61 ͒ ͒/Z. (22)
The flow rate is easily deduced from the above expressions of the friction velocity u* by integrating, respectively, the Eq. ͑18͒ and ͑20͒ over the domain height, HϪk s . As in the case of impermeable wall, the flow was also computed by using the two-layer zonal model. When suction is applied through the wall, it is necessary to compensate the suction flow rate in order to keep a constant streamwise flow rate and periodic conditions on ͑i͒ and ͑o͒͑ Fig. 1͒. This implies the introduction of a normal velocity at the boundary ͑u͒ of the flow domain. However, as mentioned before, this condition is incompatible with the previous assumption of symmetry on the upper boundary. The solution used in this simulation is then to replace the condition of symmetry on ͑u͒, by a so-called VELOCITY INLET boundary condition. This condition implies the setting of two velocity components and the turbulence quantities on ͑u͒.I n this case, it is still necessary to set the streamwise flow rate Q in order to complete the boundary conditions. In the present numerical procedure, Q was computed with the modified laws of Stevenson and Wilcox ͑Eq. ͑18͒ and ͑20͒͒. Finally, the INLET turbulent quantities introduced are the turbulence length scale and intensity usually observed in a channel flow. The computation starts with a value of 5% for the turbulence intensity on ͑u͒. This quantity is adjusted slightly during the computation by using the result obtained far from the wall.

Figure 5 compares the friction velocity u* obtained by the above computations with the results of the TLZM numerical simu-lations for a range of suction ratios of industrial interest. u* is given by Eqs. ͑13͒ and ͑14͒ and is normalized by the velocity V t at the upper boundary. There is clearly an excellent agreement of TLZM results with the modified Wilcox's law ͑Eq. ͑22͒͒ for low to moderate suction ratios. For the highest suction ratio considered in this study (V por /V t Ͼ0.01), a more pronounced difference of 9% is observed between the two results. The prediction of the modified Stevenson's law overestimates the friction velocity.

It should be noted that roughness effects are stronger than suction effects for the present conditions. In fact, the nondimensionalized friction velocity was computed by using the Blasius law for a smooth wall in a channel of half-width H and center velocity V t . Compared to this smooth wall value ͑Ϸ0.04͒, the roughness effect gives a 80% increase in u*/V t . On the other hand, the highest suction ratio gives only a 23% increase in u*/V t in comparison to the rough impermeable case.

Figures 6 and7 show mean velocity profiles over the crests computed by the different methods for two values of the suction ratio. For a moderate suction ratio (V por /V t ϭ0.0073), Fig. 6 indicates a good agreement of the TLZM results with the velocity profile given by the modified Stevenson's law ͑Eq. ͑18͒͒. Equation ͑20͒ slightly overestimates the mean velocity. The shape of the mean velocity profile as given by TLZM is considerably modified for the highest suction ratio ͑Fig. 7, V por /V t ϭ0.0104). The stiffness of the velocity profile observed for y/k s Ͼ3 is accounted for by neither the modified Wilcox's law nor by the Stevenson's law. In this case, Eq. ͑20͒͑ Wilcox's law͒ compares well with the numerical results for y/k s Ͻ2, but the corresponding velocity profile deviates significantly from the computed one in the outer part of the flow. Equation ͑18͒ does not give better agreement.

It should be remarked that the good agreement of TLZM with Eq. ͑18͒ observed in Fig. 6 deteriorates significantly when u is nondimensionalized by the friction velocity ͑Fig. 8͒. The differences in the velocity profiles reflect the discrepancy between u* as given by the two methods.

Conclusions

In the present study, the flow near a ridged wall was computed by using a two-layer zonal model. The mean velocity profiles obtained by these two-dimensional computations are in very good agreement with the semi-logarithmic law for impermeable rough walls, provided the ordinate y used in this law is counted from a virtual wall located at the crests of the wall. Moreover, the friction coefficient is computed with a good approximation by using the logarithmic law throughout the flow. These results then justify the wall-function approach for this type of two-dimensional rough- ness. The shape of the wall gives rise to a sharp separation at the crests. Moreover the reattachment zone is much smaller than the separation region. This explains most likely why the best choice for the origin of the ordinate y is at the top of the rough surface (y w ϭ0). This result, which is useful for the paper and pulp industry, may also be relevant to transverse-rib roughness, which is used in tubes for enhanced heat transfer, ͓12͔. In this case, the rib spacing will significantly affect the flow and the above conclusion should apply to small values of this parameter, which correspond to limited regions of reattachment, if any.

For three-dimensional roughness, as for Nikuradse's roughness, it is expected that the effects of separation will be weaker and consequently that the virtual origin will be located below the top of the surface (Ϫk s Ͻy w Ͻ0).

The adaptation of the laws with suction from smooth to very rough walls seems also to be satisfactory if we assume that the numerical simulations are accurate. The adapted Wilcox's law gives the best results for the friction velocity and reasonable agreement for the mean velocity profiles. The proposed modification is very simple since it consists to replacing the viscous length scale by the roughness height in the modeling of the von Kar-man's constant. Further confirmation of this model would require comparison of the present results with experiments.

The present work considered a layer of uniform thickness height, H to model the industrial problem of the flow in pressure screens. A more general case would correspond to developing turbulent boundary layers. However, longitudinal advection effects may be considered as small compared to friction effects as in Stevenson's analysis ͑Eq. ͑15͒͒ in such flows. It is then thought that the present results would be only slightly affected by slow x-variations of H.

This method seems to be successful in the present case of twodimensional roughness. It would be interesting to extend the conclusions of the present study in the following directions:

• further work could test the effect of the crest spacing on the position of the virtual wall. • three-dimensional roughness could be considered, but it would require a significantly greater effort owing to the threedimensional calculations, which would be necessary in this case. L ϭ distance between two crests or two slots p ϭ pressure Q ϭ flow rate per unit length in the normal direction q w ϭ flow rate per unit length near the wall ͑Eq. ͑6͒͒ u ϭ mean velocity u* ϭ friction velocity V t ϭ tangential velocity at the upper boundary of the computation domain V p ϭ suction velocity V por ϭ equivalent velocity for a uniformly porous wall

x ϭ abscissa counted from the wall crest y ϭ ordinate counted from the wall crest y w ϭ ordinate of the virtual wall in the logarithmic wall y 1 ϭ distance to the nearest wall Y ϭ ordinate counted from the virtual wall ϭ rate of dissipation of the turbulent kinetic energy ϭ von Karman's constant ϭ fluid density w ϭ wall shear stress ϭ dynamic viscosity ϭ kinematic viscosity

Fig. 1

 1 Fig. 1 Computation domain

Fig. 2

 2 Fig. 2 Pattern of streamlines close to the screen basket. V p ÕV t Ä1Õ3.

Fig. 3 Fig. 4

 34 Fig. 3 Influence of the position of the virtual wall on the mean velocity profile at xÄ0

Fig. 5

 5 Fig. 5 Influence of suction on the friction velocity

Fig. 6 Fig. 7

 67 Fig. 6 Mean velocity profile at xÄ0 for a moderate suction ratio. V por ÕV t Ä0.0073.

Table 1 Influence of the grid resolution on the velocity profile over the wall crest. Differences with the results of case 4 are given in %.

 1

Table 2 Influence of the position of the virtual wall on the fric- tion velocity, axis-velocity, and friction coefficient

 2

		u* ͑m/s͒	V t (H-k)/Q	C f *100
	y w ϭϪk y w ϭ0 Numerical	1.131 1.219 1.21	1.15 1.21 1.20	0.99 1.04 1.04
	computations

Acknowledgments

This research was supported by the company E.&M. LAMORT and by the French Ministry of Education, Research and Technology. Participation of Dariusz ASENDRYCH is gratefully acknowledged.