Modeling of turbulent fluid flow over a rough wall with or with outsuction

The turbulent ﬂow 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 ﬁctitious wall modeling the surface is considered at the top of the rough boundary. The model developed by Wilcox for smooth walls is modiﬁed to account for the surface roughness and gives satisfactory results, especially for the friction coefﬁcient, 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.
2 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: 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 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 ͒͒. ( The dissipation rate is modeled by introducing a second length scale l ϭk 3/2 /l .
The constants proposed by Chen and Patel ͓10͔ are with ϭ0.41, C ϭ0.09

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 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  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 %. 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͒.

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 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 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: with D͑ y,y w ͒ϭ͑ yϪy w ͒͑ln͑͑yϪy w ͒/k s ͒Ϫ1 ͒/ (11) In addition, the skin-friction coefficient was obtained from the above relations by 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͒.
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.
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 where V por is related to the suction velocity by Eq. ͑1͒. Integrating Eq. ͑15͒ in y yields 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 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 .
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: 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: with Zϭ/ϩB ϭln͑H/k s Ϫ1 ͒.
When the modified Wilcox's law ͑Eq. ͑20͒͒ is considered, u* is given by 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 and 7 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.

Nomenclature
B ϭ constant in the logarithmic law ⌬B ϭ roughness function C f ϭ skin-friction coefficient d ϭ width of the slots H ϭ height of the computation domain k ϭ turbulent kinetic energy k s ϭ roughness height 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