12th International Symposium on Turbulence and Shear Flow Phenomena (TSFP12)
Osaka, Japan (Online), July 19-22, 2022

ROUGHNESS EFFECTS MODELING WITH A DOUBLE AVERAGED
NAVIER-STOKES TURBULENCE MODEL

François Chedevergne
DMPE, ONERA, Université de Toulouse
F-31055 Toulouse, France
francois.chedevergne@onera.fr

ABSTRACTwith a mixing length model. Although limited to boundary
layer flows or channel (pipe) flows, this model was improved
over the years and successfully applied in several configurations (McClain et al. (2006); Hosni et al. (1993)). More recently, the model was extended to RANS turbulence models (Stripf et al. (2008); Hanson et al. (2019)), the equations set still being derived from budgets on a control volume. Actually, Aupoix (2016) proved that the discrete element
method can be derived through the use of a volume-averaging
technique (Whitaker (1996)) leading to the Double Averaged
Navier-Stokes (DANS) equations, initially devoted to fluids
motion in porous media.
Until recently, the initial drag force model (Taylor et al.
(1985)) was unquestioned. Yet, Chedevergne & Forooghi(2020) clearly pointed out the limit of the standard drag modeland proposed an improved version. In parallel, Kuwata &Kawaguchi (2019) also found a drag model based on theDarcy-Forchheimer law, providing similar results, and demonstrating the need for a revisited discrete element model. The
main difficulty is then to find a suitable turbulence model,
compliant with the drag model and that reproduced observations made in the roughness sublayer. Kuwata et al. (2019) derived a DANS model based on a Reynolds stress model (Craft& Launder (1996)). It can be considered as the most advanced
model of that type to deal with roughness issues in turbulent
boundary layers. However, in their model, Kuwata et al. neglected the dispersive stresses, ensuing from the use of volume
averaging, and do not consider heat transfer.
Recently, a complete DANS model dedicated to rough surfaces with heat transfer and relying on the k − ω SST modelof Menter (1994) was proposed. The model was tested on aseries of DNS databases and proved its ability to correctly reproduce the velocity profiles in the rouhness sublayer and thelogarithmic region. To drive the model, the concept of representative elementary roughness (RER) was introduced to comply with the discrete element formalism. For academic roughness distributions, the definition of the RER is straightforward,
but a process must be developed to extend the concept to randomly rough surfaces. Additionally, this DANS formulation
was only validated on low-Reynolds configurations, quite far
from most of the industrial applications. Therefore, this paper
proposes a well-defined process to determine the RER and its
application to a high-Reynolds boundary layer rough configuration Squire et al. (2017).

Recent progress of the discrete element method to account for roughness effects in wall bounded flows were obtained with the use of the double averaged Navier-Stokes equations (Chedevergne, 2021). A volume averaging technique
was set up on the Reynolds averaged Navier-Stokes equations
and a new model for the form drag was derived (Chedevergne& Forooghi, 2020). The whole model was tested on a large
panel of DNS configurations but some complementary validation on high Reynolds configurations were still required. Besides, the central question regarding the inputs of the modelfrom a randomly rough surface remained opened. The present
paper aims at answering both issues by applying the developedmodel to a set of high Renyolds number boundary layer experimental configurations (Squire et al., 2017) and by proposing
a methodology of definition of the input parameters. The resulting process provides a representative elementary roughness
from which very satisfactory agreements were obtained on thevelocity profiles for all the tested Reynolds number.

INTRODUCTION #1
The modeling of turbulent boundary layers over roughsurfaces is a long story which dates back to almost a century
with the founding works of Nikuradse (1937) and Schlichting(1937). From these pioneer studies, the concept of equivalentsand grain, whose height is denoted ks thereafter, was used inmost of the modeling approaches. In particular, a large collection of corrections (Hellsten & Laine (1998); Knopp et al.
(2009); Aupoix (2015)) dedicated to RANS turbulence models and relying on this concept were developed. However,
several well-known drawbacks come along with all these corrections such as non-physical behaviors of turbulent quantities in the viscinity of the walls and the strong dependence tohad-oc correlations to determine ks values. Additionnally, specific modifications are required to asses heat transfer since theReynolds analogy do not hold in the equivalent sand grain approach (Aupoix (2016)). Therefore, a more general approach
was also developed to account for roughness effects in turbulent boundary layers.
The idea, introduced by Schlichting (1937), was to use blockage factors to mimic roughness effects. This led to the development of the discrete element method, especially thanks
to studies by Finson (1982) Robertson (1961) Christoph &Pletcher (1983) Lin & Bywater (1982) or Taylor et al. (1985).
Source terms are introduced in the mometum and energy equations to reproduce the viscous and form drags and heat transfer due to the roughness. The standard formulation of Taylor et al. (1985) rely on a simple drag model associated

The DANS model
The model is detailed in Chedevergne (2021) and isbriefly recap thereafter. First, the RANS equations are volume

1

12th International Symposium on Turbulence and Shear Flow Phenomena (TSFP12)
Osaka, Japan (Online), July 19-22, 2022with ck = 2, cω1 = 2.7, cω1 = 0.1 and cω1 = 30. Geometricalparameters, ϕ, the effective element slope, and βm , the mean
blockage, were found to be the most influential on the turbulent properties within the roughness sublayer.
There exists no true model for the dispersive stresses −ũi u˜jthat can be used in the DANS equations (1). However, in thepresent context, the study is restricted to boundary layer configurations where only the −ũṽ has to be provided. According
to Chedevergne (2021) one writes:

averaged following Whitaker (1986). The blockage factor β ,representing the volume fraction opened to the flow, is then introduced in the equations. For a viscous incompressible fluid,
the resulting continuity and momentum equations read:∂ hβ ui i f
=0
∂ xi
!
f
∂ hui i f
1 ∂ hpi f
1 ∂∂ β hui i ff ∂ hui i
+ uj
=−
+
ν
∂t
∂xj
ρ ∂ xi
β ∂xj
∂xj
 D

f
Ef
hu
i
1 ∂
∂
β
∂
ν
i
−
β u0i u0j + β ũi u˜j −
− Fdi
β ∂xj
β ∂xj ∂xj
(1)
where the drag force Fdi is defined by:

Z
p̃
∂ ũi
nk ds
ni ds −
νIfs ρ
I f s ∂ xk
2
= (1 − β )C khui f k hui i f = fd hui i f
πd d

Fdi =

1
Sf

√
∂ k
−ũṽ = Fdisp νt
∂y


y
−1
−
kr
Fdisp = e

(6)

Z

Representative Elementary Roughness(2)

The boundary layer configurations examined in thepresent study are those of Squire et al. (2017) obtained in thelarge facility called High Reynolds Number Boundary Layer
Wind Tunnel (HRNBLWT) at the University of Melbourne.
For all these configurations, the wind tunnel floor was covered
with a sandpaper for which an optical 3D scan was performed
on a 25.4mm × 25.4mm sample. Surface statistics wereextracted from this scan and it was shown (Squire et al., 2017)that the distribution of roughness heights for the sandpaper
surface is approximately normally distributed about the meanelevation y. Figure 1 shows the scan of the sandpaper sample.
To perform computations using the DANS model describedabove, a RER must be defined as explained by Chedevergne(2021). Actually, four parameters are necessary but this figure
can be easily reduced to only two. First, considering only
a RER with circular wall parallel cross-sections, i.e. forwhich relation (3) holds, limits the number of independent
parameters. Second, by assuming an isotropic distribution ofpeaks and troughs on the surface, the spacings Lx et Lz can bereasonably taken as equals. Finally to define the RER, only
the distributions of β and d are required, the spacing Lz being
deduced from relation (3).

For years, the drag coefficient Cd was modeled following Taylor et al. (1985). Even though, some success were obtainedwith this simple model, in combination with the use of a standard turbulence model (McClain et al. (2004); Stripf et al.
(2009); Hanson et al. (2019)), a comprehensive study of thedrag force distribution based on DNS results showed the necessity of more complex modelings. Two models (Kuwata& Kawaguchi (2019); Chedevergne & Forooghi (2020)), having completely different roots, were found to adequately reproduce the drag forces obtained in the DNS. The presentdrag force (2) makes use of the formulation developedby Chedevergne & Forooghi (2020) which relies on the workof Z̆ukauskas (1972) for banks of tubes in channel flows.
Because of this link, only roughness having circular crosssections can be considered. Then,β = 1−

πd 2
4Lx Lz(3)

with d the roughness diameter and Lx and Lz the longitudinal
and transverse spacings between roughness.
The drag coefficient is given by:
Cd = 1.5 (αβ )2 ξ

(4)where α is a parameter controlling the density and ξ is adecreasing function of the Reynolds number Red , based on
the local roughness diameter d. More details can be foundin Chedevergne & Forooghi (2020) and Chedevergne (2021).
f
The modeling of the Reynolds stress − ui v j is made with amodified version of the k − ω SST model. Additional source
terms are included in the turbulent scalar equations to force thedesired behaviors in the roughness sublayer. The main idea is
to aligned these sources terms with the drag force using function fd (2). Thus, sources terms Fk and Fω , respectively applied to k and ω equations, are given by:
fk = ck fd




−1
1−ϕ
fω = cω1 + cω2
1 + cω3 exp
fdϕ
1 − βmFigure 1.

Scan of the sandpaper sample

The blockage factor β and the diameter d can becomputed from the 3D scan for each altitude y using simple
image processing. At a given altitude y, which origin y = 0 is
taken as the bottom of the sandpaper sample, a planar slice isextracted and the image is binarized. Two examples are shownon figures 2a and 2a for y = y + 0.2 mm and y = y respectively.
The blockage factor beta(y) is simply computed as the ratio of(5)

2

12th International Symposium on Turbulence and Shear Flow Phenomena (TSFP12)
Osaka, Japan (Online), July 19-22, 2022

Figure 2.

(a) y = y + 0.2(b) y = y + 0.2(c) y = y(d) y = yImages and contours of slices of the sandpaper sample of figure 1.

Figure 3. Wetted perimeter Pw as a function of the altitude y.
The vertical dashed line shows the separation between k-typeand d-type roughness behaviors located at y = y0 .

Figure 4. Blockage factor β (right) and diameter d with respect to the altitude y. The gray area represents the region
below y0 given in figure 3.

the area of the white region over the total surface, i.e. the ratioof the surface opened to the flow. To compute the diameterd(y), we make use of the definition of the hydraulic diameter4Swd=with Sw the wetted surface and Pw wetted perimeter.
Pw
Area Sw was already computed to access β . Perimeter Pw iscomputed as the sum of all the contour lengths in the image.
Examples of image of contours are plotted on figures 2band 2d. If for y = y + 0.2 mm approximating the surface
roughness distribution with circular cross-section appears tobe quite obvious, this is far from being trivial at altitude y = y.
Final distributions of β and d with respect to y are given on
figure 4, while the perimeter Pw (y) is drawn on figure 3.

instance, McClain et al. (2004) and McClain et al. (2006)
pointed out the presence of dead water zones for dense roughconfigurations and defined the notion of ”mean elevation”
from the melt-down surface. This concept was also invokedby Aupoix (2016) to develop a thermal correction in thecontext of an equivalent sand grain model. More generally,
the issue is related to the type of roughness encountered.
Following Perry et al. (1969), a distinction is made between
k-type and d-type roughness. The discrete element method,
based on the introduction of a blockage factor, is strictly
devoted to k-type roughness and is unable to predict the dragincrease due to d-type roughness. More specifically, the dragmodels Taylor et al. (1985), citeKuwata19b or Chedevergne& Forooghi (2020) are all designed to reproduce the form
drag of k-type roughness. For most of the surface of industrial
interest, k-type and d-type roughness resemble the illustrationof figure 5. For the present configuration, both types ofroughness exist and must be distinguished. Starting from the

Before finalizing the process of building the RER,a recursive issue about the discrete element method forrough surface must be addressed. For randomly roughsurfaces, some parts may not particpate to wall transfer. For3

12th International Symposium on Turbulence and Shear Flow Phenomena (TSFP12)
Osaka, Japan (Online), July 19-22, 2022
out the RER design because the spacings are not constant.
There are indeed computed from the relation (3) and, in thepresent case, Lz (= Lx ) is an increasing function of y for y > y0 .
Figure 5. sketch of typical k-type roughness (left) and d-typeroughness (right).

Acknowledgment
top of sandpaper surface, the wetted perimeter is growing as
y is decreasing. This is a typical behavior encountered withk-type roughness. However, for y below y0 = 0.54 mm, Pw isa monotonic increasing function of y, and roughness should
be considered as d-type roughness. Practically, beta and d arebounded by their lower values β (y0 ) and d(y0 ). The impactof this choice is very limited on the equations set (1) and theclosure relation (4) since the flow is almost ”blocked” in thisregion but this has a strong impact on the mean blockage
value βm that rules the source term in the ω-equation (5).

The author is very grateful to D. Chung, N. Hutchins from
Melbourne University and to M. MacDonald from the University of Auckland for sharing their data and knowledge onroughness issues. Model validation crucially depends on detailed experimental or numerical data and nothing would havebeen done without their support.

REFERENCES
Aupoix, B. 2015 Roughness Corrections for the k − ω Shear
Stress Transport Model: Status and Proposals. Journal of
Fluids Engineering 137, 021202.
Aupoix, B. 2016 Revisiting the discrete element method forpredictions of flows over rough surfaces. Journal of Fluids
Engineering p. 31205.
Chedevergne, F. 2021 A double-averaged navier-stokes kωturbulence model for wall flows over rough surfaces withheat transfer. Journal of Turbulence 0 (0), 1–22.
Chedevergne, F. & Forooghi, P. 2020 On the importance of thedrag coefficient modelling in the double averaged navierstokes equations for prediction of the roughness effects.
Journal of Turbulence 21 (8), 463–482.
Christoph, G.H. & Pletcher, R.H. 1983 Predictions of roughwall skin friction and heat transfer. AIAA Journal 21 (4),509–515.
Craft, T.J. & Launder, B.E. 1996 A Reynolds stress closure
designed for complex geometries. International Journal of
Heat and Fluid Flows 17 (3), 245–254.
Finson, M.L. 1982 A model for rough wall turbulent heating
and skin friction. AIAA Paper 82-0199 20th Aerospace Science Meeting, Orlando, Florida.
Hanson, David R., Kinzel, Michael P. & McClain, Stephen T.
2019 Validation of the discrete element roughness method
for predicting heat transfer on rough surfaces. International
Journal of Heat and Mass Transfer 136, 1217–1232.
Hellsten, A. & Laine, S. 1998 Extension of the k − ω shearstress transport turbulence model for rough-wall flows.
AIAA Journal 36 (9), 1728–1729.
Hosni, M.H., Coleman, H.W., Gardner, J.W. & Taylor, R.P.
1993 Roughness element shape effects on heat transfer andskin friction in rough-wall turbulent boundary layer. International Journal of Heat and Mass Transfer 36 (1), 147–
153.
Knopp, T., Eisfeld, B. & Calvo, J.B. 2009 A new extension
for k − ω turbulence models to account for wall roughness.
International Journal of Heat and Fluid Flow 30, 54–65.
Kuwata, Y. & Kawaguchi, Y. 2019 Direct numerical simulation of turbulence over resolved and modeled rough wallswith irregularly distributed roughness. International Journal of Heat and Fluid Flow 77, 1 – 18.
Kuwata, Y., Suga, K. & Kawaguchi, Y. 2019 An extension of
the second moment closure model for turbulent flows overmacro rough walls. International Journal of Heat and Fluid
Flow 77, 186–201.
Lin, T.C. & Bywater, R.J. 1982 Turbulence models for highspeed, rough-wall boundary layers. AIAA Journal 20 (3),
325–333.
McClain, S.T., Collins, S.P., Hodge, B.K. & Bons, J.P. 2006

Results and discussion
Until now, the model was only tested at low Reynolds
numbers from DNS configurations. The objective is to checkthe ability of the DANS model to reproduce boundary layers at
high Reynolds values and to assess its sensitivity to variation of
the Reynolds number. Among the experiments of Squire et al.
(2017), 6 cases were retained for which dimensionless boundary layer thickness δ + ranges from less than 4 000 to about
30 000. Table 1 provides the main characteristics of these configurations.
To solve the DANS equations (1) we considered the inner region of boundary layer approximations. The convective
terms were neglected and the equations were turned dimensionless using the kinematic viscosity ν and the friction velocity uτ .
The origin of the y coordinate in the experimental data must
be shifted since the reference was set to the roughness crest.
We used the values obtained from the sample of figure 1 torescale the coordinate inducing thus a slight uncertainty on theorigin y = 0. This must not affect the experimental results, except in viscnity of the wall within the roughness sublayer. In
addition, the wall correction parameter ε is taken as ε = h/2,
with h = 0.902 mm the roughness height in compliance withthe definitions given by Squire et al. (2017).
Profiles u+ = u/uτ are plotted on figure 6 for all configurations. An additional smooth profile is depicted on the figure.
The corresponding computation was performed with the standard k − ω SST model Menter (1994).
As the Reynolds number δ + is increasing, the roughness
function δ u+ increases too, shifting down the velocity
profiles. The computed roughness functions are in very good
agreement with the experimental ones with a relative error
limited to 4%, except for the second configuration where it
reached 12%. Some discrepancies occur when the profiles
enter the roughness sublayers. It is hard to say whether itcomes from the different models, the definition of the RER or
some uncertainties in the measurements regarding the locationof the probe. Nevertheless, the remarkable point of all thesecomputations is that they were made from a single RER.
The effect of the Reynolds number variations on this RER isvery well captured by the whole DANS model. Besides, themethodology followed to define the RER on such a complex
topology is efficient, even though reducing the rough surface
to a unique roughness with a circular cross-section does notappear to be self-evident. Remark that it is hard to figure

412th International Symposium on Turbulence and Shear Flow Phenomena (TSFP12)
Osaka, Japan (Online), July 19-22, 2022
Table 1.

Details of the rough wall configurations.

U∞

δ+

uτ

∆u+ (exp.)∆u+ (DANS)± error (%)

4.7

3 970

0.18

2.7

2.81

−4.1

10.1

9 730

0.40

5.2

5.83

−12.1

15.2

14 980

0.61

6.6

6.75

−2.3

20.6

20 160

0.83

7.5

7.79

−3.9

25.4

25 020

1.02

8.2

8.01

2.3

30.4

29 900

1.23

8.6

8.58

0.2

Figure 6. Comparison of streamwise mean velocity profile for all rough walls configurations. Symbols (◦) show the experimental
data and solid lines stand for the DANS computations. Gray color is used for the smooth configuration.

The importance of the mean elevation in predicting skinfriction for flow over closely packed surface roughness.
Journal of Fluids Engineering 128, 579–586.
McClain, S.T., Hodge, B.K. & Bons, J.P. 2004 Predicting skinfriction and heat transfer for turbulent flow over real gas turbine surface roughness using the discrete element method.
Journal of Turbomachinery 126, 259–267.
Menter, F.R. 1994 Two-equation eddy-viscosity turbulencemodels for engineering applications. AIAA Journal 32 (8),1598–1605.
Nikuradse, J. 1937 Laws of flows in rough pipes. Technical
Memorandum 1292. NACA, Washington.
Perry, A.E., Schofield, W.H. & Joubert, P.N. 1969 Rough wall
turbulent boundary layers. Journal of Fluid Mechanics 37,
383.
Robertson, J.A. 1961 Surface resistance as a function of theconcentration and size of roughness elements. PhD thesis,
State University of Iowa.
Schlichting, H. 1937 Experimental investigation of the problem of surface roughness. Technical Memorandum 823.
NACA, Washington.
Squire, D.T., Hutchins, N., Morrill-Winter, C., Schultz, M.P.,

Klewicki, J.C. & Marusic, I. 2017 Applicability of taylor’s hypothesis in rough- and smooth-wall boundary layers.
Journal of Fluid Mechanics 812, 398–417.
Stripf, M., Schulz, A., Bauer, H.J. & Wittig, S. 2009 Extended models for transitional rough wall boundary layerswith heat transfer – Part I: Model formulations. Journal of
Turbomachinery 131, 031016–1 – 031016–10.
Stripf, M., Schulz, A. & Bauer, H.-J. 2008 Modeling of roughwall boundary layer transition and heat transfer on turbine
airfoils. Journal of Turbomachinery 130 (2).
Taylor, R.P., Coleman, H.W. & Hodge, B.K. 1985 Predictionof turbulent rough-wall skin friction using a discrete element approach. Journal of Fluids Engineering 107, 251–
257.
Whitaker, S. 1986 Flows in porous media I: A theoreticalderivation of Darcy’s law. Transport in Porous Media 1, 3–
25.
Whitaker, S. 1996 The forchheimer equation: A theoreticaldevelopment. Transport in Porous Media 25, 27–61.
Z̆ukauskas, A. 1972 Heat Transfer from Tubes in Crossflow,
Advances in Heat Transfer, vol. 8, pp. 93–160. J.P. Hartnett
and T.F. Irvine.

5