NUMERICAL SIMULATIONS AND MODELING OF EQUILIBRIUM AND
NON-EQUILIBRIUM TURBULENT FLOWS ON ROUGH SURFACES

By

Sai Chaitanya Mangavelli

A DISSERTATION

Submitted to
Michigan State University
in partial fulfillment of the requirements
for the degree of

Mechanical Engineering—Doctor of Philosophy

2024

ABSTRACT

Wall-bounded turbulent flows in many engineering applications are subjected to spatial and/or
temporal acceleration and deceleration and surface roughness due to factors such as corrosion.
The first part of the report aims to understand: (1) the coupled effect of roughness and mean-flow
acceleration on wall-bounded turbulence and (2) to how does the turbulence response depends on
the roughness topography, especially for multi-scale roughness. To that end, direct numerical sim-
ulations (DNS) are performed on periodic half-channel flows subjected to an impulse acceleration
of the bulk velocity. Two different roughness topographies: a semi-uniform sand grain roughness
and a multiscale turbine-blade roughness are compared with a baseline smooth-wall case. The
flow developed from a transitionally rough flow at the initial equilibrium state to a fully rough
flow at the new equilibrium state. Comprehensive analyses are carried out for the single-point
statistics of all three components of an instantaneous flow quantity in the presence of a rough
wall: the space-and-phase-average, the spatial heterogeneity of the phase-averaged perturbations
(i.e. form-induced/dispersive fields), and the turbulent fluctuations. It is observed that, while the
smooth-wall case undergoes a reverse transition toward the quasi-laminar state as a results of the
strong acceleration, on the rough walls the reverse transition process is prevented.
Instead, the
turbulence near a rough wall rapidly recovers the equilibrium state, owing to an instantaneous
augmentation of the form-induced perturbations as a result of the acceleration. The roughness
texture plays a role in determining the rate of recovery of the near-wall turbulence to the new steady
state. The effect of the form-induced perturbations on turbulence production, previously thought to
be negligible for equilibrium boundary layers, is found to be significant in a strongly accelerating
flow. In addition, the approximate scaling of the form-induced velocity on the velocity at the edge
of roughness sublayer is identified.

In the second part of the report, the structures of both the turbulent and form-induced fields
are characterized; the dynamic interactions between these two components are identified and
quantitatively compared. A second mechanism through which the form-induced velocities and
stresses indirectly affect the Reynolds stress isotropy through modifying pressure fluctuations is
identified. The form-induced perturbations of the Reynolds stresses, rarely studied in the literature,
are found to play an important role in preventing the reverse transition on rough walls. In addition,
results showed that the structural response of turbulence to the acceleration on the multiscale
turbulent eddies are noticeably
roughness displayed similarities with that on a smooth wall:
elongated in the streamwise direction, indicating local process of reverse transition even on a rough
wall. Evidences suggested that this is due to strong form-induced shear, which linearly stretched
the turbulent eddies instead of non-linearly modifying the eddies through augmented turbulence
production.

The findings above suggest that a key part in the dependence of non-equilibrium turbulence

response on the roughness texture is the form-induced quantities, whose production is known to be
attributed to the work of mean flow against the roughness pressure drag (quantified by the equivalent
sandgrain height, 𝑘 𝑠). Hence, the third part of the report focuses on analyzing and modeling the
dependence of the roughness-sublayer mean-velocity profile 𝑈 (𝑦) and 𝑘 𝑠 on the roughness texture,
based on data-driven machine-learning approaches. Gaussian Process Regression (GPR) and Neural
Networks (NN) are trained on a comprehensive database of half-channel flows over different rough
surfaces consisting of new or existing numerical/experimental datasets of mean velocity profiles
or 𝑘 𝑠. Feature sensitivity analyses using both systematic selection of rough features and principal
component analyses (PCA) consistently identified an optimal (minimal) set of important features
for the prediction of 𝑈 (𝑦), including slope, magnitude and skewness of height variations, spatial
correlation, and element inclination. For 𝑈 (𝑦) prediction, the GPR model is shown capable of
accurate prediction for the majority of test cases with a maximum error (across 𝑦) less than 10%
and a 𝐿2 error less than 0.5%. For the 𝑘 𝑠 prediction, an existing NN method is extended to include
two-point features (i.e. correlation lengths), which is shown to significantly improve the prediction.
A PCA analysis showed that the majority of variation of the roughness topography in currently
available datasets is mostly captured by three principal components, which depend on a set of single-
and two-point features. In addition, the relatively small contribution from element inclination to
the principal components indicates a need for more datasets with systematically varied roughness
inclinations. These results suggest the importance of surface correlation lengths, which are not
usually used in roughness modeling, in both 𝑈 (𝑦) and 𝑘 𝑠 predictions. The sets of important features
identified for 𝑈 (𝑦) and 𝑘 𝑠 are found roughly consistent, due to similar physics (i.e. time-mean flow
separation around roughness elements) determining the two quantities.

The current work provides improved understanding of the turbulence response in non-equilibrium
accelerated flows on rough walls. Scaling relations (of total drag and time-mean velocity fields)
of the roughness-sublayer flow are identified, which will inform future physics-based roughness
models. The flow dependency on the roughness texture is characterized with machine learning tech-
niques and state-of-the-art numerical/experimental datasets. The important roughness features that
play a role in determining non-equilibrium flows response are identified for a wide range of rough
surfaces, which will provide guidance on future generation of new rough surfaces for additional
datasets and improve the physical understanding of individual roughness topographical features.
Findings of the modeling and topography-dependency analyses of the roughness-sublayer flows
carried out herein can be extended to non-equilibrium flows according to the identified scaling of
RSL mean flows. Future directions are identified to complete the understanding of non-equilibrium
wall turbulence in other roughness regimes and separating flows.

ACKNOWLEDGEMENTS

I am very grateful for the support of professors, family, friends, Office of Naval Research and
Michigan State University resources during my doctorate degree. First, I would like to thank my
advisor, Dr. Junlin Yuan, for her constant support and motivation throughout my research and
coursework at the University. The practices and guidelines are helping my future research and
I am very grateful to my committee members, Dr. Giles Brereton, Dr.
industrial endeavours.
Michael Murillo and Dr. Ricardo Mejia-Alvarez for their constant support and guidance. The
course work offered by Dr. Murillo was instrumental in the successful implementation of machine
learning techniques. The teaching assistant experience under Dr. Alvarez will be instrumental in
my future endeavours.

I would also like to express my gratitude to my research colleagues, Saurabh Pargal, Guangchen
Shen and Mostafa Aghaei Jouybari for their constant support and motivation. I would also like to
thank my friends, Karthik Krishnan, Ankit Gautam, Hitesh Ghakkar, Bharath Basti Shenoy, Vishal
Asnani, Abubakr Ayesh, Siddharth Shukla, Ronak Sripal, Abhiroop Ghosh, Nikita, Haritha Naidu,
Sunanda, Souvik, Aarzoo, Yamini, Yash Roy, Tanay, Karthik, Urvi Sethia, Ujjwal, Rohit Reddy,
Kishore Khed, Aashna Joshi, Meghana Jakkam for their constant love and support during my PhD.
Finally, I would like to thank my parents, Venugopala Chary Mangavelli, Prasanna Lakshmi
Mangavelli; uncle, Ranga Chary Mangavelli, brother, Sai Aditya Mangavelli and wife, Aditi Dhad-
waiwale; parents in law, Pradeep Dhadwaiwale, Aparna Dhadwaiwale; cousin, Vijay Chiluveru for
their infinite love and support. Without them, my PhD would have been a distant dream.

iv

.

vi

. vii

.
.
.
.

1
1
1
9

. 11
. 11
. 12
. 17
. 33

. 35
. 35
. 36
. 39
. 40
. 53

. 55
. 55
. 58
. 72
. 81

. 83
. 83
. 85

. 87

. 94

LIST OF TABLES . . . . . . . . . .

LIST OF FIGURES . . . . . . . . . .

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

TABLE OF CONTENTS

CHAPTER 1

CHAPTER 2

.

.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

INTRODUCTION .

.
.
.
1.1 Background and Motivation .
1.2 Literature Review . . . .
.
.
1.3 Research objectives and outline .

.
.
.
.
.
. . .
. .
. .
.
. .
. .
.
.
.
EFFECTS OF SURFACE ROUGHNESS TOPOGRAPHY IN SINGLE
.
POINT STATISTICS OF TRANSIENT CHANNEL FLOWS .
.
.
.
.
.
.
. .
.
.
.
.
.
.
.
.
.

Introduction . . . . . . .
2.1
2.2 Methodology . . . . . .
2.3 Results . . . . . . . . . .
2.4 Conclusions . . . . . . .

.
.
. .
. .
.
.

. .
.
. .
. .

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.

.

.
.
.
.

.
.
.
.
.

CHAPTER 3

.
3.1
Introduction . . . . . . .
.
.
3.2 Problem formulation . .
3.3 Summary of turbulence statistics .
.
3.4 Results . . . . . . . . . .
.
3.5 Conclusions . . . . . . .

EFFECTS OF FORM-INDUCED VELOCITY IN TURBULENT HALF
.
.
CHANNEL FLOWS .
.
.
. .
.
.
. .
.
.
. .
.
. .
.
. .

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.

.
.
.

.
.

.
.

.
.

.
.

.
.

.
.
.
.
.

.
.
.

.

.

CHAPTER 4

DATA-DRIVEN PREDICTIONS OF THE ROUGHNESS SUBLAYER
.
.
FLOW . . . . . .
.
. .
Introduction . . . . . . .
.
.
.

.
.
4.1
.
.
4.2 Prediction of double-averaged mean velocity profiles in RSL .
4.3 Prediction of 𝑘 𝑠
.
.
.
.
4.4 Conclusions and discussions

. .
.
.
.
.
.
.
.
.

. .
.
.
.
.
.
.
.
.

. . . . .

. .
.
.

. .
.
.

.
.
. .

. .
.

.
.
.
.
.

.
.
.
.
.

.
.

.
.

.
.

.
.

.
.

.
.

.
.

.
.

.
.

.
.

.
.

.
.

.
.

.
.

.
.

.
.

.
.

.
.

.
.

.
.

.
.

.
.

.
.

.
.

.
.

.
.

.
.
.

.

.

.
.
.

.

.

.
.
.

.

.

.
.
.

.

.

.
.
.

.

.

.
.
.

.

.

.
.
.

.

.

.
.
.

.

.

.
.
.

.

.

.
.
.

.

.

.
.
.

.

.

.
.
.

.

.

.
.
.

.

.

.
.
.

.

.

CHAPTER 5

5.1 Overall conclusions . . .
5.2 Future work . . . . . . .

CONCLUSIONS AND FUTURE WORK .
.
.
.
.
.
.

.
. .

. .
.

.
.

.
.

.
.

.
.

.
.

.
.

.
.

BIBLIOGRAPHY . . . . . . . . . . .

APPENDIX . . . . . . . . . . . . . .

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

v

LIST OF TABLES

Table 2.1 Geometrical parameters of the sandgrain (S), turbine roughness (T) and the
smooth wall (SM). 𝛿 is channel half height; 𝑅𝑎 is first-order moment of height
fluctuations; 𝑘𝑐 is maximum peak-to-trough height (crest height); 𝑘𝑟𝑚𝑠, 𝑠𝑘
and 𝑘𝑢 are the root-mean-square, skewness and kurtosis of height fluctuations;
ES𝑥𝑖 is the effective slope in 𝑥𝑖; 𝐿𝑥𝑖 is the domain size in 𝑥𝑖. . . . . . . . . . . . . 14

Table 2.2 Simulation parameters in the initial steady state (denoted by subscript ‘0’).
Superscribe ‘+’ indicates normalization by wall units (𝑢𝜏 and viscous length
scale 𝛿𝜈 = 𝜈/𝑢𝜏).

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15

.

.

.

Table 2.3 Simulation parameters in the final steady state (denoted by subscript ‘1’).

. . . . 15

Table 3.1 Simulation parameters at initial (‘0’) and new (‘1’) steady states. ‘+’ indicated
normalization in wall units (i.e. 𝑢𝜏 and viscous length scale 𝛿𝜈 = 𝜈/𝑢𝜏). Δ𝑥𝑖
is the cell size in 𝑥𝑖 direction. 𝑅𝑒𝑏 = 𝑢𝑏𝛿/𝜈 and 𝑅𝑒𝜏 = 𝑢𝜏𝛿/𝜈. 𝛿 is channel
half height.

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36

. .

. .

. .

Table 3.2 Characteristics of the rough surfaces.

. . . . . . . . . . . . . . . . . . . . . . . 38

Table 3.3 Single-point cross-correlation between F and T contours shown in Figs. 3.12

and 3.13. The correlation is defined as ⟨|F | · |T |⟩/(cid:2)⟨F 2⟩1/2 · ⟨T 2⟩1/2(cid:3).

. . . . . 52

Table 4.1 Weighted correlations between individual features with individual principal

components for the 𝑈 (𝑦) datasets.

. . . . . . . . . . . . . . . . . . . . . . . . . 67

Table 4.2 Weighted correlations between individual features with individual principal

components for the 𝑘 𝑠 datasets, when three principal components are used.

. .

. 79

Table 4.3 Weighted correlations between individual features with individual principal

components for the 𝑘 𝑠 datasets, when ten principal components are used.

. . .

. 80

vi

LIST OF FIGURES

. . . 13

. . . . . . . . . . . . . . 17

. . . . . . . . . . . . . . . . . . . . . . . . 18

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14

Figure 2.1 A fraction of sandgrain (a, case S) and turbine-blade (b, case T) rough surfaces

Figure 2.2

Figure 2.3

Figure 2.4

Figure 2.5

.

.

.

.

.

.

), S (

colored by height. (c) Variation of roughness geometry functions along 𝑦.
) and T
(a) Power spectrum of surface height fluctuations in cases S (
) with wavenumbers 𝜅 in 𝑥 (black) and 𝑧 (grey). For S, 𝜆1, 𝜆2 and 𝜆3
(
are the three semiaxes of an ellipsoid grain [1]. (b) PDFs of surface height
fluctuations for cases S and T; normalization is done such that area below
curve is one.
(a) DA streamwise velocity, (b) Reynolds stress components at 𝑡∗ = 0, for SM
) cases. 𝑑 = 0 for SM. In (a), thin dashed lines
(
are 𝑈+ = (𝑦 − 𝑑)+ and 𝑈+ = 1/0.4 log(𝑦 − 𝑑)+ + 5.0.
(a) Temporal variation of 𝐶 𝑓 for all cases; inset magnifies the plot at small 𝑡∗
values. Streamwise DA velocity normalized by instantaneous wall units in (b)
: 𝑈+ = (𝑦−𝑑)+
smooth case and cases with surface S (c) and surface T (d).
and 𝑈+ = 1/0.4 log(𝑦 − 𝑑)+ + 5.0.
(a-c) 𝑈 profiles against 𝑦 in linear scaling, for cases SM (a), S (b) and T (c).
(d) Instantaneous 𝑢𝜏 normalized by the instantaneous value of the velocity
scale at the edge of the roughness sublayer, 𝑈𝑅.

) and T (

. . . . . . . . . . . . . . . . . 20

Figure 2.6 Wall-normal profiles of normal components of the dispersive stress tensor

(normalized by 𝑢2

𝜏0) for (a-c) case S and (d-f) case T.

. . . . . . . . . . . . . . 21

Figure 2.7 Time variation of wall-normal peak value of (a) streamwise, (b) wall-normal,
(c) spanwise components of dispersive stress tensor, normalized by 𝑢2
𝜏0. (d)
Peak value of streamwise component normalized by instantaneous 𝑈2
𝑅.

. . . . . 22

Figure 2.8 Wall-normal profiles of Reynolds stresses, at representative times for S (a,c,e)

and T (b,d,f) cases.

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23

Figure 2.9 Time variation of wall-normal peak value of (a) streamwise, (b)wall-normal
and (c) spanwise normal components of Reynolds stress tensor (normalized
by 𝑢2

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24

. .

.

.

.

.

.

𝜏0). .

Figure 2.10 Variation of anisotropy parameter 𝐾 ∗ at various 𝑡∗ for (a) smooth-wall, (b) S

and (c) T cases.

Steady states.

. . . . . . . . . . . . . . . . . . . . . . . . 25

Figure 2.11 Wall-normal profiles of dispersive (black and grey) and Reynolds (red and

light red) shear stresses, normalized by 𝑢𝜏0 for cases (a) S and (b) T cases.

. . . 26

Figure 2.12 Budget balance for ⟨𝑢′2⟩𝑠 (a-d) and ⟨𝑣′2⟩𝑠 (e-h) in case S. Budget terms are
(gray) 𝜕/𝜕𝑡
(black) 𝑇𝑡,

normalized by 𝑢𝑏0 and 𝛿. Vertical line indicates 𝑦/𝑘𝑐 = 1.
(black) 𝑃𝑠; in (e-f),
Π; In (a-d),
term;
𝑇𝑝.
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28

𝑃 𝑓
𝑠 ;
.
.

. .

.

.

.

vii

Figure 2.13 Budget balance for ⟨𝑢′2⟩𝑠 (a-d) and ⟨𝑣′2⟩𝑠 (e-h) in case T. Budget terms are
(gray) 𝜕/𝜕𝑡
(black) 𝑇𝑡,

normalized by 𝑢𝑏0 and 𝛿. Vertical line indicates 𝑦/𝑘𝑐 = 1.
(black) 𝑃𝑠; in (e-f),
Π; In (a-d),
term;
𝑇𝑝.
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29

𝑃 𝑓
𝑠 ;
.
.

. .

.

.

.

Figure 2.14 Reynolds stress budget terms in the smooth-wall case, normalized by 𝑢𝑏0 and
𝛿. (a) Shear production for ⟨𝑢′2⟩𝑠 and (b) pressure work for ⟨𝑣′2⟩𝑠.
Figure 2.15 Temporal evolution of the positive peak values of Π and peak values of 𝑃 𝑓
𝑠
. . . . 31

(both normalized by 𝑢𝑏0 and 𝛿) in ⟨𝑣′2⟩ budget in cases S (a) and T (b).

. . . . . . . 30

Figure 2.16 Wall-normal profiles of 𝑃′ rms (normalized by 𝑢2
.

at different 𝑡∗.

.

.

.

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31

𝜏0) for cases (a) S and (b) T

Figure 2.17 𝑄-isosurfaces at (a-b) 𝑡∗ = 0.1 and (c-d) 𝑡∗ = 0.5 for (a,c) case S and (b,d)

case T colored by distance to wall and rough surface colored in white.

. . . . . 32

Figure 3.1

(a) Geometries of sandgrain (case S) and turbine-blade-roughness (case T)
surfaces colored by height, with zoomed-in view. (b) Power spectral density
of height fluctuations, with wavenumber 𝜅 in 𝑥 (black) and 𝑧 (gray). (c) Sketch
of a rough surface showing various length definitions.

. . . . . . . . . . . . . . 37

Figure 3.2 Flow statistics of the transient channels: (a) bulk velocity, (b) friction coef-
ficient, and wall-normal maximum values of streamwise and vertical normal
Reynolds stresses (c,d) and form-induced stresses (e,f). (b-f) are reproduced
from [2] with permission.

. . . . . . . . . . . . . . . . . . . . . . . . . . . . 39

Figure 3.3 Contours of 𝑅11 at levels 0.3(0.2)0.9 in Cases SM (a), S (b) and T(c), at four
time instances: 𝑡∗ = 0.0, 0.5, 5, and 22.
Lines fitted based on outermost
points (an example is shown with ◦ in (a)) at contour levels 0.5(0.1)0.9, to
show inclination angles (marked in red). For clarity, a shift of 4.5 and 1.5
units are used for smooth and rough cases, respectively.

. . . . . . . . . . . . . 40

Figure 3.5

Figure 3.4 Temporal variation of streamwise integral length (a,c) and inclination angle
(b,d) based on 𝑅11 in SM (◦) , S (△), T (□) cases. The correlations are
centered at 𝑦+
ref,0 = 10 for SM and 𝑦ref/𝑘𝑐 = 1 for S and T. (a,c) show linear
scalings of 𝑡∗ and (c,d) show logarithmic scalings.
1D premultiplied power spectra of 𝑢′ (a-c) and 𝑣′ (d-f) at 𝑡∗ = 0 (
) and a
, 𝑡∗ values from 0.1 to 0.9 with step 0.1, from 1 to 9
series of times (
with step 1, and from 10 to 22 with step 3, lighter colors at later times), in
smooth-wall case (a,d) at 𝑦+ = 10 and cases S (b,e) and T (c,f) at 𝑦 = 𝑘𝑐.
In (d-f), 𝐸22 at 𝑡∗ = 0 is scaled by a factor of 15 for clarity. In (a-c),
(vertical): highest physically meaningful wavenumber. In (a,d), the profile
corresponding to 𝑡∗ = 5 (discussed in text) is highlighted in red.
Figure 3.6 Variation of 𝑆∗ for (a) SM , (b) S , and (c) T cases in time. In (b) and (c),

. . . . . . . . . . . . . . . . 42

. . . . . . . . 43

vertical dotted lines indicate roughness crest height.

. . . . . . . . . . . . . . . 45

viii

Figure 3.7 Time-scale ratios of (cid:101)

𝑢𝑖 strain rates for (a-b) S and (c-d) T cases at three times:
steady state 𝑡∗ = 0 and long-time state 𝑡∗ = 25 (a,c), as well as a transient
state 𝑡∗ = 0.5 (b,d). Profiles corresponding to 𝑡∗ = 25 are offset upward by 10
units for clarify. Vertical dotted lines indicate roughness crest height.

Figure 3.8 Contours of (cid:101)𝑆12𝑞2/𝜖 at 𝑡∗ = 0.5, at respective peak elevations of 𝑆∗ 𝑓

. . . . . 46
(𝑦)
profiles for cases S (a, at 𝑦/𝑘𝑐 = 1) and T (b, at 𝑦/𝑘𝑐 = 0.75), zoomed in to
the regions shown in the full height maps in (c).

. . . . . . . . . . . . . . . . . 47

12

Figure 3.9 Cumulative distribution functions of local time-scale ratios of (cid:101)𝑆12 shown in

Figures 3.8(a,b).

.

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48

Figure 3.10 (a) RMS of the left-hand-side (red) and right-hand-side (black) of 𝑝′ Pois-
son’s equation (Equation 3.8) at 𝑡∗ = 0.5 (
) for case S.
(b,c) Steady-state RMS of source terms in cases S (b) and T (c), evaluated at
𝑡∗ = 22. All quantities are normalized using 𝑢𝑏0 and 𝛿. . . . . . . . . . . . . . . 49

) and 22 (

Figure 3.11 Temporal variations of RMS of 𝑝′ source terms at their respective wall-normal
peak elevations for the smooth case (a) and at 𝑦 = 0.7𝑘𝑐 for cases S (b) and
T (c).

. . . . . . . . . . . . . . . . . . . . . . . . . . . . .

.

.

.

.

.

.

.

.

.

. 50

Figure 3.12 Contours of F (a-c) and T (d-f) at 𝑦 = 0.7𝑘𝑐 for cases S, at 𝑡∗ = 0.1 (a,d),
0.5 (b,e), and 10 (c,f). Black contour lines mark fluid-solid interfaces. All
quantities are normalized based on 𝑢𝑏0 and 𝛿. 1/8 of the domain is shown.

. . . 51

Figure 3.13 Contours of F (a-c) and T (d-f) at 𝑦 = 0.7𝑘𝑐 for cases T, at 𝑡∗ = 0.1 (a,d),
0.5 (b,e), and 10 (c,f). Black contour lines mark fluid-solid interfaces. All
quantities are normalized based on 𝑢𝑏0 and 𝛿. 1/16 of the domain is shown. . . . 52

Figure 3.14 Temporal variations of the RMS of individual terms in F for cases S (a) and
T (b). Symbols: + sum of all terms, ⃝ 11, △ 12, ▽ 13, ⃝ 21, △ 22, ▽ 23, ⃝
31, △ 32, ▽ 33.

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53

.

.

.

.

Figure 4.1 Wall-normal profiles of double-averaged velocity inside the RSL of spatially
accelerating (a) and suction-blowing (and separating, b) boundary layer flows
[3], showing self-similarity of mean velocity regardless of pressure gradients,
when normalized by 𝑈𝑅. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57

Figure 4.2 Visualization of rough surfaces used in DNS simulations to generate data for
𝑈 (𝑦) prediction. Cases ‘Cxx’ are from [4] with a portion (𝛿 and 0.5𝛿 in 𝑥
and 𝑧, respectively) of the domain shown. Cases ‘Ax’ are newly generated
surfaces used for additional simulations herein, with the complete domain
shown (3𝛿 and 𝛿 in 𝑥 and 𝑧, respectively). Cases A1-A5 are similar to C18,
but with varying 𝐼𝑥. Cases A6 and A7 are similar to C05, with varying 𝐼𝑥.
. . . 59
Figure 4.3 Mean velocity profile as a function of wall normal distance (𝑦) in all cases. . . . 61
Figure 4.4 Clustering of rough surfaces in (𝑘𝑟𝑚𝑠,𝐸 𝑆𝑥,𝑠𝑘 ) space.
. . . . . . . . . . . . . . . 62

ix

Figure 4.5

Illustration of the quantification of sheltering, defined as regions with 𝑢 ≤ 0.
Contour shows 𝑢+. Vertical dotted lines indicate locations used to calculate
various frontal areas: leeward sheltered frontal area (𝐴𝑙), windward frontal
shelter frontal area (𝐴𝑤) and total frontal area (𝐴𝑡).

. . . . . . . . . . . . . . . 62

Figure 4.6 Pair plots of roughness features used as inputs.
Figure 4.7 Dependence of the prediction error norms (𝐿2 in orange and 𝐿∞ in blue) on
the number of features used as ML inputs. For each number of features, the
set of features that yield the minimum error among all combinations of the
same number of features is listed.

. . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . .

. 64

. 65

Figure 4.8 Dependence of 𝑈 (𝑦) prediction error on the number of principal components

used as ML inputs. .

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

. 66

Figure 4.9 PDF of 𝑈 (𝑦) prediction errors in all cases of the GPR algorithm developed
with the optimal set of inputs: (𝑅𝑙𝑥/𝑅𝑎, 𝑠𝑘 , 𝐸 𝑆𝑥, 𝐼𝑥, 𝑘𝑟𝑚𝑠): (a) (100)𝐿∞, (b)
(100)𝐿2.

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67

. .

.

.

.

.

Figure 4.10 Examples of predicted mean velocity profiles compared to their DNS-extracted
profiles: (a) low-error cases with 𝐿∞ = 0.033 (red), 0.027 (blue) and 0.028
(black); (b) relatively high-error cases with 𝐿∞ = 0.075 (red), 0.14 (blue) and
0.09 (black).

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

. .

. .

. 68

Figure 4.11 (a) Occurrences of each individual flow cases with 𝐿∞ > 0.1, when the
(b) DNS-extracted mean velocity
set of optimal input features are used.
distributions of all cases, with the three cases with the highest large-error
occurrences in (a) highlighted in black: C04 (

), C06 (

), C14 (

).

. . . 69

Figure 4.12 Effect of roughness-element inclination on 𝑢 distribution in the RSL: (a) Case
C04 with elements inclined toward the direction of flow and (b) Case C06 with
elements inclined in the opposite direction of flow. Iso-surfaces of 𝑢+ = 1.0
and 𝑢+ = −0.05 are shown in yellow and magenta, respectively, overlaid on
visualized rough surfaces.

. . . . . . . . . . . . . . . . . . . . . . . . . . . . 70

Figure 4.13 Effect of roughness porosity on 𝑢: (a) C05 case with densely distributed
elements, and (b) C18 case with the same 𝑘𝑐 as C05 but a sparser distribution.
Iso-surfaces of 𝑢+ = 1.0 and 𝑢+ = −0.05 are shown in yellow and magenta,
respectively, overlaid on visualized rough surfaces.

. . . . . . . . . . . . . . . 71

Figure 4.14 Mean velocity profiles for (a) the two cases with different 𝐼𝑥 (shown in Fig-
ure 4.12) and (b) the two cases with different porosity values (shown in
Figure 4.13).

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71

. .

.

Figure 4.15 Pair plots of roughness features and 𝑘 𝑠/𝑅𝑎. Orange and blue circles denote

roughness types: bar roughness and non-bar roughness, respectively.

. . . . . . 73

Figure 4.16 Performance of the re-calibrated NN with the same inputs as [4]. (a) Predicted
values (𝑘 𝑠,𝑝𝑟𝑒𝑑) compared to actual values (𝑘 𝑠) of the equivalent sandgrain
height. (b) Percentage error. Blue and red symbols denote numerical and
experimental measurements, respectively. . . . . . . . . . . . . . . . . . . . . . 75

x

Figure 4.17 Performance of the new NN with additional inputs of correlation lengths.
(a) Predicted values (𝑘 𝑠,𝑝𝑟𝑒𝑑) compared to actual values (𝑘 𝑠) of the equiva-
lent sandgrain height. (b) Percentage error. Blue and red symbols denote
numerical and experimental measurements, respectively. . . . . . . . . . . . . . 75

Figure 4.18 Comparison of error distributions of the NN predictions with (red) and without

(blue) correlation lengths as inputs. 𝐿′

∞ ≡ |𝑘 𝑠,𝑝𝑟𝑒𝑑 − 𝑘 𝑠 |/𝑘 𝑠.

. . . . . . . . . . . 76

Figure 4.19 (a) Variation of 𝑘 𝑠 error norm with the number of principal components used
(b) Comparison between the PDFs of 𝑘 𝑠 percentage errors

as NN inputs.
predicted with 3 and 10 principal components.

. . . . . . . . . . . . . . . . . . 77

Figure A.1 Comparison of different ML methods in 𝑈 (𝑦) prediction with different num-

bers of roughness features used as inputs, based on (a) 𝐿∞ and (b) 𝐿2.

. . . . . 94

Figure A.2 Variation of percentage errors (a) and percentage of cases with errors larger

than a threshold (b) with the test-cases (as a percentage of all datasets).

. . . . . 95

xi

CHAPTER 1
INTRODUCTION

1.1 Background and Motivation

Turbulent flows are characterized by three-dimensional (3D), chaotic motions. Such type
of flows are widely seen in applications such as blood flow, flow around turbine blades, ships
and airplanes. Early studies on wall bounded turbulent flows are focused on those bounded
by smooth walls. However, complications like corrosion, debris, cavities roughen the smooth
surfaces. Roughness is known to affect the performance efficiency of machines. For example, the
study by [5] showed a drop of 1.25% in efficiency for a Francis turbine subjected to turbulent flow.
This can have a significant economical impact since the efficiency of a typical turbine is higher than
90%. As another example, [6] studied the effect of biofilm fouling for a ship relevant conditions in
a boundary layer flow. The friction drag is shown to increase by more than three times compared to
a similar smooth-wall flow. Efforts to smooth the surfaces are both time consuming and expensive.
So, an accurate method of modeling of the roughness effects in turbulent flows is needed.

Additionally, turbulent flows in engineering applications usually displayed temporal and/or
spatial variations due to external factors like pressure gradients and surface curvature. These flows
are “non-equilibrium" if a similarity solution of the velocity cannot be found.
It is also seen
that the roughness either augmented the flow separation in APG/decelerration or acted against the
stabilizing effects of FPG/acceleration on the wall-bounded turbulent flow. Moreover, these coupled
effects are not well understood for multiscale roughness, requiring a detailed analysis using reliable
simulations. Hence, this thesis studied the wall-bounded turbulence over multi-scale roughness
subjected to temporal acceleration of mean flow using direct numerical simulations (DNS). The
study aimed to develop a physical understanding of these coupled effects of roughness texture and
unsteadiness in turbulence. The study also aimed at expanding this understanding to wide range
of rough surfaces using data-driven machine learning techniques. This is crucial in developing
models for rough-wall turbulent boundary layers and can be integrated with RANS solvers to offer
a computationally inexpensive alternatives to direct numerical simulations and time consuming and
expensive experiments. The current work is only applicable to rough surfaces that fall in the fully
rough regime (flow pattern in the vicinity of roughness is independent of Reynolds number).

1.2 Literature Review

This section summarizes the literature on wall-bounded turbulent flows subjected to wall rough-
ness and spatial or temporal acceleration or deceleration. First, the focus is given to the effect of
roughness in a fully developed turbulent flows. Next, non-equilibrium flows over roughness are
discussed, followed by a summary on strategies to treat roughness in turbulence modeling.

1

1.2.1 Turbulent flows over rough surfaces

To estimate the amount of hydrodynamic drag generated by roughness, [7] and [8] conducted
experimental studies on turbulent pipe flow with sand papers. The results are summarized in the
form of Moody diagram. Reviews of the effect of roughness on the mean flow and turbulence are
given by [9], [10] and [11] in general, as well as by [12] in atmospheric applications.

The effect of the roughness on the mean velocity is as follows. The roughness alters the

logarithmic velocity profile,

𝑈+ =

1
𝜅 ln(𝑦+) + 5.0 − Δ𝑈+,

(1.1)

where 𝜅 = 0.41 is the von Kármán constant, the subscript + denotes normalization in wall units
(friction velocity, 𝑢𝜏, and the viscous length scale, 𝛿𝜈 = 𝜈/𝑢𝜏) and Δ𝑈+ is the roughness function,
which represents the momentum deficit (or increase of wall shear stress) due to roughness. The
roughness function is dependent on the roughness Reynolds number (𝑘 +, where 𝑘 is a characteristic
roughness height) and the roughness topography. The roughness Reynolds number determines the
extent to which roughness affects the near-wall flow. If 𝑘 ≫ 𝛿𝜈 (or 𝑘 + ≫ 1), then 𝑘 is the only
important length scale near the wall. It follows that

(1.2)

Δ𝑈+ =

𝑠 ) + 𝐶,

1
𝜅 ln(𝑘 +
where 𝐶 is the intercept of logarithmic profile and 𝑘 𝑠 = 𝛼𝑘 (where 𝛼 is a constant) is the equivalent
sandgrain height of the roughness in the limit of high 𝑘 +. 𝑘 𝑠 is defined as the height of the sandgrain
roughness of [7] that would yield the same Δ𝑈+ as the roughness in question, which is usually
very different from sandgrains. The regime of 𝑘 +
𝑠 ≫ 1 is referred to as the “fully rough regime”,
where the flow pattern around the roughness protuberances are fully developed (independent of
𝑠 ∈ (5, 70), they are in the “transitionally rough regime”. In this
𝑘 +). For flows with typically 𝑘 +
region, the roughness function is not a logarithmic function of 𝑘 +
𝑠 ; it depends not only on 𝑘 𝑠 but
also on 𝛿𝜈. When 𝑘 +
𝑠 < 5, the flow is usually similar to that on a smooth wall. This is referred to
as “hydraulically smooth regime”.

The overall effect of roughness on the near wall turbulent statistics in various type of wall
bounded flows are summarized below. A more isotropic Reynolds stress tensors are observed near
the wall in the presence of roughness than on a smooth wall [13]. An intensification in the wall
normal fluctuations (𝑣′) is seen, for example, by [14]. The roughness generally affects only the
near-wall turbulent statistics upon normalizing by wall units [15; 16]. This near-wall region is
known as the roughness sublayer (RSL). According to [17], the roughness sublayer ranges from 3
to 5𝑘 𝑠, depending on the fluid quantity examined. Roughness weakly affects the turbulent structures
outside roughness sublayer [18].

Many studies showed that the roughness effects on near-wall turbulent statistics are dependent
on the roughness topography. For example, [19] studied different roughness topographies and found

2

that the ratio of 𝑘 𝑠/𝑘 depends on the roughness texture. Using different roughness topographies
(rod and mesh) with the same 𝑘 𝑠, [20] showed that the topography affects turbulent energy budget
and velocity spectra. [21] showed that the turbulent length scales on a 2D transverse square bars is
larger than a corresponding 3D roughness. As a result, the roughness sublayer on the 2D roughness
is thicker than that of the 3D one.

Most of the above studies focused on roughness made up of uniformly distributed elements.
However, in realistic flows roughnesses are typically non-uniform and multiscale, including length
scales that may be much larger than the boundary layer thickness.
[22] showed that landscape
topography and bathymetry are similar to fractal geometries, characterized by a power-law decay
in the height power spectra.
[23] and [24] conducted experiments on multiscale roughness in
turbulent boundary layers. They both observed that the large spatial heterogeneity of the rough
surfaces led to formation of secondary flows or streamwise rollers, which extended far into the
boundary layer, accompanied by low and high speed regions. Detailed analysis of the dependency
of the secondary flows on the spanwise spacing of roughness in wall bounded flows can be found
in a large number of studies, for example, [25]. These large-scale motions are also observed by
[26] in turbulence boundary layer exposed to spanwise alternating and varying height sand grit
paper roughness. These time-mean, secondary motions are observed due to alternating upward
and downward fluid motions around high and low roughness peaks, respectively. Turbulent length
scales, Reynolds stresses, and the mean velocity showed strong spatial dependencies.

With advancement in computer processing power, rough-wall flows are studied using direct
numerical simulations (DNS) and large eddy simulations (LES). DNS solves the Navier-Stokes
equations resolving all the scales of motions down to the Kolmogorov scale. LES models the small
scales with the help of a subgrid-scale model and resolved the scales inside and above the inertial
sub-range only. The modelling of roughness in numerical studies is often done using immersed
boundary (IB) method. It allows the use of a simple Cartesian grid for both the fluid and solid
domains, where a localized body force ensures that the velocities at the fluid-solid interface is set
to the solid velocity. Among others, [27] provided a review on the method.

Seminal numerical studies on how roughness topography influences fully developed wall-
bounded flows are summarized below. The early numerical studies by [28] on square rods using
DNS observed the dependency of the wall friction on the roughness topography. Further studies
by [29] on 3D square and circular roughnesses revealed a strong connection between the roughness
function and the wall normal turbulent velocity (𝑣′). The study on different arrays of cube type
roughness by [30] using LES showed that the arrangement affects the details of the organized
turbulent motions but not the types of the motions. [31] and [32] studied influences of the roughness
in a turbulent boundary layer. The roughness significantly affected the shape of near-wall streaky
motions by shortening the length of low-speed streaks. One characteristic effect of roughness is to

3

reduce the integral length scale in the streamwise direction. [33] confirmed this observation in a
turbulent boundary layer subjected to realistic (turbine-blade) roughness. Also, they showed that the
roughness led to a transfer of energy from large scales to small scales via redistribution of turbulent
kinetic energy (TKE), defined as sum of the kinetic energies of turbulent motions in all three
directions. Using DNS of turbulent channel flows, [34] compared uniform sand grain roughness
and multiscale turbine-blade roughness and observed the dependency of conditionally averaged
turbulence motions on the roughness topography. The large wavelengths in multiscale roughness
induced large quasi-streamwise vortices that shared similarities to those on a smooth wall. [35]
studied a turbulent channel flow with a random roughness. Analysis of the right-hand-side of the
turbulent pressure Poisson equation showed the dominance of the 12 component of the rapid source
term around the roughness peak and the 23 component of slow term at the roughness trough. These
terms are considered to arise due to the shear layer and quasi-streamwise vortices, respectively.
This gave an insight into the flow features that are responsible for the pressure fluctuations, which
contributes to TKE redistribution and the Reynolds stress anisotropy.

In the roughness sublayer, the presence of roughness leads to spatial heterogeneity in the time-
or ensemble-averaged variables of the turbulent flows.
[36] introduced a double average (DA)
decomposition to separate the mean-flow spatial heterogeneity from the turbulent fluctuations. The
variable representing this heterogeneity in the velocity is referred to as the “dispersive” or “form-
induced” field. The work of pressure drag against the mean velocity is considered the main source
of this field. The alteration of near-wall turbulence behavior in the presence of roughness can be
partially attributed to the form-induced velocities, which introduced additional turbulent-kinetic-
energy (TKE) production mechanisms [36; 37; 38]. Earlier observations of the form-induced
velocities are made by [39]. DNS on uniform cubical-type roughness showed dominance of the
form-induced stresses over the Reynolds stresses near the roughness. The strength of the fields
is shown to be dependent on the roughness topography [40]. Despite earlier observation that the
form-induced TKE production is negligible [41], [40] showed that the two components of the
total production (the part due to Reynolds stress, 𝑃𝑚, and that due to dispersive stress, 𝑃𝑤) are
individually significant, though their sum may be small in the canonical flows tested. It remains
unclear, however, whether in non-canonical flows such as non-equilibrium wall turbulence, the
turbulence production by form-induced velocities can still be ignored.

1.2.2 Non-equilibrium turbulent flows

Current understanding on non-equilibrium wall bounded flows over smooth and rough-walls is

summarized below. Both spatially and temporally developing flows will be discussed.

1.2.2.1 Smooth-wall flows

Non-equilibrium spatially developing turbulence are summarized below. First, the focus is given
to accelerating flows (or those with favorable pressure gradients, FPGs). Early studies showed that

4

the smooth-wall flow underwent a reverse transition toward a quasi-laminar state due to a delay in
the response of Reynolds stresses to the mean-flow acceleration [42]. The friction coefficient value
and mean velocity profiles may revert to a laminar-like state. [43] conducted experimental studies
on a flat-plate boundary layer with FPG and observed a reduction in the wall-normal momentum
transfer in the quasi-laminar state. This is associated with the sweep events (fourth quadrant) being
replaced by the third quadrant events. Numerical simulations of a spatially developing boundary
layer by [44] show dampening of 𝑝′ fluctuations and subsequent reduction in TKE redistribution in
the quasi-laminar flow. Elongation and stabilization of near-wall streamwise velocity streaks and
vortices are also observed as a consequence of lower energy transfer from 𝑢′ to 𝑣′ and 𝑤′ velocity
components [45; 46]. In addition, the increase in streamwise and spanwise extents of near-wall
turbulent structures relative to the boundary layer thickness is observed by [47] for a turbulent
boundary layer subjected to FPG. The flow recovered quickly once the FPG is removed. On the
other hand, adverse pressure gradient (APG) or free-stream deceleration may lead to reduction in
the wall shear stress and subsequent separation of the flow [48; 49; 50; 51]. In the outer layer, APG
energizes the large-scale structures unlike FPG flows (where the outer-layer response lags behind
the near-wall one) [52; 53]. Similarly, for high-Re, APG turbulent boundary layers, APG is shown
to enhance both the small-scale and large-scale energy containing structures in the outer layer.[54].
Compared to boundary layer flows subjected to spatial acceleration/deceleration, temporally
accelerated/decelerated flows share some common features, when viewed using the concept of the
equivalent convection velocity [55]. For example, [56] studied the response of a transient channel
to an impulse acceleration. The flow underwent a process similar to the reverse transition toward
the quasi-laminar state in a spatially developing flow. Various phases are identified. The initial
phase is characterized by a decrease of friction coefficient and long streaky structures near the wall.
As the flow recovered, turbulence spots are formed near the wall, which disturbed the stability of
long streaks and eventually merged and established a new steady state.

1.2.2.2 Rough-wall flows

Studies of rough-wall spatially accelerating flows showed that the roughness acts against the
effects of FPG on the turbulence. The experimental study by [57] on rough-wall quasi-equilibrium
boundary layers subjected to mild acceleration showed intensified turbulent intensity around rough-
ness elements, which prevented the reverse transition as seen in smooth wall. [44] and [58] studied
flat-plate boundary layers with sand grain roughness exposed to spatial acceleration. The more
isotropic Reynolds stress tensor in the presence of roughness is found, attributed to the addi-
tional production mechanisms due to the form-induced fields. Similar flow responses are seen in
temporally accelerated flows by [59], who studied a periodic channel in response to an impulse
acceleration with a pyramid-type roughness. The rough wall prevented the reverse transition. The
long streaky velocity structures are replaced by strong counter rotating hairpin vortices around

5

roughness elements. These vortices subsequently break down, in a rapid transition to a new equi-
librium state. This behavior is similar to roughness-induced laminar to turbulent transition. The
additional TKE production mechanisms due to the form-induced fields are also studied by [60]
and [61] in oscillatory channel flows with rough walls. The magnitude of the TKE production is
shown to be dependent on the roughness size and the phase of the oscillatory flow and could reach
significant values, unlike in canonical flows (e.g. fully developed channels) where the production
is found negligible. Recent experimental study of a 2D cylindrical rod roughness subjected to wall
suction is done by [62]. A downward shift of the velocity profile in the outer part of boundary
layer is observed. This aided the shedding of the roughness-induced vortical structures and hence
prevented the quasi-laminarization. The effect of sand grain roughness on turbulent boundary
layers subjected to APG are studied by [51]. The roughness is shown to lead to an earlier separation
and a larger separation bubble due to a higher momentum deficit.

From the above review, it is established that the roughness counteracts the FPG effects and,
to some extent, augments the APG effects in wall bounded flows. However, this understanding is
generally limited to distributed roughness elements of the same geometry or a narrow length-scale
range.
It is not clear how would turbulence respond to mean-flow acceleration/deceleration (in
either space or time) on a multiscale roughness. Also, the coupling of the effects of roughness and
unsteadiness on turbulence is not well understood. Physical knowledge of this kind is crucial for
modeling of non-equilibrium turbulence bounded by realistic roughness.

1.2.3 Modelling of rough-wall flows

Often, roughness is not resolved in engineering predictions, since it is small compared to the
length scale of the problem. Commercial codes used in practice often rely on Reynolds-averaged
Navier-Stokes (RANS) simulations and incorporate a one- or two-equation eddy-viscosity models.
The Reynolds stress tensor in the RANS equation requires closure. The approach of the eddy-
viscosity models is to model the Reynolds stress tensor as a function of mean velocity with the
help of an eddy viscosity. In the fully rough regime, the hydrodynamic drag is dependent only on a
single parameter of roughness (𝑘 𝑠). Existing roughness closures in RANS models usually involve
modifying the turbulent scales at the wall based on calibrated functions of the local 𝑘 +

𝑠 value.

A major focus of engineering roughness modeling is on the prediction of 𝑘 𝑠 based on features
of the roughness texture alone (and not on flow quantities). 𝑘 𝑠 depends on not a single geometrical
property of surface roughness, but on a collection of different properties. For example, [19]
studied the dependence of 𝑘 𝑠 on the root mean square (rms) and skewness of the roughness height
fluctuations. In addition, the dependency on the “effective slope” (i.e. the mean slope magnitude
[65] performed DNS of sinusoidal
of the surface) is studied by [63] and [64], among others.
roughnesses in a pipe flow. The roughness function is shown to be dependent on both the effective
slope and average height of roughness. Another study by [66] explored the dependency on three

6

different areas of roughness elements (total roughness reference area, frontal area and windward
wetted surface area). In addition, [64], among others, showed that the separation between the fully
and transitionally rough regimes is also topography dependent. The DNS by [67] on irregular
rough surfaces in the transitionally rough regime showed a dependency of the roughness function
on the rms of roughness height, frontal area ratio, skewness, and streamwise correlation length.
Readers are referred to [68] for a comprehensive review of the dependency of drag on roughness
topography for turbulent flow on canonical, rigid roughness. [69] also offered a comprehensive
summary on the dependency of 𝑘 𝑠 on roughness features. The study concluded that measures of
size, frontal area and coverage of roughness as minimum independent parameters required for its
prediction. However, even for the fully rough regime alone, there is no universal prediction of 𝑘 𝑠
for arbitrary surface topographies.

With the advancement in the processing powers, it is now possible to develop more universal 𝑘 𝑠
models with the help of Machine Learning (ML). [4] trained Deep Neural Networks and Gaussian
Process Regression model to predict 𝑘 𝑠 based on 𝑜(10) geometrical features for 45 different rough
surfaces. [70] trained Convolutional Neural Networks to predict momentum and heat transfer, by
employing a visual mapping between height distribution of rough surfaces and local skin friction
and Nusselt Number. 80 grit blasted rough surfaces are generated using DNS at 𝑅𝑒𝜏 = 180. The
trained algorithm predicted the skin friction and Nusselt number with maximum errors up to 8.1%
and 2.9%, respectively.

Some alternative approaches are also proposed to model the roughness effects on the mean flow
and turbulence. [71] offered a comprehensive review of modelling the force term using discrete
element approach. This force term is prescribed to be proportional to the square of double-averaged
mean velocity in the fully rough regime. Furthermore, the proportionality constant is a function of
𝑘 +
𝑠 and wall-normal distance. [72] modeled the roughness effects using an external forcing term
related to the roughness height and the density of element distribution. Other alternate approaches
include the following. [73] performed DNS of turbulent channel flow on a hexagonally packed
hemispherical roughness. The observed interactions between the large-scale outer layer structures
and the near-wall small scale fluctuations are used to build a model of the roughness effects on the
near-wall, higher-order statistics of velocity fluctuations. [74] proposed an approach that explicitly
modeled the geometry of the roughness wake. The model allowed accurate prediction of the drag
over a cuboidal roughness. Recent study by [75] over a few different rough surfaces developed a
data-informed wall function of velocity inside the roughness sublayer, to be used for equilibrium
and non-equilibrium flow RANS models when the roughness layer is resolved.

However, large uncertainties are observed when the roughness closures, typically developed for
equilibrium flows, are used for non-equilibrium turbulence. [76] compared various RANS closures
with or without wall functions for a strongly accelerating boundary layer with roughness. The wall

7

functions (without a resolved wall region) are found to perform poorly, while closures that resolved
the wall region gave considerable discrepancies from the DNS results.
[77] showed the failure
of engineering turbulence models in predicting the flow separation induced by adverse pressure
gradients and reverse transition induced by favorable pressure gradients. Furthermore, the utility
of 𝑘 𝑠 is limited to equilibrium wall turbulence where a logarithmic velocity layer is present and is
characterized by a slope similar to the smooth-wall value. In non-equilibrium flows such as those
with pressure gradients, the logarithmic layer varies or becomes non-existent and, consequently,
𝑘 𝑠 cannot be defined [78]. Additionally, [78] questioned the existence of outer-layer similarity in
non-equilibrium flows with evidences from past literature supporting both sides of the argument.
The outer-layer similarity states that the roughness affects the near-wall region only. If it does not
apply, modeling of outer-layer effect of roughness would be needed in addition to modeling of the
friction velocity or the near-wall flows. Furthermore, even though an accurate 𝑘 𝑠 model can be
built, it may yield accurate wall shear stress but does not guarantee an accurate prediction of mean
velocity in the roughness sublayer (RSL), which would have significance in, for example, prediction
of the separation point in strong-APG flows where a correct prediction of the RSL mean velocity
is important.

An alternative approach of roughness modeling is to model the RSL mean-velocity profile,
coupled with a RANS closure above the RSL. An example is the approach of [75]. [75] proposed
a composite 𝑈 (𝑦) model inside the RSL. The model consists a zero-velocity zone in the vicinity
of the wall and a upper zone that blends the lower zone with the outer layer based on a cubic
polynomial. The model is calibrated for four different rough surfaces in an open channel flow and
shown to reproduce well the mean velocity profiles. An advantage of the model is that a logarithmic
layer is not required, allowing it to be used in non-equilibrium flows. However, it is not clear how
the mean velocity profile varies as a function of roughness texture, limiting the application of this
model.

One objective of the modeling portion of this thesis is to extend the RSL 𝑈 (𝑦) modeling to a
significantly wider range of rough surfaces, based on machine learning (ML). This type of RSL
modeling not only provides a potential approach for modeling non-equilibrium turbulent boundary
layers, but also allows investigation of the dependencies of the RSL mean flow on various roughness
characteristics. This is of interest for understanding (1) form-induced stress production (as the mean
velocity contribute to one major source term) and (2) momentum and heat/mass transport within
RSL, which is important for, for example, environmental flows including those within urban or
aquatic canopies. Another objective is to extend the deep neural network used by [4] to a wider
range of rough surfaces covering bar-type roughness, pitty surfaces (with negative skewness of
height fluctuations), among others.

8

1.3 Research objectives and outline

The above mentioned limitations of the state-of-art understanding of non-equilibrium rough-

wall turbulence prompted the following research objectives.

1. Identify how roughness and its different textures influence the wall turbulence response to a

mean-flow acceleration, as an example of non-equilibrium boundary layers.

2. Perform detailed analyses of both the turbulent and the form-induced components of the

velocity in the non-equilibrium accelerating turbulence.

3. Based on the physical insights learned above, propose data-driven (machine learning) models
of the roughness sublayer, with a secondary objective to identify the roughness geometrical
features that are most important in influencing the turbulent response.

Roughness-resolved DNS will be used to provide datasets for both the fundamental analyses and
modeling. Periodic half-channel flows, either fully developed or transient accelerated, are used.
Aspects of the transient flows are similar to spatially-developing boundary layers, to which the
transient flows may be considered as a computationally low-cost alternative. In addition, the use of
periodic half-channel flows allows the inclusion of multiscale roughness with length scales much
larger than the boundary layer thickness, up to the channel length/width. The rest of the thesis is
outlined as follows.

1. In Chapter 2, DNS of half-channel flows in response to an abrupt increase in the bulk velocity
is studied. The simulation methodologies are described. Two different roughnesses are
imposed on the wall. Comprehensive analyses are carried out for single-point statistics of
all three components of an instantaneous flow quantity in the presence of a rough wall:
the space-and-phase average, the spatial heterogeneity of the phase-averaged perturbations
form-induced field), and the turbulent fluctuations. The ways how the turbulence
(i.e.
response (to acceleration) is affected by the roughness texture are identified. In addition,
RSL-flow and drag scalings are examined, to provide insights in drag and velocity modeling
in roughness-unresolved predictions.

2. Chapter 3 identifies the importance of the form-induced fields ((cid:101)

𝑢𝑖) in non-equilibrium accel-
𝑢𝑖 may
erated turbulence. This motivated a detailed study in Chapter 3 on other ways how (cid:101)
dynamically affect the turbulence structure and other mechanisms important in Reynolds
stress transport (including production of pressure fluctuations and TKE redistribution), be-
sides the additional TKE production well studied in the literature.

3. The previous chapters establish the importance of the form-induced velocity in affecting
the mean flow and turbulence response in non-equilibrium wall-bounded flows. Recall that
a main component of the production of form-induced kinetic energy is the work of mean
flow 𝑈 (𝑦) against the pressure drag (representative by 𝑘 𝑠) [38]. Hence, one motivation of
Chapter 4 is to identify the dependencies of the RSL 𝑈 (𝑦) and 𝑘 𝑠 on various roughness

9

features. In addition, ML models of the two quantities are developed based on a state-of-the-
art rough-wall-flow database to model the RSL 𝑈 (𝑦) and 𝑘 𝑠 for a wide collection of surfaces.
Machine learning techniques such as principal component analysis (PCA), Gaussian Process
Regression (GPR) and Deep Neural Network (DNN) are used, along with systematic feature
elimination studies, to identify the most important roughness features for RSL 𝑈 (𝑦) and 𝑘 𝑠,
as well as to predict these quantities.

4. Chapter 5 summaries the findings of the thesis, followed by proposed areas of future work.

10

CHAPTER 2
EFFECTS OF SURFACE ROUGHNESS TOPOGRAPHY IN SINGLE POINT 
STATISTICS OF TRANSIENT CHANNEL FLOWS

2.1  Introduction

This  chapter  of  the  thesis  will  be  also  referenced  in  the  text  as  [2].  Wall roughness  plays 
an important role in many engineering and environmental applications and a large body of work 
had been carried out to gain understanding of the dynamics of turbulent flows over rough walls. 
In  addition,  wall-bounded  flows  are  often  subjected  to  external  disturbances  such  as  freestream 
pressure  gradients  and  unsteadiness,  some  of  which  have  been  studied  extensively.  Reviews  of 
the literature in these topics are provided in Chapter 1.

Many  aspects  of  turbulent  flows  in  the  vicinity  of  walls  with  different  kinds  of  roughness 
have been  examined.  However,  understandings  are  generally  limited  to  effects  of  roughness  in 
equilibrium  wall  turbulence  (e.g.  fully  developed,  steady  channels  and  ZPG  boundary  layers). 
The  effects  in  time-  and/or  space-varying  flows  have  received  limited  attention.  This  study 
focused  on  an  important  transient-flow  phenomenon—the  suppression  of  quasi-laminarization 
by surface roughness—with emphasis on the role played by the roughness geometry.  To this end, 
we characterized the response of wall-bounded turbulence to an abrupt mean-flow acceleration in 
a  channel  flow,  using  a  configuration  similar  to  that  of  [59].  Unlike  a  spatially  accelerating 
boundary  layer, where the roughness wavelengths are constrained by the length scale of the mean 
velocity variation, a transient channel-flow simulation can accommodate roughness wavelengths as 
long as the channel itself, because of the use of a periodic streamwise boundary condition.

Direct  numerical  simulations  are  conducted  in  transient  half-height  channel  flows,  in  which 
abrupt  increases  in  mass  flow  rate  are  imposed  on  an  initially  steady-state  flow.  The  behavior 
of turbulent flow over two surfaces with different roughnesses is compared to that in the baseline 
case  of  smooth-wall  flow.  In  this  chapter,  Sec.  3.2  introduced  the  governing  equations  and 
details  of the  numerical  solution,  roughness  geometries,  and  simulation  parameters.  Section  3.4 
compared the results of the smooth-wall and the two rough-wall cases for both the initial steady 
state  (before  the  transient  started,  Sec.  2.3.1)  and  the  recovery  after  the  transient  is  imposed. 
Transient  flow  comparisons  are  carried  out  for  the  mean  velocity  (Sec.  2.3.2),  form-induced 
velocities (Sec. 2.3.3),  the  Reynolds  stresses  (Sec.  2.3.4),  and  the  Reynolds  stress  budgets  (Sec. 
2.3.5).

11

2.2 Methodology

2.2.1 Governing equations

The incompressible flow of a Newtonian fluid is governed by the equations of conservation of

mass and momentum:

𝜕𝑢𝑖
𝜕𝑥𝑖
𝜕𝑢𝑖𝑢 𝑗
𝜕𝑥𝑖

= 0,

= −

𝜕𝑃
𝜕𝑥 𝑗

𝜕𝑢 𝑗
𝜕𝑡

+

+ 𝜈∇2𝑢 𝑗 + 𝐹𝑗 ,

(2.1)

(2.2)

where: 𝑥1, 𝑥2 and 𝑥3 (or 𝑥, 𝑦 and 𝑧) are, respectively, the streamwise, wall-normal and spanwise
directions, and 𝑢 𝑗 (or 𝑢, 𝑣 and 𝑤) are the velocity components in those directions; 𝑃 = 𝑝/𝜌 is the
modified pressure, 𝜌 the density and 𝜈 the kinematic viscosity. The term 𝐹𝑗 in Equation (2.2) is a
body force of the immersed boundary (IB) method used to impose non-slip/penetration boundary
conditions at the rough surface, which is well-resolved by the grid. The method is based on the
volume-of-fluid approach [1]; its detailed implementation in the in-house fluid solver and validation
are described in [38] and [79]. Specifically, the volume fraction occupied by the fluid 𝜙(𝑥, 𝑦, 𝑧)
in each grid cell is calculated during pre-processing. It is used to calculate the force 𝐹1(𝑥, 𝑦, 𝑧, 𝑡)
at each time step, to decrease the momentum in each grid cell by a fraction of (1 − 𝜙). The 𝐹1
values are non-negligible in the boundary cells of roughness, and are negligible elsewhere. The
simulations are performed using a well-validated code that solved the governing equations (2.1)
and (2.2) on a staggered grid using second-order central differences for all terms, second-order
accurate Adams-Bashforth semi-implicit time advancement, and MPI parallelization [80].

In the roughness sublayer, the presence of roughness elements led to spatial heterogeneity in
ensemble-averaged variables. Such heterogeneity is separated from turbulent fluctuations using a
double-averaging (DA) decomposition [36],

𝜃 (𝑥, 𝑦, 𝑧, 𝑡) = ⟨𝜃⟩(𝑦, 𝑡) + (cid:101)𝜃 (𝑥, 𝑦, 𝑧, 𝑡) + 𝜃′(𝑥, 𝑦, 𝑧, 𝑡),

(2.3)

∫

𝐴 𝑓

where 𝜃 is a local, instantaneous flow variable, ⟨𝜃⟩ is the intrinsic spatial average in the (𝑥, 𝑧)-
𝜃𝑑𝐴 (where 𝐴 𝑓 is the area occupied by fluid at elevation 𝑦), 𝜃 denoted
plane, ⟨𝜃⟩(𝑦, 𝑡) = 1/𝐴 𝑓
the ensemble average (the ‘mean’ component), 𝜃′ = 𝜃 − 𝜃 is the turbulent fluctuation from 𝜃, and
(cid:101)𝜃 = 𝜃 − ⟨𝜃⟩ is the dispersive (or form-induced) fluctuation. A superficial spatial averaging is also
used; it is defined as ⟨𝜃⟩𝑠 (𝑦, 𝑡) = 1/𝐴 𝑓
𝜃𝑑𝐴, where 𝐴𝑜 = 𝐿𝑥 𝐿𝑧 is the total horizontal area of the
simulation domain and 𝐿𝑥𝑖 is the domain size in 𝑥𝑖.

𝐴𝑜

∫

The streamwise component of the IBM body force, 𝐹1(𝑥, 𝑦, 𝑧, 𝑡), can be integrated in the wall-
normal direction to yield the total hydrodynamic drag at a (𝑥, 𝑧) location and a time 𝑡, including both
viscous and pressure drags. The wall shear stress 𝜏𝑤 is then calculated as the ensemble-averaged

12

total drag per domain horizontal area:

𝜏𝑤 (𝑡) =

𝜌
𝐿𝑥 𝐿𝑧

∫

V

𝐹1(𝑥, 𝑦, 𝑧, 𝑡) d𝑥 d𝑦 d𝑧,

(2.4)

where V represented the volume of the simulation domain. The friction velocity is then obtained as
𝑢𝜏 (𝑡) = √︁𝜏𝑤/𝜌. This method of 𝜏𝑤 calculation is validated by [38] in a fully developed channel flow
over sandgrain roughness by comparing with the total drag calculated from the mean momentum
balance.

2.2.2 Rough surfaces

Figure 2.1 A fraction of sandgrain (a, case S) and turbine-blade (b, case T) rough surfaces colored
by height. (c) Variation of roughness geometry functions along 𝑦.

Two quite different surfaces roughnesses are considered: one synthetic sandgrain roughness ‘S’
(Fig. 2.1a), and one replicated from a surface scan on a hydraulic turbine blade ‘T’ (Fig. 2.1b). Flow
over these two surfaces is studied in fully-developed half-height channels by [40], who compared
details of their surface-geometric parameters.

Figure 2.1(c) compared the roughness geometry functions between the two roughnesses.
It
increased with 𝑦 and reached its maximum value of 1 at 𝑦 = 𝑘𝑐. The power spectrum and
the probability density functions (PDFs) of local surface height fluctuations 𝑘 (𝑥, 𝑧) − ⟨𝑘 (𝑥, 𝑧)⟩𝑜
are compared in Fig. 2.2. Here, operator ⟨·⟩𝑜 denoted averaging across the wall-parallel plane.
The ‘sandgrain roughness’ surface S is created as a uniform distribution of randomly oriented
ellipsoids of the same geometry [1]. Its roughness height distribution is characterized by dominant
wavelengths much smaller than the channel half-height 𝛿; the spectral decay occured only for scales
smaller than the grain size. In contrast, surface T yielded an approximate power-law decay of the
power spectrum with a slope of −2, indicating similarity to a fractal geometry. The initial turbine
roughness scan is mirrored in both 𝑥 and 𝑧 to produce the final T surface that satisfied periodic
boundary conditions in 𝑥 and 𝑧. When these roughness height distributions are compared with
𝛿, surface S may be considered a homogeneous roughness, whereas surface T is heterogeneous
(with important features at scales larger than 𝛿). The PDFs of height fluctuations show that both

13

roughness geometries are characterized by rather evenly distributed height values with respect to
the average height. Case S shows slight bias toward lower height values.

Figure 2.2 (a) Power spectrum of surface height fluctuations in cases S (
) with
wavenumbers 𝜅 in 𝑥 (black) and 𝑧 (grey). For S, 𝜆1, 𝜆2 and 𝜆3 are the three semiaxes of an ellipsoid
grain [1]. (b) PDFs of surface height fluctuations for cases S and T; normalization is done such that
area below curve is one.

) and T (

Surface 𝑅𝑎/𝛿

𝑘𝑐/𝛿

𝑘𝑟𝑚𝑠/𝑅𝑎

𝑠𝑘

𝑘𝑢

ES𝑥

ES𝑧

(𝐿𝑥, 𝐿𝑧)/𝛿

SM
S
T

0
0.014
0.014

0
0.09
0.13

0
1.05
1.17

0
0
0.48
2.97
0.20 3.49

0

0

0.43 0.44
0.08
0.10

(12, 6)
(12, 6)
(12, 12)

Table 2.1 Geometrical parameters of the sandgrain (S), turbine roughness (T) and the smooth wall
(SM). 𝛿 is channel half height; 𝑅𝑎 is first-order moment of height fluctuations; 𝑘𝑐 is maximum peak-
to-trough height (crest height); 𝑘𝑟𝑚𝑠, 𝑠𝑘 and 𝑘𝑢 are the root-mean-square, skewness and kurtosis of
height fluctuations; ES𝑥𝑖 is the effective slope in 𝑥𝑖; 𝐿𝑥𝑖 is the domain size in 𝑥𝑖.
Geometrical parameters of the two surface roughnesses are compared in Table. 3.2. Both surfaces
shared the same first-order moment of height statistics (𝑅𝑎), and have very similar (with 10%
difference) RMS roughness heights (𝑘𝑟𝑚𝑠), but quite different values of maximum peak-to-trough
height (𝑘𝑐, also called crest height), skewness (𝑠𝑘 ) and effective slope (ES) in 𝑥 or 𝑧. 𝑘𝑐/𝛿 ≤ 0.13
for both cases. A smooth-wall baseline case ‘SM’ is used for purposes of comparison.

2.2.3 Simulation Parameters

Simulation parameters for all cases are listed in Tables 2.2 and 2.3. The initial and final steady
states are denoted using subscripts ‘0’ and ‘1’, respectively. Parameters such as the Reynolds
numbers and roughness crest heights are chosen to ensure that the transient and the final states of
each flow corresponded to the Reynolds-number-independent fully rough regime of flow for both
rough surfaces. The final steady state is reached at different times for different flow statistics; after
the steady state is reached, the value of a flow quantity does not vary in time any more. Note
that during the transient some quantities cross the final steady-state values before the steady state

14

is reached; these quantities include the frictional coefficient, dispersive stresses, and streamwise
Reynolds stress, among others.

Surface 𝑅𝑒𝑏0 𝑅𝑒𝜏0 𝑅𝑎/𝛿

𝑘𝑐/𝛿

SM
S
T

4,000
4,000
4,000

244
320
294

0
0.019
0.019

0
0.08
0.12

, Δ𝑧+)

min

𝑠,∞ (Δ𝑥+, Δ𝑦+
𝑘 +
0
21
10

(4.0, 0.9, 2.0)
(5.0, 1.2, 2.5)
(3.4, 1.1, 3.4)

Table 2.2 Simulation parameters in the initial steady state (denoted by subscript ‘0’). Superscribe
‘+’ indicates normalization by wall units (𝑢𝜏 and viscous length scale 𝛿𝜈 = 𝜈/𝑢𝜏).
, Δ𝑧+)
𝑅𝑒𝜏1 𝑅𝑎/𝛿
626
1023
858

𝑠,∞ (Δ𝑥+, Δ𝑦+
𝑘 +
0
76
23

(9.8, 2.3, 4.8)
(15.3, 3.7, 7.6)
(10.1, 3.3, 10.1)

𝑅𝑒𝑏1
12,000
12,000
12,000

0
0.019
0.019

0
0.08
0.12

SM
S
T

Surface

𝑘𝑐/𝛿

min

Table 2.3 Simulation parameters in the final steady state (denoted by subscript ‘1’).
Here, 𝑘 𝑠,∞ is the corresponding equivalent Nikuradse sandgrain height in the fully rough regime,
defined by [10]. The values of 𝑘 +
𝑠,∞ in Tables 2.2 and 2.3 for both surfaces are deduced from the
ratio 𝑘 𝑠,∞/𝑘𝑐 determined in the study of [64], where a Reynolds number sweep is performed for
each surface to identify the critical roughness Reynolds number 𝑘 +
𝑠,𝑐𝑟𝑖 corresponding to the lower
limit of the fully rough regime. In the final steady state, the values of 𝑘 +
𝑠,∞ are 76 and 23, for S and
T respectively, while the respective values of 𝑘 +
𝑠,𝑐𝑟𝑖 are found to be 60 and 20 [64]—lower than the
steady-state values of 𝑘 +
𝑠,∞—indicating that fully-rough flow had been achieved. Considering that
both rough surfaces share the same value of 𝑅𝑎 and slightly different 𝑅𝑒𝜏, the observation that case
S yields a much higher value of 𝑘 +
𝑠,∞ indicated that it generates significantly more hydrodynamic
drag than case T.

The domain sizes in 𝑥 and 𝑧 are (12𝛿, 6𝛿) for cases SM and S, and (12𝛿, 12𝛿) for case T. For case
T, a larger spanwise domain is required to accommodate larger spanwise roughness wavelengths.
For all cases, 𝐿𝑥 are chosen similar to that used by [56] for a similar smooth-wall transient flow
setup. In the study by [56], a similar 𝑢𝑏1/𝑢𝑏0 ratio of 2.6 is used, but the Reynolds numbers is
much lower (𝑅𝑒𝑏 = 2825 compared to 7404 herein). They found that 𝐿𝑥/𝛿 ≈ 12 is sufficient for
accommodating large-scale motions according to the two-point velocity correlations. For the rough
cases, the near-wall 𝑢′
𝑖 coherence lengths are expected to be shorter than on a smooth wall. It is
validated for cases S and T that the two-point correlation of 𝑢′ at streamwise separation of a half of
the domain size falls below 0.1 at 𝑦/𝛿 = 0.3 for all 𝑡∗.

The numbers of grid points are (768, 273, 768), (768, 263, 768) and (1024, 259, 1024) in
(𝑥, 𝑦, 𝑧) for cases SM, S and T, respectively. The mesh is uniform in 𝑥 and 𝑧 while refined
near the wall in 𝑦. At the final steady state (the more critical one for spatial resolution between
the two steady states), Δ𝑦min < 0.3 for SM and Δ𝑦+ falls between 0.6 and 1.7 in the layer below

15

roughness crest for the rough cases. The stretching ratio of the non-uniform 𝑦 grid is kept below
0.03. The grid sizes normalized by the local Kolmogorov lengthscale 𝜂(𝑦, 𝑡) are, for SM, S, and T
cases respectively, 3-7, 4-9, 2-6 in 𝑥 and 1-2, 2-4, 2-6 in 𝑧, throughout the transient. Δ𝑦/𝜂 values
are much smaller than those in 𝑥 and 𝑧. For a curved channel flow, scales smaller than 15𝜂 are found
to contribute to most of the dissipation [81]. Therefore the spatial resolution herein is believed to
be sufficient to resolve the dissipative scales of turbulent motions. This is also supported by the
small residuals found in Reynolds stress balances as discussed later in Sec. 2.3.5.

Aside from the requirement to resolve the small-scale motions, a rough-wall simulation also
required resolution of the roughness geometry. For the sandgrain roughness, the spatial resolution
of the grain geometry is 8 and 16 points per grain in 𝑥 and 𝑧, which is slightly higher than the
resolution used by [38] (case R1 therein). For the turbine blade roughness, there is no single surface
length scale and so an equivalent ‘roughness element’ size is defined based on the Taylor microscale
of surface height fluctuations 𝜆 as in [64]. 𝜆 is resolved by 10 points in both 𝑥 and 𝑧, same as in
[40]. For the rough-wall cases, 𝑦 = 0 is chosen at the lowest trough elevation of the surface. The
virtual origin, 𝑑—defined as the centroid of the 𝑦-profile of averaged drag distribution [82]—is
used to offset the wall-normal profiles when comparing results between different cases.

The flow configurations used are similar to those of [56]. A strong temporal acceleration is
achieved by imposing a temporally varying streamwise pressure gradient, which is adjusted at each
time step to produce a prescribed ramp-like linear increase in the mass flux from 𝑢𝑏0 to 𝑢𝑏1 (where
𝑢𝑏 is the instantaneous bulk velocity). At time 𝑡∗ ≡ 𝑡𝑢𝑏0/𝛿 = 0, the mass-flow rate is subjected to
a ramp increase (during a short time interval Δ𝑡𝑢𝜏0/𝛿 = 0.08) to achieve 𝑅𝑒𝑏1 = 3𝑅𝑒𝑏0. The flow
is then allowed to develop until a new equilibrium state is reached. This ramp-like acceleration is
strong enough to cause departure from equilibrium above a smooth wall. A half-height channel
is simulated, with symmetric boundary conditions applied at the top boundary, to reduce the
computational cost as compared to a full channel. Periodic boundary conditions are applied in 𝑥
and 𝑧.

For each case, the transient simulations are repeated with different initial conditions (obtained
from a separate fully developed half-height channel simulation at 𝑅𝑒𝜏0 for each case). The transient
simulation is repeated 20 times for SM and S and 40 times for T. The flow statistics at each 𝑡∗
are calculated from ensemble averaging. The adequacy of the number of ensembles from which
averages are constructed is tested by halving it, which led to changes of around 0.5% in the mean
streamwise velocity and less than 4.5% in the Reynolds stresses.

16

2.3 Results

2.3.1

Initial steady state flows

Figure 2.3 (a) DA streamwise velocity, (b) Reynolds stress components at 𝑡∗ = 0, for SM (
S (
𝑈+ = 1/0.4 log(𝑦 − 𝑑)+ + 5.0.

),
) cases. 𝑑 = 0 for SM. In (a), thin dashed lines are 𝑈+ = (𝑦 − 𝑑)+ and

) and T (

At 𝑡∗ = 0, the wall-normal profile of the DA streamwise velocity 𝑈 ≡ ⟨𝑢⟩ in wall units is shown
in Fig. 2.3(a) for all cases. The smooth-wall profile compares well with the classical law-of-the-
wall, except that the logarithmic intercept 𝐵 is measured as 5.5 (assuming the von Kármán constant
is 0.4), within the range of values that are typically observed (within 5% of 5.2 according to [83]).
[84] fitted parameters for the logarithmic law based on available data in the literature and observed
that 𝐵 falls between 5 and 6 for channel flows at displacement thickness Reynolds numbers (𝑅𝑒𝛿∗)
between 300 and 2,000; the present smooth-wall case (SM) gives a 𝑅𝑒𝛿∗ of 715, which falls in this
range. Compared to case SM, the presence of roughness yielded a velocity deficit (Δ𝑈+) in the
logarithmic region as expected. Case S gives a higher Δ𝑈+, consistent with the knowledge that
surface S led to higher drag generation than surface T [40].

Figure 2.3(b) compared all non-zero components of the Reynolds stress tensor among all cases.
The outer-layer similarity is demonstrated for all tensor components. The near-wall layer (where
no matching is observed between rough and smooth cases) corresponds to the RSL, defined herein
as the near-wall region where the dispersive perturbations are significant with respect to the local
mean momentum (see Sec. 2.3.3).
In the RSL, similar observations are made as those by [40]
for the same rough surfaces at a higher Reynolds number of 𝑅𝑒𝜏 = 1000, including a stronger 11
component and weaker 22 and 33 components of the Reynolds stress tensor in Case T compared to
S.

The Reynolds stresses at the channel centerline depart slightly from those at the centerline in a
full channel simulation. This is due to the symmetric boundary condition at the top boundary used
herein to simulate a half-height channel. Specifically, here 𝑣′ is imposed to be zero at 𝑦/𝛿 = 1. The
turbulence kinetic energy is redistributed to 𝑢′ and 𝑤′, which is probably the reason for the slightly

17

elevated 11 and 33 Reynolds stress components at 𝑦/𝛿 = 1. Similar phenomena are observed by
[1].

2.3.2 Temporal variation of ensemble-averaged velocities

Figure 2.4 (a) Temporal variation of 𝐶 𝑓 for all cases; inset magnifies the plot at small 𝑡∗ values.
Streamwise DA velocity normalized by instantaneous wall units in (b) smooth case and cases with
surface S (c) and surface T (d).

: 𝑈+ = (𝑦 − 𝑑)+ and 𝑈+ = 1/0.4 log(𝑦 − 𝑑)+ + 5.0.

First, the friction coefficient is obtained by normalizing the instantaneous 𝜏𝑤 with the initial

bulk velocity,

𝐶 𝑓 (𝑡) =

𝜏𝑤 (𝑡)
𝜌𝑢2
/2
𝑏0

,

(2.5)

and plotted in Fig. 2.4(a). In all cases 𝐶 𝑓 underwent a sudden increase immediately after 𝑡∗ = 0
in response to the impulse increase of 𝑢𝑏. Following this jump, 𝐶 𝑓 decreased rapidly as the flow
starts to recover. In the smooth-wall case the 𝐶 𝑓 (𝑡) curve showed a dip (lower than the long-time
value) with a minimum 𝐶 𝑓 found at around 𝑡∗ = 5. For this flow in the corresponding 𝑡∗ range,
[56] observed elongated low-speed streaks, laminar-like mean velocity profiles (also shown in
Fig. 2.4b), more anisotropic Reynolds stress tensor (also shown later in Fig. 2.9), and increased
TKE production with a frozen pressure strain term (also shown later in Fig. 2.14). These phenomena
are indicative of a reverse transition process toward a quasi-laminar state discussed for spatially
developing turbulence under a strong freestream acceleration (see e.g. [42] and [44]). The new

18

steady state is reached at 𝑡∗ ≈ 10 for the smooth case. In the presence of roughness, there is no
appreciable dip in 𝐶 𝑓 regardless of roughness geometry, indicating that the flows do not undergo
reverse transition. Absence of quasi-laminarization on a rough-walled boundary layer subjected to
strong streamwise acceleration is also reported by [58]. Also notice that 𝐶 𝑓 is lower for case T
compared to case S throughout the transient.

The DA streamwise velocity profiles 𝑈 against 𝑦 (offset by the virtual origin 𝑑), both normalized
by instantaneous wall units, are plotted in Fig. 2.4(b-d). The value of 𝑑 stays almost invariant in
time for both rough cases throughout the transient: 𝑑/𝛿 ≈ 0.045 and 0.058 for cases S and T,
respectively. On the smooth wall, 𝑈+ departs significantly from the logarithmic law-of-the-wall
during reverse transition (𝑡∗ < 5). Initially, the 𝑈 profile is characterized by a thickened viscous
sublayer and an increase of the von Kármán constant 𝜅, as shown at 𝑡∗ = 1. At 𝑡∗ = 3, the profile
became laminar-like, indicating that the flow is in the quasi-laminar state. The logarithmic profile
is recovered after the re-transition. For both rough cases, 𝜅 also increased temporarily. For case
T, this is shown at 𝑡∗ = 1, while for case S, this happened between 𝑡∗ = 0 and 1 and is not shown
in the figure. On T roughness, the increase of 𝜅 is shown at 𝑡∗ = 1. A higher 𝜅 is also observed
for spatially developing boundary layer on the sandgrain roughness [58]. Such phenomenon and
the difference between cases S and T are explained in Sec. 2.3.4. Lastly, we observed that the
logarithmic region recovered faster in case S than in case T. These observations on 𝜅 in cases S and
T will be explained in Sec. 2.3.4.

The 𝑈 profiles with linear 𝑦 axis are compared in Fig. 2.5 in the near-wall region. For all cases,
immediately after the impulse acceleration the near-wall rate of shear (𝜕𝑈/𝜕𝑦) increased rapidly to
values similar to the final steady state. [85] showed that a constant roughness geometry function,
Φ(𝑦) ≡ 𝐴 𝑓 (𝑦)/𝐴𝑜 led to an exponential 𝑈 profile up to the crest height, while a monotonically
increasing Φ(𝑦) as the crest is approached—as in the present cases—lead to a linear velocity
profile. Quasi-linear 𝑈 profiles below 𝑘𝑐 are indeed observed for both cases S and T throughout
the transient. Notice that for case T the sampling of large surface wavelengths is not as sufficient as
the small wavelengths in case S, leading to slight undulation of the quasi-linear 𝑈 (𝑦) profile below
𝑘𝑐 for case T.

In a developing flow, the hydrodynamic drag is more directly connected to the near-wall velocity
than to the bulk velocity [58]. A measure of the near-wall velocity scale is the velocity at the edge
of the RSL, 𝑈𝑅 (𝑡) ≡ 𝑈 (𝑦𝑅, 𝑡), where 𝑦𝑅 is the thickness of the RSL, to be further discussed in
𝑢2⟩1/2/𝑈 < 0.1,
Sec. 2.3.3. The RSL is defined herein for all 𝑡∗ as the near-wall region where ⟨(cid:101)
similar to the approach used by [86] (except for a different threshold value used). Throughout
the course of the transient, 𝑦𝑅/𝑘𝑐 takes values (1.10 ± 0.13) and (1.05 ± 0.10) for cases S and T,
respectively. In other words, the RSL thickness varied weakly in time.

19

Figure 2.5 (a-c) 𝑈 profiles against 𝑦 in linear scaling, for cases SM (a), S (b) and T (c).
(d)
Instantaneous 𝑢𝜏 normalized by the instantaneous value of the velocity scale at the edge of the
roughness sublayer, 𝑈𝑅.
Figure 2.5(d) showed the ratio [𝑢𝜏 (𝑡)/𝑈𝑅 (𝑡)]2 for the two rough cases. After 𝑡∗ ≈ 3, both cases gave
a constant ratio, indicating that the near-wall flows are fully developed even though the 𝑈 profile
shape farther from the wall still varied in time as shown in Fig. 2.5(b-c). For 𝑡∗ < 3, the variation of
𝑢𝜏 (𝑡)/𝑈𝑅 (𝑡)—including the initial spike at 𝑡∗ = 0 and the later dip—showed the non-equilibrium
response of the near-wall 𝑢((cid:174)𝑥, 𝑡). The minimum 𝑢𝜏 (𝑡)/𝑈𝑅 (𝑡) occurred at 𝑡∗ ≈ 0.5 and 1.0 for cases
S and T, respectively. These time instances are much later than the mean mass flow rate ramp up
which ends at 𝑡∗ = 0.08. The observation that the minimum 𝑢𝜏 (𝑡)/𝑈𝑅 (𝑡) is reached at an earlier 𝑡∗
in case S than in case T is consistent with the earlier recovery of Reynolds stresses inside the RSL
discussed in Sec. 2.3.4. The ratio 𝑢𝜏 (𝑡)/𝑈𝑅 (𝑡) depended on individual variations of 𝑢𝜏 and 𝑈𝑅 in
time. The 𝐶 𝑓 (𝑡) variation (Fig. 2.4(a)) described the variation of 𝑢𝜏 magnitude (as 𝑈∞ is constant
for 𝑡 > 0). The 𝐶 𝑓 (𝑡) comparison showed that 𝑢𝜏 recovers more slowly on S than on T. Yet, the 𝑈𝑅
variation is faster on S than on T, contributing to earlier recovery of RSL turbulence on S.

20

2.3.3 Temporal variations of dispersive stresses

𝜏0) for (a-c) case S and (d-f) case T.

Figure 2.6 Wall-normal profiles of normal components of the dispersive stress tensor (normalized
by 𝑢2
In both
Figure 2.6 compared normal dispersive stresses for cases S and T at different times.
cases, significant magnitudes of normal dispersive stresses are confined in a near-wall layer. All
components increased immediately after the impulse acceleration, rapidly reaching the maximum
values before 𝑡∗ = 0.5. These maximum values are considerably higher than the values in the final
steady state. After 𝑡∗ = 0.5, all components rapidly decreased, reaching the steady-state values
at 𝑡∗ = 3. Notice that the stress profile shapes do not change significantly during the transient;
the peak elevations vary only weakly. Figure 2.6 also showed that the RSL thickness, 𝑦𝑅, is not
significantly affected during the transient.

Figure 2.5(a) showed that the peak values of ⟨(cid:101)

𝑢2⟩ slight increased at large times (𝑡∗ > 10)
for both cases. This is probably connected to the development of large scale secondary motions
(high and low momentum regions) induced by spanwise roughness heterogeneity [23; 24], toward a
higher-Reynolds-number new equilibrium state. Our earlier DNS study [87] of open channel flows
at a 𝑅𝑒𝜏 of 1000 showed that the present roughness geometries with similar 𝑘𝑐/𝛿 lead to large
secondary motions of the time-averaged flow. These motions are reflected in the non-zero values
𝑢 in the outer layer. These motions developed even for case S, which does not show dominant
of (cid:101)
spanwise wavelength larger than 𝛿.

The main difference between the two rough cases is that case T is characterized by a higher
𝑤). The temporal variations
dispersive stress anisotropy (with more intense (cid:101)
of peak values of the wall-normal profiles of dispersive stress are compared in Fig. 2.7(a)-(c). It
is shown that the rapid initial increase and subsequent decrease in peak values take place rather

𝑢 and weaker (cid:101)

𝑣 and (cid:101)

21

synchronously for both cases, except for ⟨(cid:101)
𝑡∗ ≈ 3 in case S.

𝑣2⟩, which recovers at 𝑡∗ ≈ 1 in case T, compared to

Figure 2.7 Time variation of wall-normal peak value of (a) streamwise, (b) wall-normal, (c)
spanwise components of dispersive stress tensor, normalized by 𝑢2
𝜏0. (d) Peak value of streamwise
component normalized by instantaneous 𝑈2
𝑅.
For both rough cases, 𝑢𝜏 (𝑡)2/𝑈𝑅 (𝑡)2 in Fig. 2.5(d) and near-wall dispersive stresses in Fig. 2.7(a-c)
approached steady-state values at almost the same time of 𝑡∗ = 3. This may be explained by the
dependence of 𝑢2
𝜏 (or the total drag) on roughness-scale flow separations, which are associated with
𝑢2⟩ peak value
intense (cid:101)
normalized by 𝑈2
𝑅 is plotted in Fig. 2.7(d). One observed that the large variation (72% of the
𝑢2⟩max scaled using 𝑢𝜏0 in Fig. 2.7(a) became much weaker (34% of
maximum value in time) of ⟨(cid:101)
the maximum value in time) when scaled using the instantaneous 𝑈𝑅. The residual time variation
in Fig. 2.7(d) is attributed to the non-equilibrium transient process, and possibly to the transition
from the initial transitionally rough flow to a later fully rough flow.

𝑢𝑖 perturbations. To demonstrate such dependence, the variation of the ⟨(cid:101)

𝑖 ⟩(𝑦) profile shape varied little throughout (Fig. 2.6), normalizing the profile
𝑢2
Given that the ⟨(cid:101)
values using 𝑈2
𝑅 would yield dispersive stress curves that are rather invariant in time. This suggested
𝑖 ⟩(𝑦) is absorbed by 𝑈𝑅 (𝑡)2.
𝑢2
that most of the temporal variations of ⟨(cid:101)

In other words, the mean flow pattern 𝑢𝑖 ( (cid:174)𝑥)/𝑈𝑅 inside the sublayer appears to stay close to

equilibrium during the transient.

22

2.3.4 Temporal variations of Reynolds stresses

Figure 2.8 Wall-normal profiles of Reynolds stresses, at representative times for S (a,c,e) and T
(b,d,f) cases.

For a transient accelerating channel on a smooth wall, [56] showed that ⟨𝑣′2⟩ and ⟨𝑤′2⟩ are
frozen in the initial stage of transient, leading to a quasi-one-dimensional turbulence. Similar
observations on a smooth wall are also made by various researchers from spatially developing
boundary layers under strong acceleration.

The variation of Reynolds stresses on the two rough walls are shown in Fig. 2.8. The Reynolds
stresses are normalized with 𝑢𝜏0, to show variations of the stress magnitude. Unlike the smooth
case described by [56], on the present rough walls all components increase rapidly towards the
new steady state, starting from the near-wall region. During the rapid augmentation of Reynolds
stresses, the peak elevations of all normal components in case S appear to coincide at 1.0𝑘𝑐. This

23

indicated that the new turbulence generation immediately after the acceleration occurs at the crest
height and is closely linked to the dynamics of flow around tall roughness protuberances. Similarly,
for case T, the augmented streamwise Reynolds stress peaked constantly at 0.8𝑘𝑐 for T roughness.
However, for this roughness the peaks of wall-normal and spanwise Reynolds stresses moved from
0.8𝑘𝑐 toward the crest height at late times, probably due to the wide variation in local roughness
peak heights brought by the surface heterogeneity. Between the two rough cases, close to the wall
all normal Reynolds stresses increased faster in the S case than in T. In both cases, the development
of the outer layer profiles lags significantly behind that of the near-wall profiles; the equilibrium
states in the outer layer are achieved at around 𝑡∗ = 22.

Figure 2.9 tracked the time variation of the peak value of Reynolds stress profiles. A dependency
on the roughness geometry is evident. For the smooth case, ⟨𝑣′2⟩ and ⟨𝑤′2⟩ peak values remain
almost frozen for a significant period of time till 𝑡∗ ≈ 5 (corresponding to the stage when 𝐶 𝑓
decreases), while ⟨𝑢′2⟩ peak increased. For 𝑡∗ > 5, ⟨𝑣′2⟩ and ⟨𝑤′2⟩ rapidly increase, marking the
phase of retransition. The streamwise fluctuations overshoot over the steady-state value before
the new equilibrium state is reached. This is associated with both the newly generated TKE
predominantly resting in 𝑢′ and the retransition process. Overshoots (though weaker) are also
observed for other normal Reynolds stress components. Theses observations are consistent with
those described by [56]. A drastic difference brought by the presence of roughness is the rapid
increase of all normal Reynolds stress components starting at the beginning of the transient. All
normal components increase much faster for surface S. The combination of Figs. 2.8 and 2.9

Figure 2.9 Time variation of wall-normal peak value of (a) streamwise, (b)wall-normal and (c)
spanwise normal components of Reynolds stress tensor (normalized by 𝑢2
reveals two stages of Reynolds-stress recovery from an abrupt acceleration on rough walls: (1) an
initial stage of rapid recovery of Reynolds stresses inside the RSL as shown in Fig. 2.8, together
with a rapid variation of 𝑢𝜏/𝑈𝑅 (Fig. 2.5(d)) toward the steady state (till around 𝑡∗ = 1 in case S
and 𝑡∗ = 3 in case T); and (2) a second stage (𝑡∗ > 1 for case S and 𝑡∗ > 3 for T) marked mainly by
the variation of the outer-layer Reynolds stresses.

𝜏0).

Figure 2.10 showed the variation of anisotropy parameter (or the streamwise energy partition

24

parameter) defined by [88] as

𝐾 ∗(𝑦, 𝑡) =

2⟨𝑢′2⟩
⟨𝑣′2⟩ + ⟨𝑤′2⟩

.

(2.6)

As an overall measure of the Reynolds stress anisotropy, 𝐾 ∗ quantified the distribution of TKE
among streamwise fluctuations and those in the other two directions. Profiles of 𝐾 ∗ for all cases
are compared in Fig. 2.10 at several 𝑡∗. For the smooth case in the early stage (𝑡∗ < 3), 𝐾 ∗
increased rapidly near the wall from its initial steady-state profile, indicating flow development
toward quasi-one-dimensional turbulence. On the rough walls, in contrast, the Reynolds stress
anisotropy monotonically decreased, recovering to the long time values well before 𝑡∗ = 3; the
subsequent slow recovery of outer-layer Reynolds stresses does not alter the anisotropy. The
difference between the initial and final steady states for the rough cases is because the flows are not
fully rough at 𝑅𝑒𝜏0. Compared to case S, case T yielded higher Reynolds stress anisotropy for all
𝑡∗.

Figure 2.10 Variation of anisotropy parameter 𝐾 ∗ at various 𝑡∗ for (a) smooth-wall, (b) S and (c) T
cases.

Steady states.

Figure 2.11 compared the dispersive shear and the Reynolds shear stresses. For both roughnesses,
𝑣⟩ is significantly weaker than ⟨𝑢′𝑣′⟩. During the transient,
𝑢
in the initial and final steady states ⟨(cid:101)
(cid:101)
𝑣⟩ becomes significant below 𝑘𝑐, sometimes surpassing the Reynolds shear stress. Thus,
𝑢
however, ⟨(cid:101)
(cid:101)
the dispersive shear stress can be dynamically important for mean-flow response in a developing
flow. We note also that the dispersive shear stress becomes positive (with small magnitudes) for a
brief time at the beginning of the transient. It remains an interesting question for future work as to
what contributes to this phenomenon.

We now explain the earlier observations on the logarithmic slope of 𝑈 profiles shown in Fig. 2.4,
namely, the increase of 𝜅 in rough-wall cases at early 𝑡∗ and the faster recovery of 𝜅 for S than for
T. On a smooth wall, the increase of 𝜅 under a favorable pressure gradient is explained by [89] as a
result of departure of the velocity scale of attached eddies from the scale of 𝑢𝜏, i.e.

𝜅
𝜅𝑜

=

𝑢𝜏
𝑈𝑇

,

25

(2.7)

where 𝜅𝑜 is a universal von Kármán constant measured from canonical flows with zero pressure
gradient, 𝑦𝑐 is the critical 𝑦 at which the near-wall sublayer ‘becomes unstable’ [90], 𝑈𝑇 ≡ √︁𝜏(𝑦𝑐)/𝜌
is the critical velocity scale of the sublayer (i.e. velocity scale of the attached eddies), and 𝜏(𝑦) is the
total shear stress. By fitting smooth-wall data with a widerange of streamwise pressure gradients,

Figure 2.11 Wall-normal profiles of dispersive (black and grey) and Reynolds (red and light red)
shear stresses, normalized by 𝑢𝜏0 for cases (a) S and (b) T cases.
[89] identified 𝑦+
𝑐 to be around 10 to 14. With the presence of rough wall herein, we take 𝑦𝑐 = 𝑘𝑐 for
simplicity; a supporting argument for this choice is that the negative peak of combined Reynolds
and dispersive shear stress is located at around 𝑘𝑐 for both S and T for all 𝑡∗ (Fig. 2.11). This stress
is then considered to set the velocity scale of the attached eddies. After the impulse acceleration,
𝑣⟩(𝑘𝑐) is
𝑢
𝑢𝜏 is higher than its steady-state value (Fig. 2.4(a)), while 𝜏(𝑘𝑐) ≈ 𝜌⟨𝑢′𝑣′⟩(𝑘𝑐) + 𝜌⟨(cid:101)
(cid:101)
lower than its steady state value. Since in the steady state, 𝜏(𝑘𝑐) ≈ 𝑢2
𝜏 [58], this means 𝜅/𝜅𝑜 > 1
according to Eq. (2.7). In addition, the faster recovery of 𝜅 on S roughness is due to the fast recovery
of its 𝜏(𝑘𝑐) to the steady state value—as early as 𝑡∗ ≈ 1.0 (Fig. 2.11(a)); at this point Eq. (2.7)
yields 𝜅/𝜅𝑜 ≈ 1.

Results in this section show that turbulence inside the RSL develops rapidly to the new steady
state, while the outer layer varied slowly. The role of roughness geometry came in by affecting the
rate of response of turbulence inside the RSL and, in doing so, affecting the rate of recovery farther
from the wall.

2.3.5 Temporal variations of Reynolds stress budget balances

The budget equation of the normal Reynolds stresses ⟨𝑢′

𝛼⟩𝑠 (no summation on Greek indices)
in a transient half-height channel flow on a rough wall is extended from the equation for a fully-

𝛼𝑢′

26

developed open-channel flow on a rough wall [37] and reads

𝛼 ⟩𝑠

𝜕⟨𝑢′2
𝜕𝑡

=












𝜕⟨𝑢𝛼⟩
𝜕𝑦
(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)

𝛼𝑣′⟩𝑠

−2⟨𝑢′
(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)

(cid:124)

(cid:123)(cid:122)
𝑃𝑠
(cid:42)

𝜕𝑢′
𝛼
𝜕𝑥𝛼
(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)

−2

𝑃′

(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)

(cid:124)

(cid:123)(cid:122)
Π





+ 𝑃𝑤 + 𝑃𝑚

(cid:32)(cid:32)(cid:32)(cid:32)
(cid:32)(cid:32)(cid:32)(cid:32)

(cid:125)
(cid:123)(cid:122)

𝑃 𝑓

𝑠


𝛼⟩𝑠

(cid:124)

−

(cid:124)

𝜕⟨𝑃′𝑢′
𝜕𝑥𝛼
(cid:123)(cid:122)
𝑇𝑝

(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)

(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)

−2

𝑠
(cid:125)

(cid:124)

(cid:125)

(cid:43)

(cid:28) 𝜕
𝜕𝑥 𝑗
(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)

(cid:157)𝑢′
𝛼𝑢′
𝑢 𝑗
𝛼(cid:101)
(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)

(cid:123)(cid:122)
𝑇𝑤
(cid:42) 𝜕2𝑢′2
𝛼
𝜕𝑥 𝑗 𝜕𝑥 𝑗
(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)

(cid:43)

+𝜈

(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)

(cid:125)

(cid:124)

(cid:123)(cid:122)
𝑇𝜈

(cid:28) 𝜕
𝜕𝑥 𝑗
(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)

(cid:29)

−

(cid:124)

𝑠
(cid:125)

(cid:29)

𝑠
(cid:125)

𝑢′
𝛼𝑢′
𝛼𝑢′
𝑗
(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)

(cid:123)(cid:122)
𝑇𝑡
(cid:42) 𝜕𝑢′
𝛼
𝜕𝑥 𝑗

(cid:123)(cid:122)
𝜖

(cid:43)

𝜕𝑢′
𝛼
𝜕𝑥 𝑗
(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)

𝑠
(cid:125)

−2𝜈
(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)

𝑠
(cid:125)

(cid:124)

,

(2.8)

𝛼𝑢′

𝑗 ⟩⟨𝜕

where the terms on the right-hand-side are, respectively, shear production due to shear of DA
velocity (𝑃𝑠), additional shear production due to shear of the form-induced velocities (𝑃 𝑓
𝑠 , as
the sum of the wake production 𝑃𝑤 = −2⟨(cid:157)𝑢′
𝑢𝛼/𝜕𝑥 𝑗 ⟩𝑠 [36] and a correction term 𝑃𝑚 =
𝑗 𝜕
𝛼𝑢′
(cid:101)
−2⟨𝑢′
𝑢𝛼/𝜕𝑥 𝑗 ⟩𝑠 [37]), transport due to form-induced velocities (𝑇𝑤), turbulent transport (𝑇𝑡),
(cid:101)
viscous transport (𝑇𝜈), pressure strain rate (Π), pressure transport (𝑇𝑝), and viscous dissipation
(𝜖). The interface cells are not included in performing the spatial averaging in the present IBM
framework. If the interface cells are also included, an additional term 2⟨𝐹′
𝛼⟩𝑠—with a small
magnitude of 4% of 𝑃𝑠 peak value appears due to the IBM body force, but the nature of the balance
is not affected. This point is discussed by [40]. For the smooth case, Eq. (2.8) still applies with
the form-induced terms 𝑃 𝑓
𝑠 and 𝑇𝑤 equal to zero. For all cases the residuals of Eq. (2.8) obtained
in the steady states are up to 1% of the peak value of 𝑃𝑠. Here, the budget balance refered to the
left-hand-side and right-hand-side of Eq. (2.8) being equal. The balance applied to all 𝑡∗.

𝛼𝑢′

The budget terms in the smooth case are discussed first. Figure 2.14 showed the sources terms
of 𝑢′ energy (𝑃𝑠,11) and of 𝑣′ energy (pressure work, Π22 + 𝑇𝑝22) at several time instances: (1) the
initial state, (2) 𝑡∗ = 3 (quasi-laminarization stage) and (3) 𝑡∗ = 5 (near the onset of retransition).
Figure 2.14(a) showed that 𝑃𝑠,11 increased steadily throughout the quasi-laminarization stage,
leading to intensification of 𝑢′ (i.e. 𝜕⟨𝑢′2⟩𝑠/𝜕𝑡 > 0 ) immediately after 𝑡∗ = 0. In comparison,
the pressure work remained almost frozen in the quasi-laminarization stage. Consequently, the 𝑣′
intensity remained around the initial level in this stage, as shown at 𝑡∗ = 3 in Fig. 2.9. A similar
slow response of pressure work for 𝑤′ energy is also observed (not shown). The delay in increase
of 𝑣′ and 𝑤′ intensities together with the immediate augmentation of 𝑢′ led to the much increased
Reynolds stress anisotropy and quasi-one-dimensional turbulence during quasi-laminarization.

27

Figure 2.12 Budget balance for ⟨𝑢′2⟩𝑠 (a-d) and ⟨𝑣′2⟩𝑠 (e-h) in case S. Budget terms are normalized
by 𝑢𝑏0 and 𝛿. Vertical line indicates 𝑦/𝑘𝑐 = 1.
Π; In (a-d),
𝑇𝑝.

(black) 𝑃𝑠; in (e-f),

(gray) 𝜕/𝜕𝑡 term;

(black) 𝑇𝑡,

𝑃 𝑓
𝑠 ;

28

Figure 2.13 Budget balance for ⟨𝑢′2⟩𝑠 (a-d) and ⟨𝑣′2⟩𝑠 (e-h) in case T. Budget terms are normalized
by 𝑢𝑏0 and 𝛿. Vertical line indicates 𝑦/𝑘𝑐 = 1.
Π; In (a-d),
𝑇𝑝.

(black) 𝑃𝑠; in (e-f),

(gray) 𝜕/𝜕𝑡 term;

(black) 𝑇𝑡,

𝑃 𝑓
𝑠 ;

29

In the presence of roughness, a different balance from that in the smooth case is observed
immediately after 𝑡∗ = 0. Figure 2.12 showed the budget for 𝑢′ and 𝑣′ energy for case S, at
four characteristics time instances: 𝑡∗ = 0, 𝑡∗ = 0.1 and 0.5 (two instances in the stage of rapid
development of RSL turbulence), and 𝑡∗ = 1 (after the RSL is fully developed). Comparison with
the balances at 𝑡∗ = 25 (not shown) suggested that the near-wall balance recovers to the steady state
at 𝑡∗ ≈ 1.0. For the streamwise Reynolds stress, 𝑃 𝑓
𝑠 took positive values with significant magnitudes
in both the transient and the final steady state, though the magnitudes are smaller compared to the
magnitudes of 𝑃𝑠 and Π in most of the RSL. The time-rate-of-change term remained positive
(indicating augmentation of 𝑢′) till 𝑡∗ ≈ 1.0. For the wall-normal Reynolds stress, 𝑃 𝑓
𝑠 also took
significant magnitudes. At 𝑡∗ = 0.1, 𝑃 𝑓
𝑠 is of similar magnitude as Π. At this time, the time-rate-of-
change of ⟨𝑣′2⟩ is positive and peaks roughly at the elevation of the 𝑃 𝑓
𝑠 peak. At later times (shown
at 𝑡∗ = 0.5 as an example), the Π magnitude overtook that of 𝑃 𝑓
𝑠 . The steady state is recovered at
𝑡∗ = 1 for the near-wall balance of 𝑣′ energy also. The subsequent recovery of outer-layer budget
terms did not modify the nature of the balance.

Figure 2.14 Reynolds stress budget terms in the smooth-wall case, normalized by 𝑢𝑏0 and 𝛿. (a)
Shear production for ⟨𝑢′2⟩𝑠 and (b) pressure work for ⟨𝑣′2⟩𝑠.

The same budget terms are shown in Fig. 2.13 for case T, also at four characteristic time instances.
The first three of the four 𝑡∗ instances are chosen to be the same as in Fig. 2.12, while the last time
is chosen at a slightly later time of 𝑡∗ = 3, at which RSL became fully developed for case T. For the
streamwise Reynolds stress, 𝑃 𝑓
𝑠 showed a rather complex pattern, as its sign changed for different 𝑦
and 𝑡∗. Recall that 𝑃 𝑓
𝑢𝑖 kinetic energy
(or ‘wake kinetic energy’, WKE, termed by [36]) at scales larger than turbulent eddies to the TKE
by shear production. Such complex dependence of the sign of 𝑃 𝑓
𝑠 on the wall-normal distance
and time may suggest a competition between TKE-to-WKE and WKE-to-TKE transfers. On a
𝑢𝑖 scales exist,
multiscale roughness, it is possible that both very large 𝑢′ scales and very small (cid:101)
contributing to local energy transfer from TKE to WKE. For the wall-normal Reynolds stress, case

𝑠 represented the amount of kinetic energy transferred from (cid:101)

30

𝑠 in the initial stage of RSL development, although 𝑃 𝑓
T also yielded large magnitudes of 𝑃 𝑓
𝑠 ≪ Π
in the steady states. In both cases, 𝑃 𝑓
𝑠 remained positive for all 𝑡∗ and all elevations for the ⟨𝑣′2⟩
𝑣 to 𝑣′ motions. This is probably because intense 𝑣′ values
budget, indicating energy transfer from (cid:101)
(associated with small-scale active motions) occur typically at a scale smaller than the characteristic
lengthscale of the time-mean velocity variation.

In summary, the similarities between the two rough cases included: (1) the immediate augmen-
𝑠 , especially important for ⟨𝑣′2⟩ production; (2) the subsequent increase of Π, which
𝑠 weakened in time; and that (3) the ⟨𝑢′2⟩ budget

tation of 𝑃 𝑓
later became the dominant source of ⟨𝑣′2⟩ as 𝑃 𝑓
balance is also modified by 𝑃 𝑓

𝑠 , which took significant magnitudes during the transient.

Figure 2.15 Temporal evolution of the positive peak values of Π and peak values of 𝑃 𝑓
normalized by 𝑢𝑏0 and 𝛿) in ⟨𝑣′2⟩ budget in cases S (a) and T (b).

𝑠 (both

Figure 2.16 Wall-normal profiles of 𝑃′ rms (normalized by 𝑢2
𝑡∗.

𝜏0) for cases (a) S and (b) T at different

The main effect of the roughness geometry appeared to be reflected in different rates of recovery of
the budget balance. Here, we focused on the ⟨𝑣′2⟩ budget, as the temporal variation of its balance
is more significant than that for ⟨𝑢′2⟩. Figure 2.15 compared the temporal variation of the peak

31

values of Π and 𝑃 𝑓
𝑠 for cases S and T. The sign of Π changes inside the RSL; the positive peak of
Π is shown here to better represent the nature of this term being an overall source for ⟨𝑣′2⟩. It is
shown that for both surfaces 𝑃 𝑓
𝑠 peaks at a similar 𝑡∗ of around 0.3-0.4. This is because that 𝑃 𝑓
𝑠
𝑢𝑖/𝜕𝑥 𝑗 , and that the time variations of the dispersive stresses between cases S and T
depended on 𝜕
(cid:101)
𝑢𝑖 on 𝑈𝑅. In contrast, Π increases and
are almost synchronous due to the approximate scaling of (cid:101)
peaks at 𝑡∗ ≈ 0.5-0.7 for case S, while it peaked much later—at 𝑡∗ ≈ 3—for case T. To provide
further explanation to the different timing of Π term recovery, Fig. 2.16 compared the time variation
of the pressure fluctuation RMS profiles between S and T cases; a faster recovery of near-wall 𝑃′
intensity is seen for case S. Specifically, the 𝑃′ RMS values below 𝑘𝑐 in case S almost reached their
steady state values at 𝑡∗ = 0.5, while in case T they continued to increase in large time.

Figure 2.17 𝑄-isosurfaces at (a-b) 𝑡∗ = 0.1 and (c-d) 𝑡∗ = 0.5 for (a,c) case S and (b,d) case T
colored by distance to wall and rough surface colored in white.

To provide insights as to why the 𝑃′ intensity varied differently between cases S and T, Fig. 2.17
compared the coherent vortical turbulent motions visualized using the second-invariant of the
velocity gradient tensor,

𝜕𝑢𝑖
1
𝜕𝑥 𝑗
2
The generation of new turbulence of small scales during 0.1 < 𝑡∗ < 0.5 are evident. These motions
are located at 𝑦 of 0.1-0.2𝛿 (in the vicinity of 𝑘𝑐). On surface S, these newly generated eddies
are widespread in the horizontal plane, as the roughness protuberances are closely packed. On

𝜕𝑢 𝑗
𝜕𝑥𝑖

𝑄 = −

(2.9)

.

32

surface T, in comparison, the distribution of these motions are more spotty, associated with more
sparsely distributed high-slope roughness peaks. It is conjectured that the large wavelengths (with
mild slopes) of surface T did not lead to flow separation or strong 𝜕
𝑢𝑖/𝜕𝑥 𝑗 that generated 𝑣′ energy
(cid:101)
through 𝑃 𝑓
𝑠 , while the sparsely distributed small-wavelength peaks produced turbulence, which then
slowly spread to cover the whole horizontal domain. As the 𝑃′ distribution is related to vortical
motions, the 𝑃′ RMS developed more gradually in case T.
2.4 Conclusions

Direct numerical simulations of the impulse response of half-height channel flows on a smooth
wall and two rough surfaces to a step increase of bulk velocity are carried out. Two roughness
geometries are included: a densely distributed sandgrain roughness (case S) and a multiscale
turbine-blade roughness (case T). Both cases shared the same mean roughness height and developed
from a transitionally rough flow to a fully rough one during the transient. Different flow responses
between the rough and smooth cases, as well as between the two rough cases, are discussed. As seen
in the existing literature of temporally and spatially accelerated flows, turbulence over a smooth
wall underwent a quasi-laminarization process, characterized by high Reynolds stress anisotropy
caused by a delayed response of the TKE redistribution to the acceleration. These phenomena are
not observed in the presence of roughness.

On rough walls, the roughness sublayer layer (RSL) thickness stayed almost constant during
the transient. Inside the RSL, the temporal variations of the form-induced fluctuation are largely
absorbed when normalized by the instantaneous mean velocity at the edge of the RSL (i.e., 𝑈𝑅).
This indicated that the ensemble-averaged flow (𝑢/𝑈𝑅) in the sublayer stayed close to equilibrium
during the transient.

Two stages of the turbulence recovery from an abrupt acceleration on a rough wall toward the
new steady state are observed, including: (1) a short, initial stage of rapid recovery of hydrodynamic
drag and increase of RSL turbulent intensities till their early recovery to the long-time values, and
(2) a slower second stage marked mainly by the variations of the outer-layer Reynolds stresses. In
the second stage, the Reynolds stress balance is similar to the steady state, while during the first
stage a fundamental change of Reynolds stress balance is observed. In particular, the magnitude of
the form-induced production 𝑃 𝑓
𝑠 became much larger than in the steady state, strikingly so for ⟨𝑣′2⟩.
The result is a rapid augmentation of TKE and a decrease of Reynolds stress anisotropy, manifested
by the generation of small-scale new turbulence in regions of intense form-induced shear. Also, in
the first stage the total shear stress at roughness crest is significantly lower than the wall shear stress,
leading to a reduced logarithmic slope of the mean velocity. The roughness geometry affects the
rate of recovery of the RSL turbulence. A slower recovery of TKE redistribution (due to a slower
amplification of pressure fluctuations) is seen in case T. This is associated with sparsely distributed
small-wavelength roughness peaks, producing new turbulence which took time to spread and cover

33

the whole domain.

This work provides detailed explanations on how roughness suppressed reverse transition in
a non-equilibrium accelerated turbulence. Results indicated that, unlike fully-developed wall
bounded turbulence where outer-layer similarity applied and the roughness sublayer flow affected
the value of 𝑢𝜏 only, in the present non-equilibrium transient flows the sublayer details played a
more important role by determining the rate of flow recovery throughout the boundary layer.

34

CHAPTER 3
EFFECTS OF FORM-INDUCED VELOCITY IN TURBULENT HALF CHANNEL
FLOWS

3.1

Introduction

This chapter of the thesis will be also referenced in the text as [91]. The existing studies
summarized in Chapter 1 showed that the form-induced field played a role in affecting the turbulence
response in the presence of roughness. Chapter 2 characterized the overall responses of the turbulent
flows to the acceleration, focusing on the mean flow and single-point turbulence statistics.
In
𝑢𝑖) are found to be much more important for the dynamics of
addition, the form-induced fields ((cid:101)
a transient accelerating channel, compared to fully-developed channel flows or boundary layers
as studied by [38], [40] and [37]. The resultant form-induced production becomes the dominant
source of 𝑣′ at the early stage of the transient, which played a major role in preventing reverse
transition. Furthermore, the apparently synchronous growths of both the TKE pressure-strain term
and the form-induced production at the early stage suggested a possible role of the form-induced
quantities in affecting the TKE redistribution through the pressure-strain term.

To fully understand how the roughness geometry determined the particular turbulence response
𝑢𝑖 affected
in a non-equilibrium transient flow, it is necessary to examine (1) how roughness and (cid:101)
turbulence structure (besides single-point turbulent statistics studied in Chapter 2) and (2) vari-
ous mechanisms through which the form-induced quantities modified Reynolds stress transport.
Specific questions of this chapter included the following.

1. How do the length scales of turbulent and form-induced fields vary in time? How are they

affected by roughness topography?

2. How do the near-wall strain rate of the double-averaged velocities and that of the form-induced
velocities vary in time? are they strong enough to linearly modify the shape of the turbulent
eddies without affecting production (i.e. rapid distortion), similar to the quasi-laminarization
process on a smooth wall?

3. How much do the turbulent and form-induced velocities affected the turbulent pressure
fluctuations and, consequently, the TKE redistribution process as an important mechanism
in the Reynolds stress transport?

The DNS data obtained in Chapter 2 are used for analyses herein. The organization of this chapter
is as follows. The problem formulation is summarized in 3.2, an overview of turbulent statistics
from 2 are shown in 3.3. Temporal variations of the length and time scales of turbulence and the
form-induced velocities are studied in Sections 3.3.1 and 3.4.1, respectively. Various sources of
the pressure fluctuations are compared in Section 3.4.2. Conclusions are provided in Section 4.4.

35

3.2 Problem formulation

Case 𝑅𝑒𝑏0 𝑅𝑒𝜏0
244
SM 4,000
320
4,000
S
294
4,000
T

𝑘 +

𝑠∞,0
0
21
10

(Δ𝑥+, Δ𝑧+)0
(4.0, 2.0)
(5.0, 2.5)
(3.4, 3.4)

𝑅𝑒𝑏1
12,000
12,000
12,000

𝑅𝑒𝜏1
626
1023
858

𝑘 +

𝑠∞,1
0
76
23

(Δ𝑥+, Δ𝑧+)1
(9.8, 4.8)
(15.3, 7.6)
(10.1, 10.1)

Table 3.1 Simulation parameters at initial (‘0’) and new (‘1’) steady states. ‘+’ indicated normal-
ization in wall units (i.e. 𝑢𝜏 and viscous length scale 𝛿𝜈 = 𝜈/𝑢𝜏). Δ𝑥𝑖 is the cell size in 𝑥𝑖 direction.
𝑅𝑒𝑏 = 𝑢𝑏𝛿/𝜈 and 𝑅𝑒𝜏 = 𝑢𝜏𝛿/𝜈. 𝛿 is channel half height.

A summary of the methodologies and simulation setups as reported in [2] is provided here. The

readers are referred to that publication for further details.

The incompressible flow of a Newtonian fluid is simulated by solving the equations of conserva-
tion of mass and momentum. 𝑥, 𝑦 and 𝑧 are, respectively, the streamwise, wall-normal and spanwise
directions of the half-channel flow, and 𝑢, 𝑣 and 𝑤 are the velocity components in those directions.
𝑝 is the pressure, 𝜌 is the density and 𝜈 is the kinematic viscosity. Periodic boundary conditions
are applied in 𝑥 and 𝑧 domain boundaries, and symmetric and no-slip conditions are applied to the
top and bottom boundaries, respectively. In the rough cases, an immersed boundary method is used
to impose the no-slip and no-penetration boundary condition on the rough walls. It is based on the
volume-of-fluid approach; its detailed implementation in the in-house fluid solver is described by
[38] and [79]. The governing equations are solved on a staggered grid using second-order central
differences for all terms, second-order accurate Adams-Bashforth semi-implicit time advancement,
and Message-Passing Interface parallelization.

To conduct the present transient channel simulations, first a separate simulation is carried out for
each case at the initial steady state with a bulk Reynolds number 𝑅𝑒𝑏0, to generate statistics at this
state and be used as initial conditions for the transient simulations. The transient simulations are
then carried out, based on independent initial conditions, for around 20 times for each case to collect
data for ensemble averaging at each time 𝑡. During a transient simulation, the variation of 𝑅𝑒𝑏 (𝑡)
is imposed, equivalent to a rapid three-fold linear increase of the bulk velocity (𝑢𝑏) that started at
time 𝑡 = 0 and lasted for a duration of 𝑡∗ = 𝑡𝑢𝑏0/𝛿 = 0.08 (see Figure 3.2a). 𝑢𝑏 remained constant
thereafter. This simulation setup is similar to that of [56], though with quantitative differences in
the Reynolds number and the amount and rate of the velocity ramp-up.

Simulation parameters for all cases are listed in Table 3.1. The initial and final steady states are
denoted using subscripts “0” and “1”, respectively. “SM”, “S”, and “T” represent the smooth-wall,
sandgrain, and a multiscale turbine-blade-roughness cases, respectively. The same time-variation
of 𝑅𝑒𝑏 is imposed in all cases. Both rough cases started in transitionally rough regime at the initial
steady state and became fully rough at the new steady state, as shown by the 𝑘 +
𝑠∞ values. Here,
𝑘 𝑠∞ is the equivalent Nikuradse sandgrain height in the fully rough regime of a rough surface; the

36

critical values of 𝑘 +
regime are reported by [64].

𝑠∞ for the present roughnesses corresponding to the lower limit of the fully rough

Figure 3.1 (a) Geometries of sandgrain (case S) and turbine-blade-roughness (case T) surfaces
colored by height, with zoomed-in view. (b) Power spectral density of height fluctuations, with
wavenumber 𝜅 in 𝑥 (black) and 𝑧 (gray). (c) Sketch of a rough surface showing various length
definitions.
On a rough wall, the wall shear stress 𝜏𝑤 results from both the viscous and pressure drag due to
roughness. It is calculated as

𝜏𝑤 (𝑡) =

𝜌
𝐿𝑥 𝐿𝑧

∫

V

𝑓

1(𝑥, 𝑦, 𝑧, 𝑡) d𝑥 d𝑦 d𝑧 + 𝜇

𝜕⟨𝑢⟩𝑠
𝜕𝑦

(cid:12)
(cid:12)
(cid:12)
(cid:12)𝑦=0

,

(3.1)

37

where 𝑓1(𝑥, 𝑦, 𝑧, 𝑡) is the streamwise component of the immersed boundary method body force,
V is the volume of the simulation domain, and 𝐿𝑥𝑖 is the domain length in 𝑥𝑖 direction. For the
two rough surfaces considered herein, the second term on the right-hand-side of Equation (3.1) is
negligible. The friction velocity is then obtained as 𝑢𝜏 (𝑡) = √︁𝜏𝑤/𝜌. The corresponding friction
Reynolds numbers are tabulated in Table 3.1; the sandgrain roughness (case S) produced higher
wall friction than the multiscale roughness (case T). Grid sizes of the uniform mesh in 𝑥 and 𝑧 are
shown in wall units in Table 3.1. In 𝑦, the grid is refined near the wall. At the final steady state (i.e.
the more critical one for spatial resolution between the two steady states), Δ𝑦+
< 0.3 for SM and
between 0.6 and 1.7 in the layer below roughness crest for the rough cases.

min

Surface 𝑅𝑎/𝛿

𝑘𝑐/𝛿

𝑘𝑟𝑚𝑠/𝑅𝑎

𝑠𝑘

𝑘𝑢

S
T

0.014
0.014

0.09
0.13

1.05
1.17

0.48
2.97
0.20 3.49

𝑦𝑅/𝑘𝑐

ES𝑧

ES𝑥
0.43 0.44 1.10 ± 0.13
1.05 ± 0.10
0.10

0.08

Table 3.2 Characteristics of the rough surfaces.

Two different surfaces roughness, one synthetic sandgrain roughness ‘S’ and one replicated from a
surface scan on a hydraulic turbine blade ‘T’, are used in the simulations. Figures 3.1(a) and (b)
compare the roughness geometries and the power spectral density of their height fluctuations. Case
S showed a prominent spectral peak at the separation length between neighboring grains in 𝑥 and
𝑧, denoted by 𝜆𝑠, while the T roughness resembled a fractal roughness with a spectral power decay
of −2, without a dominant peak at a particular length scale. For further discussion regarding the
height power spectral densities, see [2].

Parameters of the two rough surfaces are compared in Table 3.2. Both surfaces shared the same
first-order (𝑅𝑎) and similar second-order moments (𝑘𝑟𝑚𝑠) of height statistics, but quite different
values of maximum peak-to-trough height (𝑘𝑐, also called crest height), skewness (𝑠𝑘 ) and effective
slope (ES) in 𝑥 or 𝑧. Both surfaces are positively skewed (i.e. peaky), with kurtosis values of
around 3. Compared to S, the T roughness displayed sparser distribution of peaks, leading to a
lower skewness value and milder effective slopes. The thickness of the RSL (i.e. 𝑦𝑅, defined as the
𝑢2⟩1/2(𝑦)/⟨𝑢⟩(𝑦) > 0.1) is around 1.1 times of 𝑘𝑐 in both cases. 𝑦𝑅 varied mildly
layer in which ⟨(cid:101)
(within around 10%) throughout the transient process. The domain lengths in 𝑥 and 𝑧 are 12𝛿 and
6𝛿, respectively, for cases SM and S, and both 12𝛿 in case T to accommodate its large in-plane
roughness length scales in both 𝑥 and 𝑧. For case S, the spatial resolution of the grain geometry
is 8 and 16 points per grain in 𝑥 and 𝑧, respectively. For case T, there is no dominant roughness
length scale. The Taylor micro-scale (𝜆𝑇 ) [64; 4] is thus used to quantify the length of an equivalent
roughness element. 𝜆𝑇 is resolved by around 10 points in both 𝑥 and 𝑧.

38

3.3 Summary of turbulence statistics

Some observations on flow statistics are summarized here to provide a context for the structural

and other analyses reported in Section 3.4. For comprehensive discussions of the flows, see [2].

(a) bulk velocity, (b) friction coefficient,
Figure 3.2 Flow statistics of the transient channels:
and wall-normal maximum values of streamwise and vertical normal Reynolds stresses (c,d) and
form-induced stresses (e,f). (b-f) are reproduced from [2] with permission.

39

3.3.1 Structural characteristics of turbulent velocity

Figure 3.3 Contours of 𝑅11 at levels 0.3(0.2)0.9 in Cases SM (a), S (b) and T(c), at four time
instances: 𝑡∗ = 0.0, 0.5, 5, and 22.
Lines fitted based on outermost points (an example is shown
with ◦ in (a)) at contour levels 0.5(0.1)0.9, to show inclination angles (marked in red). For clarity,
a shift of 4.5 and 1.5 units are used for smooth and rough cases, respectively.
Figure 3.2(a) showed the acceleration of 𝑢𝑏 that is imposed in all cases. The friction coefficient
𝐶 𝑓 = 2𝜏𝑤 (𝑡)/(𝜌𝑢2
𝑏0) is compared in Figure 3.2(b). 𝐶 𝑓 underwent a sudden increase for all
cases immediately following the 𝑢𝑏 ramp-up, and decreased rapidly as the flow starts to recover.
Figures 3.2(c-f) highlighted the differences in Reynolds and form-induced stress developments in
the three cases.

On the smooth wall, turbulence undergoes reverse transition toward a quasi-laminar state,
characterized by a dip of 𝐶 𝑓 below its long-time value and high Reynolds stress anisotropy (with
TKE mostly resting in 𝑢′), caused by a delayed response of the TKE redistribution to the acceleration.
On the rough walls, however, a clear dip of 𝐶 𝑓 is absent, and Reynolds stress anisotropy stayed
almost unchanged. As opposed to the monotonic increase of 𝑢′
𝑖 fluctuations to a higher-Reynolds-
𝑢𝑖 fluctuations displayed rapid augmentation following the 𝑢𝑏 ramp-up and then
number state, the (cid:101)
𝑢𝑖 intensity is shown to be a result of the
decreased toward the long-time values. This variation of (cid:101)
variation of ⟨𝑢⟩ at the edge of the RSL; the ratio between the root-mean-square (RMS) of (cid:101)
𝑢𝑖 and
⟨𝑢⟩(𝑦𝑅) is shown to be almost constant for each rough case.

3.4 Results

First, the two-point auto-correlation of the turbulent velocity 𝑢′

𝑖 centered at an elevation 𝑦𝑟𝑒 𝑓 is

calculated as

𝑅11(𝑟𝑥, 𝑟𝑦, 𝑡) =

⟨𝑢′(𝑥, 𝑦𝑟𝑒 𝑓 , 𝑧, 𝑡)𝑢′(𝑥 + 𝑟𝑥, 𝑦𝑟𝑒 𝑓 + 𝑟𝑦, 𝑧, 𝑡)⟩

⟨𝑢′2⟩(𝑦𝑟𝑒 𝑓 , 𝑡)

40

,

(3.2)

where 𝑟𝑥 and 𝑟𝑦 are separations in streamwise and wall-normal directions, respectively. The intrinsic
averaging operator in the nominator on the right-hand-side of Equation (3.2) indicated that only
contributions from products with both points ((cid:174)𝑥 and (cid:174)𝑥 + (cid:174)𝑟) in the fluid domain are accounted for in
calculating the correlation. Contours of 𝑅11 in the (𝑥, 𝑦) plane are shown in Figures 3.3 for all cases
at 𝑡∗ = 0.0, 0.5, 5, and 22. In the smooth case, the correlation is centered at 𝑦+
ref,0 = 𝑢𝜏0𝑦ref/𝜈 = 10,
which corresponded to the peak elevation of the TKE production at 𝑡∗ = 0, while in the rough
cases 𝑦ref/𝑘𝑐 = 1, which is near the peak elevations of normal Reynolds stresses [2]. Following
[18], the inclination angle 𝜃 is calculated using the best fitted line by the linear least square method,
based on the points farthest away from the self-correlation point at (𝑟𝑥, 𝑟𝑦) = (0, 0), both upstream
and downstream on the contour lines with levels from 0.5 to 0.9. These outermost points used for
line fitting are shown in the 𝑅11 contours at 𝑡∗ = 0.0 in Figure 3.3(a), as an example. The values
of 𝜃 are listed in the figure. To compare quantitatively the size and inclination of the two-point
correlations, Figure 3.4 showed temporal variations of the streamwise integral length (𝐿11) and 𝜃
in all cases. Here, 𝐿11 is calculated as the integral of 𝑅11 along 𝑟𝑥 from 𝑟𝑥 = 0 to the 𝑟𝑥 value at
which 𝑅11 first decreased below 0.3. A number of 𝑦ref elevations (not shown) are used to evaluate
the dependencies of both 𝐿11 and 𝜃 on 𝑦ref; no significant differences in the comparison among
cases and times are observed.

For case SM, in the reverse transition stage (shown at 𝑡∗ = 5), 𝐿11 increased significantly,
accompanied by a much reduced 𝜃, before both quantities recover toward the long-time values
at 𝑡∗ ≈ 13. The elongation of coherent 𝑢′ motions (i.e.
low speed streaks) is also observed
for equilibrium accelerating turbulence such as the sink flow [79] and turbulent boundary layer
subjected to FPG [47]. In both rough cases, on the other hand, 𝐿11 reduced rapidly, immediately
after the step acceleration, reaching a minimum at 𝑡∗ ≈ 0.4 for case S and at a later time, 𝑡∗ ≈ 2.0, for
case T. After then, 𝐿11 slowly increased toward the long-time values. During this time, an overall
monotonic increase of 𝜃 is observed on the rough walls. The change in roughness texture only
affected quantitative aspects of the large 𝑢′ motions, but not their trends: the multiscale roughness
yielded slower variations of both 𝐿11 and 𝜃 following the acceleration at 𝑡∗ = 0.

Now focus on the initial (𝑡∗ = 0) and late-time states (𝑡∗ = 22). In the presence of sandgrain
roughness, the 𝑥 extent is shorter than that in the smooth case at both times. This is because
the sandgrain elements are much smaller than 𝛿; it imparts its length scale to the near-wall flow
by breaking up long streaky motions near the wall. On the multiscale roughness, however, the
initial 𝑥 extent is closer to that in the sandgrain case, while at late times it approaches the smooth-
wall values. This suggested that the streaks are broken up only by the sparsely distributed peaky
elements, which are separated by large roughness scales. As the Reynolds number increased after
the acceleration, the streaks appeared less affected, as the streamwise distances separating the peaky
elements approached the order of magnitude of smooth-wall streak lengths, 𝑜(103)𝛿𝜈.

41

Figure 3.4 Temporal variation of streamwise integral length (a,c) and inclination angle (b,d) based
on 𝑅11 in SM (◦) , S (△), T (□) cases. The correlations are centered at 𝑦+
ref,0 = 10 for SM and
𝑦ref/𝑘𝑐 = 1 for S and T. (a,c) show linear scalings of 𝑡∗ and (c,d) show logarithmic scalings.
To characterize turbulent motions of all sizes, the premultiplied one-dimensional (1D) power spectra
(𝜅1𝐸 (𝜅1)) of 𝑢′ and 𝑣′ with streamwise wavenumber 𝜅1 = 2𝜋/𝜆1 (where 𝜆1 is the streamwise
wavelength) are shown in Figures 3.5. First, a remark is given on the DNS spatial resolution,
which is deemed sufficient if the wavenumber of the largest energy-containing motions and the
highest physically “meaningful” wavenumber are two or more decades apart [92]. The highest
meaningful wavenumber is the 𝜅1 value at which the effective resolution, quantified by an effective
1,effΔ𝑥/𝜋 [93], reached its maxima. For the present second-order central-
wavenumber 𝜅∗
difference scheme, the effective wavenumber of the first derivative 𝜅∗
𝜋)/𝜋 [93] obtained
a maximum at 𝜅1 = 𝜋/(2Δ𝑥). Since the effective wavenumber droped for 𝜅1 higher than this
value, the simulated fluid motions at scales smaller than 𝜋/(2Δ𝑥) are not physically meaningful.
Figures 3.5(a-c) showed that 𝜋/(2Δ𝑥) is around two decades higher than the 𝛿 scale in all cases,
satisfying the criterion of [92] for DNS spatial resolution.

= sin(𝜅∗
1

= 𝜅

1,eff

1,eff

42

Figure 3.5 1D premultiplied power spectra of 𝑢′ (a-c) and 𝑣′ (d-f) at 𝑡∗ = 0 (
) and a series of
, 𝑡∗ values from 0.1 to 0.9 with step 0.1, from 1 to 9 with step 1, and from 10 to 22 with
times (
step 3, lighter colors at later times), in smooth-wall case (a,d) at 𝑦+ = 10 and cases S (b,e) and T
(c,f) at 𝑦 = 𝑘𝑐. In (d-f), 𝐸22 at 𝑡∗ = 0 is scaled by a factor of 15 for clarity. In (a-c),
(vertical):
highest physically meaningful wavenumber. In (a,d), the profile corresponding to 𝑡∗ = 5 (discussed
in text) is highlighted in red.

Figures 3.5(a,d) showed that, for the smooth-wall case, in the reverse-transition phase (e.g. at
𝑡∗ = 5 marked in red) the energy of 𝑢′ motions is concentrated at large scales of order 𝛿, with
only slight changes in the peak location of the 𝑣′ spectrum. This is consistent with the streamwise
elongation of turbulent coherent motions observed in a quasi-laminar flow [44; 58]. On the other
hand, case S (Figures 3.5 (b,e)) does not show an increase of the low-wavenumber 𝑢′ content.
Instead, the peaks of both 𝑢′ and 𝑣′ power spectra shift monotonically toward the dominant length
scale of the roughness geometry (i.e. distance between centers of neighboring grains, 𝜆𝑠, shown in
Figure 3.1(b) as 𝜅𝜆𝑠 ≈ 50). One can infer that the 𝑣′ associated with the new turbulence generated
around the sandgrain roughness had a characteristic streamwise length of the grain separation, and
that the 𝑢′ motions are of lengths two to three times as long, which probably reflects channeling
phenomena between grains. In comparison to the sandgrains, the T roughness (Figures 3.5 (c,f))—
characterized by a wide spectrum of length scales up to around 6𝛿—led to a large portion of 𝑢′ being
distributed in the scales between 𝑜(𝛿) and 𝑜(0.1𝛿) during the transient process. One similarity can
be seen between case T and the smooth case: both cases showed an amplification of the large-scale
𝑢′ content at early times (though to a lesser extent in case T), Suggesting that, similar to the smooth
case, streamwise-elongated eddies are also present around T roughness during the response to

43

acceleration.

3.4.1 Time scales of mean-flow strain rates

The previous section showed that time-variation of the structure of near-wall turbulence de-
pended on the in-plane scales of roughness, with flow over the multiscale T roughness displaying
some similarities with that on a smooth wall in terms of velocity coherence length and spectra.
Past studies on smooth-wall turbulence attributed the observed streamwise elongation of inclined
near-wall turbulent eddies to the stretching due to mean shear, which is quantified using a time-scale
ratio (also called “shear rate parameter”), between the turbulence time scale and the time scale of
the mean shear. For a transient channel, the time-scale ratio is

𝑆∗ =

𝑆(𝑦, 𝑡) 𝑞2(𝑦, 𝑡)
2𝜖 (𝑦, 𝑡)

,

where

𝑆(𝑦, 𝑡) = 2(𝑆𝑖 𝑗 𝑆𝑖 𝑗 )1/2,
1
2

(cid:18) 𝜕⟨𝑢𝑖⟩
𝜕𝑥 𝑗

𝑆𝑖 𝑗 (𝑦, 𝑡) =

+

𝜕⟨𝑢 𝑗 ⟩
𝜕𝑥𝑖

(3.3)

(3.4)

(3.5)

(cid:19)

,

𝑞2 is twice the TKE, and 𝜖 is the TKE dissipation rate. 𝑞2/𝜖 represented roughly the time required
for the energy of an energy-containing eddy to fully break down to motions in the dissipative
subrange through the energy cascade process. [94] showed that, for homogeneous turbulence under
a strong uniform mean shear, the turbulent eddies are distorted by the mean flow so rapidly that
they are simply stretched, other than passing down the effect into smaller-scale eddies through the
cascade.
In other words, the turbulence time scale is much larger than that of the mean shear,
[88] showed that high 𝑆∗ values between 30 and 40 imposed in initially
reflected by 𝑆∗ ≫ 1.
homogeneous turbulence (without presence of a wall) can produce turbulent velocity and vorticity
anisotropies similar to those observed near the wall in a channel flow. [58] demonstrated that in
the region of quasi-laminar flow in a strongly accelerated smooth-wall boundary layer, near-wall 𝑆∗
values reached up to around 60, compared to a value around only 10 on a rough wall where reverse
transition is prevented. In the logarithmic region of the flow, 𝑆∗ is around 4 to 5, a result of locally
equilibrium turbulence.

The 𝑦 profiles of 𝑆∗ in all cases are compared in Figure 3.6. Here, the only non-zero component
of the 𝑆𝑖 𝑗 tensor is the 12 component, 𝜕⟨𝑢⟩/𝜕𝑦. In the smooth case, high magnitudes of 𝑆∗ near 70
are achieved at 𝑡∗ ≈ 5 near the wall, expected for the reverse transition process. In the presence of
roughness, the wall-normal peak of 𝑆∗(𝑦) occured near the roughness crest, where the 𝑦 gradient
of ⟨𝑢⟩ is the strongest. The 𝑆∗ peak values stayed between 5 and 15 throughout the transient, much
lower than those in the smooth case, suggesting that the rough-wall DA flow is not strong enough to
linearly reshape turbulent eddies. Far from the wall (𝑦/𝛿 > 0.1, or 𝑦+ > 100 herein), 𝑆∗ values are

44

at around 5 at all times, regardless of wall condition. In case T near 𝑦 = 0, one observed significant
fluctuations of 𝑆∗ values. This is probably due to fluctuations in the 𝑦-gradient of the roughness
geometry function (or fluid area fraction), Φ(𝑦), in case T, through its involvement in the intrinsic
averaging process, ⟨𝜃⟩ = 1/(Φ𝐴𝑜) ∫
𝜃d𝐴 (where 𝐴𝑜 is the total planar area). See details of Φ in
the Appendix. Using superficial averaging, however, would mask these fluctuations as ⟨·⟩𝑠 = Φ⟨·⟩
and Φ ≈ 0 near the wall.

𝐴 𝑓

Figure 3.6 Variation of 𝑆∗ for (a) SM , (b) S , and (c) T cases in time. In (b) and (c), vertical dotted
lines indicate roughness crest height.

The discussions above, as well as those in [58], focus on the ratio between global (space-averaged)
time scales of turbulence and double-averaged mean strain, instead of the local values. However,
𝑢𝑖 inside the RSL may yield local time scale ratios
in rough-wall cases the strongly heterogeneous (cid:101)
much higher than the global value. In particular, Figure 3.2 showed that, within a short time after
the impulse acceleration (at 𝑡∗ ≈ 0.3), the peak of ⟨(cid:101)
𝑢2⟩ reached one order-of-magnitude higher than
that of ⟨𝑢′2⟩. This is expected to yield strong local gradients of form-induced velocity inside the
RSL, which is quantified next.

We defined a time-scale-ratio tensor, 𝑆∗ 𝑓

𝑢𝑖.
𝑖 𝑗 , to quantify the various strain rate components of (cid:101)
It is similar to the definition of 𝑆∗, but based on the form-induced velocity instead of the DA one,

𝑆∗ 𝑓
𝑖 𝑗 (𝑦, 𝑡) =

(cid:42)

(cid:101)𝑆𝑖 𝑗 𝑞2(𝑥, 𝑦, 𝑧, 𝑡)
𝜖 (𝑥, 𝑦, 𝑧, 𝑡)

(cid:43)

,

𝑢𝑖,
where (cid:101)𝑆𝑖 𝑗 is the rate-of-strain tensor of (cid:101)

(cid:101)𝑆𝑖 𝑗 (𝑥, 𝑦, 𝑧, 𝑡) =

(cid:18) 𝜕
𝑢𝑖
(cid:101)
𝜕𝑥 𝑗

1
2

+

(cid:19)

,

𝑢 𝑗
𝜕
(cid:101)
𝜕𝑥𝑖

(3.6)

(3.7)

and the superscript “ 𝑓 ” indicated that the time-scale ratio is calculated using the form-induced
velocity. Note that, different from the 𝑆∗ definition in Equation (3.3), local values of three-
dimensional ensemble-averaged 𝑞2 and 𝜖 are used here to evaluate the local time scale ratio. 𝑆∗ 𝑓
𝑖 𝑗
components are negligible for the smooth-wall flow.

45

𝑢𝑖 strain rates for (a-b) S and (c-d) T cases at three times: steady
Figure 3.7 Time-scale ratios of (cid:101)
state 𝑡∗ = 0 and long-time state 𝑡∗ = 25 (a,c), as well as a transient state 𝑡∗ = 0.5 (b,d). Profiles
corresponding to 𝑡∗ = 25 are offset upward by 10 units for clarify. Vertical dotted lines indicate
roughness crest height.
Wall-normal profiles of all independent components of the 𝑆∗ 𝑓
tensor in both rough cases are
𝑖 𝑗
compared in Figures 3.7. The r.m.s. values of (cid:101)𝑆𝑖 𝑗 are used in the Equation (3.6) to calculate 𝑆∗ 𝑓
𝑖 𝑗 .
The regions of negative (cid:101)𝑆𝑖 𝑗 (indicating contraction of time-mean coherent motions) are expected to
be limited to small regions near the small wavelength rough peaks, indicating a lower contribution
to the overall r.m.s values. Among all, the 12 component is the dominant component for both cases,
at all times. Within the 12 component, we found that 𝜕
𝑣/𝜕𝑥,
(cid:101)
𝑢 is the dominant form-induced velocity component. Overall, the plane-averaged form-induced
as (cid:101)
strain rate is similar or weaker compared to the strain rate of ⟨𝑢𝑖⟩ during the transient process, with
the exception of the 12 component in case T at early stage of the transient (shown at 𝑡∗ = 0.5),
which increased to around 18 below the crest elevation. This is roughly 2.5 times the shear rate
of ⟨𝑢𝑖⟩ at this time (Figure 3.6c). The higher plane-averaged form-induced time-scale ratio on T
𝑢 on T
roughness than on S can be attributed to a stronger 𝜕

𝑢/𝜕𝑦 is significantly larger than 𝜕
(cid:101)

𝑢/𝜕𝑦, resulting from more intense (cid:101)
(cid:101)

46

roughness than on S (Figure 3.2e), together with almost identical RSL thicknesses (Table 3.2). In
addition, the comparison between the time-scale ratios at the initial and long-time states (𝑡∗ = 25)
in Fig. 3.7(a,c) showed that they are similar, with the main difference being slightly higher values
of 𝑆∗ 𝑓
𝑖 𝑗 of all components in the outer region at 𝑡∗ = 25 than those at 𝑡∗ = 0. This is probably due to
more intense high- and low-momentum pathways (as discussed by, for example, [23] and [95]) in
the outer region under the higher Reynolds number. These pathways are discussed, for the present
roughnesses and at a Reynolds number similar to the present long-time state, by [87]. On the T
𝑢𝑖
strain rates than on the smooth wall in the spatial-average
roughness, despite the much lower (cid:101)

12

Figure 3.8 Contours of (cid:101)𝑆12𝑞2/𝜖 at 𝑡∗ = 0.5, at respective peak elevations of 𝑆∗ 𝑓
(𝑦) profiles for
cases S (a, at 𝑦/𝑘𝑐 = 1) and T (b, at 𝑦/𝑘𝑐 = 0.75), zoomed in to the regions shown in the full height
maps in (c).
sense as seen by comparing Figures 3.6(a) and 3.7(d), locally the time scale ratio can achieve a
level similar to that on the smooth-wall, due to heterogeneous nature of this roughness. Figure 3.8
compares wall-parallel contours of local contribution to 𝑆∗ 𝑓
12 on both rough walls at 𝑡∗ = 0.5. The
𝑦 elevations of the contours are selected as the respective peak elevations of the average time-scale
ratios in the two cases according to Figure 3.7. Only a portion of the domain is shown for clarity.
In case S, local time-scale ratio rarely exceeds 20, and thus the augmented form-induced straining
during the transient process acted mainly to promote turbulence production. However, in case T
high time-scale ratios are achieved locally in the regions between roughness peaks, with a large
area of flow reaching a ratio of 60, similar to the values on a smooth wall.

To quantitatively compare the time-scale ratios in the two rough cases, the cumulative distribu-
tion functions (CDF) of (cid:101)𝑆12𝑞2/𝜖 in the whole (𝑥, 𝑧) domain at the same 𝑦 elevations are compared
in Figure 3.9. For case T, around 20% of the fluid domain is subjected to local time-scale ratios
higher than 30 (i.e. level seen in smooth-wall fully developed channels near the wall) and 5%

47

Figure 3.9 Cumulative distribution functions of local time-scale ratios of (cid:101)𝑆12 shown in
Figures 3.8(a,b).
corresponded to ratios higher than 60. Therefore, in the presence of large roughness length scales
𝑢/𝜕𝑦 became locally
(with case T as an example), it is possible that strong positive values of 𝜕
(cid:101)
important in early stage of the transient as it led to streamwise elongation of inclined turbulent
eddies (indicated by 𝑅11 contours in figure 3.3c) near the wall. This is consistent with the changes
in integral length scale shown in Figure 3.4 and in the 𝑢′ spectra shown in Figure 3.5. Also, addition
of 𝑆∗ to 𝑆∗ 𝑓
𝑖 𝑗 (𝑦, 𝑡) is expected to further increase the magnitude and occurrence of the positive values
of 𝜕

𝑢/𝜕𝑦, further emphasizing the above findings.
(cid:101)

3.4.2 Sources of turbulent pressure fluctuations

Beside the more widely studied form-induced production, we suggested a second mechanism
through which the form-induced velocities modified the Reynolds stress balance: indirect augmen-
tation of TKE redistribution (i.e. Π in Equation 2.8) by intensifying turbulent pressure fluctuations,
𝑝′. One supporting evidence is that both form-induced velocity intensities and the TKE redistri-
bution to 𝑣′ through pressure-strain term increased rapidly at early times in the rough cases, as
compared to a slow increase of TKE redistribution on a smooth wall [2].

The following is a selective summary of current understanding of how roughness affects 𝑝′.
Based on DNS of rough-wall channel flows, [96] observed that roughness augmented pressure
fluctuations throughout the boundary layer, and that the RSL thickness defined based on 𝑝′ RMS
is thicker that those based on 𝑢′
𝑖 or vorticity RMS. [97] characterized the wall-pressure spectrum
[98] showed with
and analyzed its scaling based on high-Reynolds-number experimental data.
large-eddy simulations that roughness geometry and distribution modify wall-pressure fluctuations.
Specially, intense pressure fluctuations occur on the front surface of roughness, instead of the back
surfaces, suggesting that instantaneous impinging flow structures are more important than vortex
shedding to wall-pressure RMS and generation of self noise. [35] examined the two 𝑝′ source terms
that are identified as the dominant ones in fully developed smooth-wall flows: the 12 component

48

) and 22 (

Figure 3.10 (a) RMS of the left-hand-side (red) and right-hand-side (black) of 𝑝′ Poisson’s equation
(Equation 3.8) at 𝑡∗ = 0.5 (
) for case S. (b,c) Steady-state RMS of source terms
in cases S (b) and T (c), evaluated at 𝑡∗ = 22. All quantities are normalized using 𝑢𝑏0 and 𝛿.
of the mean shear term (“MS12” therein) and the 23 component of the turbulent-turbulent term,
and concluded that the MS12 term is a more important source between the two on the rough wall.
They also explained the links between these terms and local structure of the mean and turbulent
velocities.

Due to non-zero (cid:101)
Poisson’s equation of 𝑝′,

𝑢𝑖 in rough-wall flows, an additional term (represented by F ) appears in the

where the source terms are

1
𝜌

∇2 𝑝′ = M (𝑥𝑖, 𝑡) + T (𝑥𝑖, 𝑡) + F (𝑥𝑖, 𝑡),

,

𝜕𝑢′
𝑗
𝜕𝑥𝑖
(cid:16)

𝑖𝑢′
𝑢′

𝑗 − 𝑢′

𝑖𝑢′
𝑗

(cid:17)

,

and

M = −2

T = −

F = −2

𝜕⟨𝑢𝑖⟩
𝜕𝑥 𝑗
𝜕2
𝜕𝑥𝑖𝜕𝑥 𝑗
𝜕
𝑢𝑖
(cid:101)
𝜕𝑥 𝑗

𝜕𝑢′
𝑗
𝜕𝑥𝑖

.

(3.8)

(3.9)

(3.10)

(3.11)

Here M is the source due to gradients of double-averaged velocity (i.e. the ‘mean shear term’), T
𝑢𝑖. It is unclear how
is the turbulent-turbulent term, and F is the source due to the gradients of (cid:101)
these terms compare, since, to our knowledge, quantitative comparison of all terms had not been
reported in the context of rough-wall flows.

In the following, we compare the magnitude of each source term as a function of time, for
all cases. First, to estimate the accuracy of the source term calculation in the context of the
immersed boundary method, the RMS of the left-hand-side and right-hand-side of Equation (3.8)
are compared at various time instances in Figure 3.10(a), taking case S as an example. The RMS
is calculated based on plane averaging using one snapshot at each time (i.e., without ensemble
averaging) as (·)rms = ⟨(·)2⟩1/2

, where (·) represented a local instantaneous flow quantity. The

𝑠

49

Figure 3.11 Temporal variations of RMS of 𝑝′ source terms at their respective wall-normal peak
elevations for the smooth case (a) and at 𝑦 = 0.7𝑘𝑐 for cases S (b) and T (c).
difference between the two sides of the equation is shown to be up to 5-8% of the local RMS values.
This difference is attributed to discretization error of the finite difference scheme used to calculate
the terms in Equation (3.8) and to the exclusion of grid cells that are cut by the fluid-solid interfaces
when calculating these terms.

1/2

The 𝑦-profiles of the RMS of all terms near the new steady state (𝑡∗ = 22) are compared in
Figure 3.10(b,c) between the two rough cases. Here both ensemble and plane averaging are used
in calculating the RMS, obtained as (·)rms = ⟨(·)2⟩𝑠
. [35] observed that the M12 term is more
important than the T23 term, in connection with the presence of maximum 𝜕⟨𝑢⟩/𝜕𝑦 inside the
RSL and the 𝜕𝑣′/𝜕𝑥 augmented by roughness. Surprisingly, the results here show that the least
important source is M in both rough cases, when all components of all term are considered in the
calculation. This is probably due to the relatively weak 𝜕⟨𝑢⟩𝑖/𝜕𝑥 𝑗 compared to the form-induced
shear 𝜕
𝑖 components in the RSL. In addition, T is shown
here to be the dominant source term regardless of roughness texture. All terms reach a peak slightly
below the roughness crest. This is because that the peaks of form-induced stresses [40], Reynolds
stresses, and DA mean shear [9; 99; 100] all reach their respective peak values near the roughness
crest. In case T, the peaks of the source terms are located lower (toward the middle of the RSL),
because of the rarer attainment of the maximum crest height by local roughness peaks due to the
presence of large surface scales.

𝑢𝑖/𝜕𝑥 𝑗 , as well as strong gradients of all 𝑢′
(cid:101)

Next, the RMS values of each of the three source terms at a 𝑦 elevation are evaluated for all times
in Figure 3.11. For the smooth-wall case, as the peak elevations of different terms are not the same
and vary with time, the source terms are evaluated at the instantaneous peak location of each term.
In the rough-wall cases, however, all terms are evaluated at 𝑦 = 0.7𝑘𝑐, near the peak elevations of
all source terms, independent of time. In the rough cases, the M term is the least important one
not only at steady states, but also during the transient, as turbulence and form-induced velocities in
the RSL adjust to the accelerated mean velocity rapidly. All terms increase instantly following the

50

Figure 3.12 Contours of F (a-c) and T (d-f) at 𝑦 = 0.7𝑘𝑐 for cases S, at 𝑡∗ = 0.1 (a,d), 0.5 (b,e),
and 10 (c,f). Black contour lines mark fluid-solid interfaces. All quantities are normalized based
on 𝑢𝑏0 and 𝛿. 1/8 of the domain is shown.
step acceleration on rough walls. Specifically, F becomes particularly significant in the early stage,
𝑢𝑖 (Figures 3.2(e,f)). The T term increased synchronously with F ,
due to the intensification of (cid:101)
probably due to the intensification of TKE associated with stronger form-induced production. On
𝑢𝑖, the T term is not augmented immediately
the smooth wall, due to the absence of roughness and (cid:101)
following the bulk acceleration. Instead, the turbulence is “frozen”, signified by the persistence of
source terms at the initial-state levels till 𝑡∗ ≈ 3. The decorrelation between mean flow acceleration
and the developments of M and T terms are also observed in a smooth-wall boundary layer subject
to strong spatial acceleration [44]. These results confirmed that characteristics of the form-induced
𝑢𝑖
velocity as contained in F are important for the 𝑝′ evolution. They also suggested that the (cid:101)
distribution determines how turbulence responses to strong acceleration (with more evidence to be
shown next).

Spatial distributions of F and T terms in both rough cases are compared in Figures 3.12 and
3.13, at three time instances: 𝑡∗ = 0.1, 𝑡∗ = 0.5, and 𝑡∗ = 10.
In both cases, regions with strong
F and those with strong T are well correlated at all three time instances, as shown in Table 3.3
by the single-point cross-correlation between the two sources of values around 0.5 to 0.6. This
suggested that early generation of new turbulence depends primarily on form-induced production.
𝑢𝑖 strain rates form not just an important source of 𝑝′, but also the
In other words, it appears that (cid:101)
origin of the other important source, T .

In addition, comparison between Figures 3.12 and 3.13 highlighted the effect of roughness
texture on the rate of increase of 𝑝′. First, focus on the largest source term, T . At early times this

51

Figure 3.13 Contours of F (a-c) and T (d-f) at 𝑦 = 0.7𝑘𝑐 for cases T, at 𝑡∗ = 0.1 (a,d), 0.5 (b,e),
and 10 (c,f). Black contour lines mark fluid-solid interfaces. All quantities are normalized based
on 𝑢𝑏0 and 𝛿. 1/16 of the domain is shown.
term takes high magnitudes in the vicinity of peaky roughness elements (e.g. see Figures 3.12(e)
and 3.13(e)).
In case S, the peaky elements associated with the tips of the grains are densely
distributed, while in case T the roughness peaks are sparser, separated by large roughness scales.
As a result, at early times (0 < 𝑡∗ < 0.3) the rate of increase of T in case T is slower compared
to that in case S, as shown in Figure 3.11. In addition, Figure 3.11 showed that, in case S, after
F reached its peak in time at 𝑡∗ ≈ 0.3 and started to decrease toward the steady state, the T
term followed the same trend and decreased from its peak (another evidence of the importance of
form-induced production of TKE). In comparison, in case T the T term displayed only a slightly
decrease following the decline of F at 𝑡∗ ≈ 0.3. Examination of Figure 3.13(e,f) explained that
this almost constant T level in case T for 𝑡∗ > 0.3 is due to a combination of both the weakening
of local T values (due to a decrease in F ) and the delayed activation of other regions in producing
intense turbulent fluctuations.

Case

𝑡∗ = 0.1 𝑡∗ = 0.5 𝑡∗ = 10

S
T

0.52
0.53

0.56
0.48

0.56
0.50

Table 3.3 Single-point cross-correlation between F and T contours shown in Figs. 3.12 and 3.13.
The correlation is defined as ⟨|F | · |T |⟩/(cid:2)⟨F 2⟩1/2 · ⟨T 2⟩1/2(cid:3).

To further characterize the F term, individual components −2(𝜕

pared in Figure 3.14. For case T, the 12 term −2(𝜕
terms at all times. This is because that the large roughness scales induced significant (cid:101)

𝑢 in the

𝛽/𝜕𝑥𝛼) are com-
𝑢𝛼/𝜕𝑥𝛽) (𝜕𝑢′
(cid:101)
𝑢/𝜕𝑦) (𝜕𝑣′/𝜕𝑥) is significantly higher than other
(cid:101)

52

Figure 3.14 Temporal variations of the RMS of individual terms in F for cases S (a) and T (b).
Symbols: + sum of all terms, ⃝ 11, △ 12, ▽ 13, ⃝ 21, △ 22, ▽ 23, ⃝ 31, △ 32, ▽ 33.
RSL [40]. In case S, however, the 32 and 23 terms are the largest ones, while other terms are also
𝑣 and
significant. The difference between the two cases highlighted the dynamical importance of (cid:101)
𝑢 in the case of a multiscale roughness.
𝑤 around densely distributed roughness, versus that of (cid:101)
(cid:101)
3.5 Conclusions

This work evaluated various pathways through which the form-induced velocity (cid:101)

𝑢𝑖 affected
turbulent flows in equilibrium and non-equilibrium states, based on DNS data of transient turbulent
half channels in response to a three-fold step acceleration of bulk velocity at time 𝑡∗ = 0. Two
roughness geometries—a densely distributed synthesized sandgrain roughness (case S) and a
multiscale turbine-blade roughness (case T)—are imposed on the channel wall; a baseline smooth-
wall case (case SM) is also compared. The two roughness geometries shared similar first and
second moments of roughness height fluctuations, but differ in in-plane scales.

Analyses of two-point turbulent velocity correlations in the transient flow showed that, opposed
to the smooth-wall cases where energy-containing motions are elongated in 𝑥 and decrease in
inclination angle as signatures of the reverse transition process, on the rough walls the streamwise
coherence decreases rapidly after the acceleration and the inclination angle increases sharply.
Energy spectra, on the other hand, highlight a difference between the rough cases: while in case
S the spectral peaks of 𝑢′ and 𝑣′ motions monotonically shift toward characteristic scales of the
sandgrains, in case T the 𝑢′ motions are initial augmented at large scales; the latter suggests local
streamwise elongation of velocity streaks similar to the smooth-wall case to some extent.

We show that the dependence of turbulent structure on the in-plane scales of roughness described
𝑢𝑖 straining and turbulence.
above may be explained by the local time-scale ratio between the (cid:101)
Although both strain rate of double-averaged velocity ⟨𝑢𝑖⟩ and the spatial-averaged local strain rate
𝑢𝑖 is shown to
of (cid:101)
give local time-scale ratios above 60 in case T at early times, in large areas between sparse peaks of

𝑢𝑖 yielded much smaller time-scale ratios compared to values on a smooth wall, (cid:101)

53

roughness. These values are similar to the level reached when the turbulence is relaminarized on
a smooth wall. This result indicated that, in the presence of in-plane roughness scales of order 𝛿,
𝑢/𝜕𝑦 may be important during the transient response and result in linear stretching of turbulent
𝜕
(cid:101)
eddies.

Besides the additional form-induced production known to be important in transient channel
flows [2], evidence of a second mechanism through which the form-induced velocity indirectly
affected the Reynolds stress balance by augmenting 𝑝′ is presented. We compared all three 𝑝′
source terms in the roughness sublayer and found that, surprisingly, the least important one at all
times is the source due to ⟨𝑢𝑖⟩ shear. Both the source due to (cid:101)
𝑖 are
important, with the spatial distribution of the latter correlated with the former, probably due to the
production of TKE from (cid:101)

𝑢𝑖 gradients and that due to 𝑢′

𝑢𝑖 shear.

This work provided evidences that the form-induced velocity inside the roughness sublayer is
𝑢𝑖 gradients form
important for turbulence statistics and structure in wall-bounded turbulent flows. (cid:101)
an important source of 𝑝′ in an equilibrium flow, while in a non-equilibrium accelerating flow
they generate significant amount of TKE in the vicinity of roughness peaks and reshaped turbulent
eddies in the regions between peaks.

54

CHAPTER 4
DATA-DRIVEN PREDICTIONS OF THE ROUGHNESS SUBLAYER FLOW

4.1

Introduction

Chapters 2 and 3 demonstrated the influence of roughness geometry on the turbulence evolution
in a type of non-equilibrium wall turbulence and the important role played by the form-induced ve-
𝑢𝑖) in shaping the turbulence response, including evolution of turbulent eddies and Reynolds
locity ((cid:101)
stress transport. The next question is to understand how roughness texture affects the form-induced
stress. [40] studied the form-induced velocity distribution on rough surfaces with different length
𝑢 to be dependent on the
scales. They showed the occurrence of both positive and negative values of (cid:101)
wavelengths of the surface. Additionally, the long-wavelength surface components in a multiscale
𝑤, leading to negligible form-induced shear stress
surface led to lower occurrences of intense (cid:101)
unlike on a homogeneous sandgrain surface. Studies on regular arrays of cubes and on gravel-
𝑣 leading to
bed roughness by [101] and [37], respectively, observed large magnitudes of (cid:101)
strong negative form-induced shear stress near re-circulation regions. These studies, among others,
showed the dependencies of various components of the form-induced velocities on the roughness
texture. However, past analyses were typically done for single types of rough surfaces only. The
current work is focused on expanding the understanding to a wide range of rough surfaces.

𝑢 and (cid:101)

𝑣 and (cid:101)

Instead of analyzing the 3D (cid:101)

low-dimension quantities that are essential for the transport of (cid:101)
𝑢𝑖⟩/2) [36; 37; 38], in a periodic half-channel flow reads
𝑢𝑖 (cid:101)
kinetic energy (WKE, i.e. ⟨(cid:101)

𝑢𝑖 (𝑥, 𝑦, 𝑧) field itself, it is proposed to address, in this chapter,
𝑢𝑖. The budget equation for the wake

𝑣⟩
𝑢𝑖(cid:101)
−⟨(cid:101)

𝜕𝑈𝑖
𝜕𝑦

(cid:28)

+

𝑖𝑢′
(cid:103)𝑢′
𝑗

(cid:29)

𝜕
𝑢𝑖
(cid:101)
𝜕𝑥 𝑗

1
2

−

𝑠

𝑣⟩𝑠
𝑢𝑖(cid:101)
𝑢𝑖 (cid:101)
𝜕⟨(cid:101)
𝜕𝑦
(cid:68)

𝜕

(cid:69)
(cid:101)𝑃
𝑣
(cid:101)
𝜕𝑦

+

(cid:42) 𝜕 (cid:101)𝑃
(cid:43)
𝜕𝑥

−𝑈

(cid:68)

𝜕

𝑖𝑣′(cid:69)
𝑢𝑖 (cid:103)𝑢′
(cid:101)
𝜕𝑦

𝑠

−

𝑠

+ 𝜈 (cid:10)

𝑢𝑖∇2
(cid:101)

𝑢𝑖(cid:11)
(cid:101)

𝑠 = 0.

(4.1)

(4.2)

𝑢𝑖 fluctuations are: (i) the production due to the mean shear (i.e.
The two sources of production of (cid:101)
𝜕𝑈/𝜕𝑦, with 𝑈 ≡ ⟨𝑢⟩𝑠), as shown by the first term on the left-hand-side in Equation (4.1), and (ii)
the work of pressure drag of roughness (i.e. −⟨𝜕 (cid:101)𝑃/𝜕𝑥⟩) against the mean velocity 𝑈, as shown
by the fifth term on the left-hand-side. Therefore, understanding how both 𝑈 (𝑦) (and 𝜕𝑈/𝜕𝑦 by
extension) and the pressure drag depend on the roughness geometry is essential in understanding
how the roughness geometry controls WKE development and, in turn, the turbulence response in
non-equilibrium flows.

First, past studies on the modeling of 𝑈 (𝑦) based on roughness geometry are summarized below.
[85] provided analytical derivation of 𝑈 (𝑦) profiles for 2D, steady rough-wall turbulent flows based
on the flow conditions, roughness geometrical function 𝜙(𝑦) (i.e. porosity) and the relative
submergence of roughness. Specifically, roughness with constant 𝜙 produced an exponential

55

velocity profile, while a linear velocity profile is applicable to roughness with monotonically
increasing 𝜙(𝑦). Another study by [100] over varying height, aspect ratio and distribution of
rectangular-prism roughness elements predicted an exponential 𝑈 (𝑦) profile with an attenuation
[75] proposed a composite
factor obtained based on wake interaction within rough surfaces.
𝑈 (𝑦)-profile model based on observations of the 𝑈 (𝑦) of four different rough surfaces: sandgrain
roughness, half ellipsoids, 2D wavy surface and a random surface generated with spatial Fourier
modes. The composite model consists of a zero-velocity layer near the wall, combined with an
upper layer with 𝑈 (𝑦) prescribed based on a polynomial expression, which smoothly blends the
zero-velocity layer and the logarithmic layer. The model coefficients are determined by fitting to
DNS data of the four rough surfaces. The model is successfully integrated into channel-flow RANS
simulations that did not resolve the RSL. However, the composite model has limited application
as it is fitted and tested for a few surfaces only. Recent study on spatially developing flat-plate
boundary layers under acceleration and deceleration [3] showed that 𝑈 (𝑦), when appropriately
normalized by 𝑈𝑅, becomes self-similar inside the RSL regardless of pressure gradients, as far as
the flow is attached and under high 𝑘 +. This is shown in Figure 4.1, where 𝑦𝑅 is the roughness
sublayer thickness and 𝑈𝑅 is the mean velocity at 𝑦 = 𝑦𝑅. This is an important finding as it
indicates that a 𝑈 (𝑦) model developed for equilibrium flows may be directly applicable to a subset
of non-equilibrium flows.

Next, the problem of modeling the pressure drag based on the roughness topography can be
reduced to the modeling of the topographical dependencies of the equivalent sandgrain height (𝑘 𝑠,
as defined in Section 1.2.1) instead. A detailed explanation of the relation between 𝑘 𝑠 and the
logarithmic shift in the mean velocity profile can be found in the Section 1.2.1. This logarithmic
shift measures the higher total drag on a rough wall than on a smooth wall, which is mostly pressure
drag under sufficiently high 𝑘 +
𝑠 .

A large body of literature explored the modeling of 𝑘 𝑠 for fully rough flows based on the
roughness topography alone. However, there is no universal formula of 𝑘 𝑠 for arbitrary roughness
geometries. Recent reviews are provided by [69] and [68]. In most of these studies, empirical
correlations are developed based on single or a small number of surface features (slope, height
fluctuation intensity, skewness, etc.) and fitted to 𝑘 𝑠 data of a few surfaces [102; 103; 104; 105].
So far, the most advanced existing approach is based on machine learning (ML), e.g. [4] enabled
by large 𝑘 𝑠-geometry datasets. In the work of [4], 45 datasets with different roughness geometries

56

Figure 4.1 Wall-normal profiles of double-averaged velocity inside the RSL of spatially accelerating
(a) and suction-blowing (and separating, b) boundary layer flows [3], showing self-similarity of
mean velocity regardless of pressure gradients, when normalized by 𝑈𝑅.
are used to train deep neural network (DNN) and Gaussian process regression models that contained
more than 10 surface features as inputs. The ML model is shown to give significantly better overall
𝑘 𝑠 predictions for a large range of surfaces than existing empirical correlations refitted on these
surfaces. [106] optimized the training of neural networks with the help of Particle Swarm Optimized
Backpropagation technique to predict 𝑘 𝑠 for 86 surfaces. The data included 2D and 3D roughness
type with both regular and irregular arrangements. The roughness types include sandgrain, grit-
blasted, truncated cone and multiscale roughness. With the help of roughness height statistics (r.m.s.
of height fluctuations, skewness and kurtosis) and effective slope, the optimized algorithm had a
prediction accuracy compared to other existing empirical models with a mean absolute percentage
error of 28.1%. In addition, feature sensitivity analysis showed r.m.s. of height fluctuations to be
the most important feature followed by skewness, effective slope and kurtosis.

The limitation of ML models is, to a large extent, controlled by the generality of the datasets.

57

[107] showed that the DNN model of [4], as is, yielded large errors in a test with larger datasets
containing types of roughnesses not well represented during the ML model development. However,
after it is retrained on all data (including new ones), the DNN model outperformed physics-
based or empirical 𝑘 𝑠 correlations. Currently, the largest collection of 𝑘 𝑠 data on different rough
surfaces is the Southampton Roughness Database (https://roughnessdatabase.org). As
data is important for 𝑘 𝑠 prediction and can be costly to obtain, insights are needed on limitations
of currently available data and on the minimum set of important features affecting 𝑘 𝑠 (which
allows reduction of inputs and thereby data needs for ML training, as well as gaining fundamental
knowledge).

This chapter aimed to establish relations between geometrical features of roughness and (i) the
double-averaged mean velocity profile 𝑈 (𝑦) near the wall, as well as (ii) the equivalent sandgrain
height 𝑘 𝑠, for a wide range of roughnesses. For the 𝑘 𝑠 prediction, the complete Southampton
Roughness Database will be used; the database is not directly applicable to 𝑈 (𝑦) prediction,
however, as 𝑈 (𝑦) profiles are typically not reported; a subset of that database will be used instead.
Specific objectives include (1) extracting the important roughness features influencing 𝑈 and
𝑘 𝑠, and (2) using the above information to build data-driven regression models to predict 𝑈 (𝑦)
and 𝑘 𝑠 for a wide range of equilibrium, rough-wall turbulent channel flows, based on existing
and new DNS datasets of half-channel flows over rough surfaces. Since the datasets and detailed
algorithms used are different for the predictions of 𝑈 (𝑦) and 𝑘 𝑠, the chapter is divided into two
parts: Section 4.2 focused on roughness features important for 𝑈 (𝑦) and its prediction based
on these features, while Section 4.3 employed similar approaches (though varying in details) to
𝑘 𝑠 prediction and analyses.
In Section 4.2, Section 4.2.1 surveyed the roughness datasets used
for training the regression models, Section 4.2.2 introduced the machine learning methodologies,
and Section 4.2.3 identified the important roughness features and summarized the performance
of the chosen algorithms. Similar outline followed in Section 4.3 for 𝑘 𝑠 prediction and analyses.
Conclusions are provided in Section 4.4.

4.2 Prediction of double-averaged mean velocity profiles in RSL

𝑢 is not negligible compared to 𝑈 (e.g. ⟨(cid:101)

For the 𝑈 (𝑦) prediction, only the RSL region is considered. The RSL is typically defined as the
𝑢2⟩1/2/𝑈 > 0.1, similar to [86]). Under such
layer in which (cid:101)
definition, the sublayer thickness (𝑦𝑅) is usually slightly higher than the crest height (for example,
see [2; 108]. In this work, 𝑈 in the region 𝑦/𝑘𝑐 = 1 is examined, as 𝑘𝑐 is an approximate estimate
of the RSL thickness. Using 𝑘𝑐 instead of the exact 𝑦𝑅 removes the uncertainty of RSL-thickness
identification brought by the threshold value.

58

Figure 4.2 Visualization of rough surfaces used in DNS simulations to generate data for 𝑈 (𝑦)
prediction. Cases ‘Cxx’ are from [4] with a portion (𝛿 and 0.5𝛿 in 𝑥 and 𝑧, respectively) of the
domain shown. Cases ‘Ax’ are newly generated surfaces used for additional simulations herein,
with the complete domain shown (3𝛿 and 𝛿 in 𝑥 and 𝑧, respectively). Cases A1-A5 are similar to
C18, but with varying 𝐼𝑥. Cases A6 and A7 are similar to C05, with varying 𝐼𝑥.
4.2.1 RSL flow datasets

A total of 21 half-channel flows with different wall roughness topographies are used in training
the regression techniques for the prediction of 𝑈 (𝑦). The surfaces are shown in Figure 4.2. They
included a subset of the cases simulated by [4], indicated as ‘Cxx’ (using the same case number as
in [4]). Seven additional surfaces (A1 to A7) are synthesized to investigate the effects of element
inclination angle and porosity. All surfaces are used in DNS simulations of turbulent half-channel
flows with friction Reynolds Number 𝑅𝑒𝜏 = 1000 and domain dimensions (𝐿𝑥, 𝐿 𝑦, 𝐿𝑧) = (3𝛿, 𝛿, 𝛿),
in streamwise, wall normal and spanwise directions, respectively, where 𝛿 is the half-channel height.
𝐿𝑥 and 𝐿𝑧 satisfy the rough-wall minimal-span channel constraints defined by [109] and [110]. All
the cases are considered to fall in the fully rough regime based on the values of 𝑘 +
𝑐 ≫ 1, where 𝑘𝑐
is the roughness crest height measured from the trough and + denotes normalization by wall units
(𝑢𝜏 and 𝛿𝜈). The rough surfaces are formed by distributed ellipsoidal sandgrains similar to those of
[1]. The distribution, inclinations and size of the ellipsoids are systematically varied. Both regular
arrangement (with every element having the same rotation and size) and random arrangements
(with random orientations and random size) are included.

Based on the literature regarding the effect of roughness texture on mean momentum (predom-
inantly in terms of the hydrodynamic drag, 𝑘 𝑠) , a list of surface features are considered as relevant
for affecting the RSL mean velocity profile. Readers are referred to [69] for a comprehensive

59

summary of dependency of 𝑘 𝑠 on different roughness features. They included the average peak-
to-trough height (𝑘𝑡), effective slope (𝐸 𝑆𝑥), various moments of local height fluctuations from the
mean height (𝑅𝑎, 𝑘𝑟𝑚𝑠, 𝑠𝑘 , and 𝑘𝑢), roughness-element inclination angle (𝐼𝑥), porosity (𝑃𝑜), and
spatial correlation length (𝑅𝑙𝑥). They are defined as,

𝐸 𝑆𝑥 =

1
𝐴1

∫

𝑥,𝑧

(cid:12)
(cid:12)
(cid:12)

𝜕𝑘
𝜕𝑥

(cid:12)
(cid:12)
(cid:12)

𝑑𝐴,

𝑘𝑟𝑚𝑠 =

(cid:20) 1
𝐴1

∫

𝑥,𝑧

(𝑘′)2𝑑𝐴

(cid:21) 1/2

,

𝑠𝑘 =

𝑘𝑢 =

1
𝐴1𝑘 3

𝑟𝑚𝑠

1
𝐴1𝑘 4

𝑟𝑚𝑠

∫

𝑥,𝑧

∫

𝑥,𝑧

(𝑘′)3𝑑𝐴,

(𝑘′)4𝑑𝐴,

𝐼𝑥 = tan−1(cid:110) 1
2

𝑠𝑘

(cid:19) (cid:111)

,

(cid:18) 𝜕𝑘
𝜕𝑥

𝑃𝑜 =

1
𝐴1𝑘𝑐

∫ 𝑘 𝑐

0

𝐴2𝑑𝑦,

and

𝑅𝑙𝑥 = 𝑟𝑥 𝑅𝑥=0.37,

(4.3)

(4.4)

(4.5)

(4.6)

(4.7)

(4.8)

(4.9)

𝑥,𝑧

∫
𝑥,𝑧

𝑘 𝑑𝐴, and 𝑅𝑥 (𝑟𝑥) = ∫

𝑘′(𝑥, 𝑧)𝑘′(𝑥 + 𝑟𝑥, 𝑧) 𝑑𝐴/(𝑘 2

where 𝑘 (𝑥, 𝑧) is the distribution of roughness height (with 𝑘 = 0 at the roughness trough), 𝐴1(𝑦) is
the total planar area at a 𝑦 elevation, 𝐴2(𝑦) is the fluid area, 𝑘′ = 𝑘 − 𝑘 is the local height fluctuation
𝑟𝑚𝑠 𝐴1) is the
from the mean height, 𝑘 = 1/𝐴1
two-point correlation of 𝑘′ with streamwise separation 𝑟1. The exact threshold of 𝑅𝑥 = 0.37 is used
to remain consistent with the correlation-length calculations used in the Southampton Roughness
Database. The features with dimension of length (e.g. 𝑘𝑟𝑚𝑠, 𝑅𝑙𝑥 and 𝑘𝑡) are normalized by the
∫
first moment of 𝑘′, 𝑅𝑎 = 1/𝐴1
𝑥,𝑧 |𝑘′|𝑑𝐴. The list of features analyzed are similar to those studied
by [4], except for (i) the additional two-point characteristics, 𝑅𝑙𝑥, and (ii) considering only the
streamwise characteristics for the features that are directional-dependent (i.e. slope, inclination and
correlation length) for an arbitrary roughness. The latter change is made to reduce the number of
model inputs in light of the limited train-test data available.

To allow the ML scheme to efficiently detect the 𝑦-pattern of the velocity profile, the 𝑦 values
(normalized by 𝑘𝑐) at which 𝑈 values are evaluated, 𝑦/𝑘𝑐 ∈ [0, 1], are also used as input variables
along with the surface and flow features mentioned above. The DNS used refined 𝑦 mesh near

60

Figure 4.3 Mean velocity profile as a function of wall normal distance (𝑦) in all cases.
𝑦 = 0. To remove bias of datasets due to more data originating from the finer-resolution region
near 𝑦 = 0, the 𝑈 (𝑦) profiles are interpolated to an equally-spaced 𝑦 grid in the layer 𝑦/𝑘𝑐 ∈ [0, 1].
Figure 4.3 showed the 𝑈 (𝑦) profiles of all cases in the datasets. The normalization of 𝑈 (𝑦) leads
to all profiles crossing the point (1, 1); on the other hand, at 𝑦 = 0 (i.e. trough elevation) 𝑈 = 0
due to the no-slip condition. The profiles in blue are those that display near-wall flow reversal
(i.e. 𝑈 < 0). These cases corresponded to surfaces with relatively dense distribution of rough
elements with steep slopes. The rough surfaces in the datasets are surveyed in Figure 4.4 in terms
of the values of three features that are generally considered the most crucial for the roughness
effects on hydrodynamic drag by the literature (see reviews by [69; 68; 107]: 𝑘𝑟𝑚𝑠/𝑅𝑎, 𝑆𝑘 and 𝐸 𝑆𝑥.
The optimum number of clusters for data classification is chosen based on varying the number of
clusters from 1 to 7 and selecting the one with the highest silhouette score. Silhouette score is a
measure of how well a data point is represented by one cluster in comparison to its nearest cluster.
The two classes corresponded well with the two types of flows distinguished in Figure 4.3 based on
whether or not mean-flow reversion occurred inside the RSL. Blue circles corresponded to surfaces
with mean-flow reversions as shown in Figure 4.3 (as blue curves therein). It is shown here that
they are likely characterized by higher 𝐸 𝑆𝑥 and lower 𝑠𝑘 . On the other hand, red circles represented
surfaces without mean-flow reversions, corresponding to overall lower 𝐸 𝑆𝑥 and higher 𝑠𝑘 . Also,
it is shown that almost all present surfaces are characterized by 𝑠𝑘 > 0. In other words, very few
‘pitty’ surface is included. Figure 4.4 highlighted the need for 𝑈-profile data on pit-dominated
surfaces, as well as those with mild skewness and high r.m.s. (e.g. surfaces with similar peak and
dip geometries), and high skewness and low r.m.s (e.g. sparse distribution of peaky elements).

Besides the geometry-based features mentioned above, an additional flow-based feature is also

61

Figure 4.4 Clustering of rough surfaces in (𝑘𝑟𝑚𝑠,𝐸 𝑆𝑥,𝑠𝑘 ) space.
proposed, to examine whether inclusion of some degree of flow information significantly improves 𝑈
or 𝑘 𝑠 prediction. Given the importance of the sheltering phenomenon in roughness drag generation
(i.e. the reduction of drag generated by a roughness element due to its location inside the wake of
an upstream element; see, for example, [10; 100], the additional feature estimated the significance
of the sheltering phenomenon.

Figure 4.5 Illustration of the quantification of sheltering, defined as regions with 𝑢 ≤ 0. Contour
shows 𝑢+. Vertical dotted lines indicate locations used to calculate various frontal areas: leeward
sheltered frontal area (𝐴𝑙), windward frontal shelter frontal area (𝐴𝑤) and total frontal area (𝐴𝑡).
It is termed the sheltering parameter and is defined as

𝛼 =

(cid:205) 𝐴𝑙 − (cid:205) 𝐴𝑤
(cid:205) 𝐴𝑡

,

(4.10)

where 𝐴𝑙 and 𝐴𝑤 are the leeward and windward frontal areas, respectively, of one roughness
element and 𝐴𝑡 is the total frontal area of the element. These areas are illustrated in Figure 4.5.

62

Summation is then performed for each area across all roughness elements. The sheltered region
of rough surface is defined as that of the surface on which the time-averaged 𝑢 velocity is zero
or negative. In the limiting case of fully sheltered roughness elements, 𝐴𝑙 = 𝐴𝑤, rending 𝛼 = 0.
Therefore, 𝛼 is a measure of the inverse of sheltering.

4.2.2 ML techniques

The problem solved is a multi-output regression one from 𝑅𝑛 to 𝑅𝑚, where 𝑛 is the number of
roughness features plus the number of grid points of 𝑦 in the RSL and 𝑚 is the number of discrete
𝑈 (𝑦) values. For the prediction of 𝑈 (𝑦), both the non-dimensional input and output data are further
scaled by removing the mean and dividing by the standard deviation to uniformize the range of
dataset.

As a part of feature selection, the correlation between each set of two out of all feature are
calculated and shown in Figure 4.6. The diagonal plots show instead the p.d.f. of the distribution
of each individual feature. Significant correlations are found for the pairs of (𝑃𝑜, 𝑠𝑘 ), (𝑘𝑟𝑚𝑠, 𝑘𝑢) and
(𝑠𝑘 , 𝑘𝑢). The strong correlation within the pairs (𝑃𝑜, 𝑠𝑘 ) and (𝑘𝑟𝑚𝑠, 𝑘𝑢) as highlighted in the black
box may be explained by the fact that all surfaces are peaky ones (i.e. surfaces characterized by
𝑠𝑘 > 0 and formed by predominantly peaky elements with 𝑘′ > 0 instead of pits with 𝑘′ < 0) in the
datasets. The strong correlation between 𝑠𝑘 and 𝑘𝑢 may be attributed to database limitations. Based
on the correlations, seven input features from the initial list of ten are considered as independent
features describing the differences of present surfaces and are used to train the algorithm: 𝑅𝑙𝑥, 𝑠𝑘 ,
𝐸 𝑆𝑥, 𝑘𝑡, 𝑘𝑟𝑚𝑠, 𝛼, and 𝐼𝑥.

The rough surfaces are randomly divided into two sets: 85% data (18 cases) for training and
15% (3 cases) for testing. Gaussian Process Regression (GPR)—a non-parametric, probabilistic
prediction machine learning technique—is used. Two other regression techniques, KernelRidge
and Decision Tree Regressor are also used for comparison; in comparison, GPR is shown (in
Appendix 9) to give the lowest 𝐿∞ error (defined by Equation (4.11)). Also, GPR had the
capability of identifying non-polynomial shapes of 𝑈 (𝑦), which can be important in a multi-output
regression problem as the current one. The functioning of the algorithm is dependent on its hyper
parameters that include different kernels (quadratic or exponential squared), over-fitting/under-
fitting penalty parameter and the number of iterations. These parameters are fine-tuned using
the present datasets, by running the algorithm for 𝑁 = 1300 times over different training-testing
combinations for the chosen input feature combination. The optimal set of hyper parameters are
chosen based the maximum 𝐿∞ error across all flow cases. The optimized GPR algorithm used an
exponential squared kernel (i.e. radial basis function) to fit the data. It used 200 iterations to reach
convergence while fitting the data. After considering various values for over-fitting/under-fitting
penalty parameter in the range (10−3, 103), a value of 10−2 is found to offer the best performance.
Along with GPR, principal component analysis (PCA) (which will be introduced in the next section)

63

is also used to identify the minimum required number of input features needed for 𝑈 (𝑦) prediction.

The error of the mean velocity prediction using GPR is measured by two norm metrics,

Figure 4.6 Pair plots of roughness features used as inputs.

𝐿∞ = max
𝑦∗∈[0,1]

|𝑈∗

𝐷𝑁 𝑆 (𝑦∗) − 𝑈∗

𝐺𝑃𝑅 (𝑦∗)|

and

𝐿2 =

∫ 1

0

[𝑈∗

𝐷𝑁 𝑆 (𝑦∗) − 𝑈∗

𝐺𝑃𝑅 (𝑦∗)]2𝑑𝑦∗,

(4.11)

(4.12)

where 𝑈∗ = 𝑈/𝑈𝑘 𝑐 and 𝑦∗ = 𝑦/𝑘𝑐 are the dimensionless mean velocity and elevation. 𝑈𝐷𝑁 𝑆
and 𝑈𝐺𝑃𝑅 are, respectively, the velocity extracted from DNS and that predicted by GPR. The
respective averages of 𝐿∞ and 𝐿2 (calcualted across both train and test data) for the 𝑁 = 1300

64

random combinations of train-test split are used as the overall prediction errors. The values
of parameters for all surfaces and the code for regression analysis are made publicly available
(https://github.com/Mangavelli94/Prediction-of-wake-velocity-using-ML).

4.2.3 Results

Figure 4.7 Dependence of the prediction error norms (𝐿2 in orange and 𝐿∞ in blue) on the number
of features used as ML inputs. For each number of features, the set of features that yield the
minimum error among all combinations of the same number of features is listed.

Initially, the prediction of optimized GPR on training datasets with all input rough features is
measured using 𝐿∞. The prediction yielded 𝐿∞ < 0.1 for all 𝑁 random runs, indicating the ability
of the chosen algorithm in successfully learning the different patterns in the training datasets. Next,
to analyze how the number of roughness features used for GPR training and testing affects the
performance of the ML prediction, subsets of the 7 features are tested. The number of features in
these subsets are varied from 1 to 7, for each of which all possible combinations of features (out
of the total 7) are evaluated for all 𝑁 random train-test data splits. The minimum error norms and
the ‘best’ feature combinations (associated with the minimum 𝐿∞) as a function of the number of
features are shown in Figure 4.7(a). Both 𝐿∞ and 𝐿2 showed a decrease with the number of features
increasing from 1 to 5. This is expected as the complex dependency of the mean velocity profile on
roughness geometry is not captured by only a very small number of inputs; increasing the number
of inputs allowed capturing additional flow physics. Both norms reached a minima when four or
five features are used for training and testing. When more than five features are used, the error is
shown to slightly increase with the number of features. This is probably due to over-fitting, which
is expected with 21 datasets only. The best 4-feature combination is found to be (𝑅𝑙𝑥, 𝑆𝑘 , 𝐸 𝑆𝑥, 𝐼𝑥),
while the best 5-feature one contained an additional 𝑘𝑟𝑚𝑠. The finding suggested that 𝑅𝑙𝑥, 𝑘𝑟𝑚𝑠 (or
𝑆𝑘 ), 𝐸 𝑆𝑥 and 𝐼𝑥 are the most important features from the total 10 initial ones and that they each
contribute independently to the RSL mean velocity, for the present collection of rough surfaces. It

65

had been established that 𝑘𝑟𝑚𝑠, 𝑆𝑘 and 𝐸 𝑆𝑥 are important rough-wall characteristics that controls
the hydrodynamic drag imposed by a rough surface [68]. The present result suggested that these
characteristics are also important for RSL mean velocity distribution. In addition, the streamwise
inclination and two-point statistics (such as 𝑅𝑙𝑥) appeared to be also important. Note that 𝑘𝑟𝑚𝑠 and
𝑆𝑘 are expected to affect the flows differently and, as such, introducing an additional 𝑘𝑟𝑚𝑠 should
significantly reduce 𝐿∞ and 𝐿2 errors, as opposed to the insignificant changes observed herein.
This is due to limitations of the present database: the fact that (i) surfaces with 𝑆𝑘 < 0 are not
included and that (ii) 𝑘𝑟𝑚𝑠 and 𝑆𝑘 are strongly correlated by chance.

Figure 4.8 Dependence of 𝑈 (𝑦) prediction error on the number of principal components used as
ML inputs.
In order to verify the above analysis on the smallest number of features needed to model 𝑈 (𝑦),
another method, the principal component analysis (PCA), is used to reduce the dimensionality of
the feature space to extract the uncorrelated principal components(PCs). These components are
obtained using linear mapping of input features that best describes the variation of 𝑈 (𝑦) in the
datasets. Similarly, the number of principal components is allowed to vary from 1 to 7. The
GPR algorithm is re-trained for 𝑁 combinations of train-test split for varying number of principal
components as inputs. The maximum value of error norms (calculated as the maximum value
of errors for all rough cases) as a function of the number of principal components are shown in
Figure 4.8. Both error norms decreased overall from 1 to 5 principal components, reaching a
minimum at 5. This is consistent with the findings of the lowest errors observed for the use of 4 or
5 features for training/testing, as shown in Figure 4.7.

Now focus on the GPR algorithm trained with five principal components that yields the minimum
𝐿∞ error. To provide insights on what features are the major contributors to each principal
component (i.e. PC1 to PC5), Table 4.1 tabulated the weighted contribution, calculated as a product
of squared correlation of the feature and the fraction of variance of its principal component. The
correlation between a rough feature and a principal component measures how much a particular
roughness feature influenced the principal component. To quantify the importance of each of the

66

five PCs, the fractions of variance for each of the 5 principal components are listed in the caption
of Table 4.1. The fractions of variance in [PC1,PC2,PC3,PC4,PC5] are [0.5,0.22,0.16,0.06,0.025].
They quantified the amount of variance of input features captured by the corresponding principal
component, which sums up to 1 for all five principal components. Most of the variance of input
features is absorbed in the first three principal components (i.e. PC1 to PC3). Results suggested
that 𝑅𝑙𝑥, 𝑠𝑘 , 𝑘𝑟𝑚𝑠 and 𝐸 𝑆𝑥 captured the majority of the 𝑈 (𝑦) variations. This is consistent with the
observations made from Figure 4.7. In addition, 𝑘𝑡 and 𝐼𝑥 are also shown to play a role. Results
of both analyses described above (i.e. feature elimination and PCA) demonstrated that each of the
features 𝑅𝑙𝑥, 𝑆𝑘 , 𝐸 𝑆𝑥, 𝐼𝑥 and 𝑘𝑟𝑚𝑠 independently affected 𝑈 (𝑦) in the RSL.

𝑅𝑙𝑥/𝑅𝑎

𝑘𝑟𝑚𝑠/𝑅𝑎

𝑘𝑡

𝑆𝑘

ES𝑥

I𝑥

PC1
PC2
PC3
PC4
PC5

0.00735
0.0726
0.056
0.00192
0.00325

0.095
0.0242
0.00912
0.00216
0.000975

0.125
0.001078
0.00176
0.00018
0.00525

0.095
0.0000792
0.0128
0.0252
0.00005

0.09
0.000242
0.0128
0.0276
0.0019

0.008
0.0638
0.0624
0.00096
0.0035

𝛼

0.065
0.0572
0.00272
0.00138
0.00975

Table 4.1 Weighted correlations between individual features with individual principal components
for the 𝑈 (𝑦) datasets.

Figure 4.9 PDF of 𝑈 (𝑦) prediction errors in all cases of the GPR algorithm developed with the
optimal set of inputs: (𝑅𝑙𝑥/𝑅𝑎, 𝑠𝑘 , 𝐸 𝑆𝑥, 𝐼𝑥, 𝑘𝑟𝑚𝑠): (a) (100)𝐿∞, (b) (100)𝐿2.
The analyses above (Figures 4.7 and 4.8) discussed the error norms that are averaged across both
train and test datasets and across all 𝑁 = 1300 random train/test data splits (also called ‘runs’).
Below we discuss the distribution of the errors for test datasets only, across the 1300 runs, when
the 5 most important input features identified above (𝑅𝑙𝑥, 𝑠𝑘 , 𝐸 𝑆𝑥, 𝑘𝑟𝑚𝑠 and 𝐼𝑥) are used to train
the model. Here, the exclusive use of test data gives an estimation of the model’s performance
for predicting 𝑈 (𝑦) in flows with roughness not included in the training datasets. The probability
density function (PDF) of the error norms (in terms of their maximum and mean values among the
test datasets) are displayed in Figure 4.9 for 𝐿∞ (a) and 𝐿2 (b). Among the 1300 runs, only 26%
of the runs yielded a maximum 𝐿∞ higher than 0.1, with the peak of PDF at around 𝐿∞ = 0.06.

67

Similarly, 22% of the runs gave maximum 𝐿2 higher than 0.005, with the PDF peak at around
0.001. This result confirms that a majority of the test cases are well predicted.

Next, the sensitivity of the GPR prediction to the dataset size is tested by adding six more
cases (by varying the 𝑃𝑜 and 𝐼𝑥 to include new surfaces that are sparser or varying within a wider
range of inclination angles). These additional surfaces are shown in Figure 4.2 starting with ‘A’.
Among 2600 runs (i.e. 17% of all possibilities of picking 4 test cases from the total of 27 cases
with different surfaces; not all possibilities are tested due to due to long computational time), the
probabilities of the maximum of 𝐿∞ being higher than 0.1 and the maximum of 𝐿2 being higher
than 0.005 increased to 31% and 23% from 26%, 22% respectively. The slight increase of the
occurrences of high-error cases indicates that the current algorithm is only weakly sensitive to the
addition of new datasets generated using the identified important features.

Figure 4.10 Examples of predicted mean velocity profiles compared to their DNS-extracted profiles:
(a) low-error cases with 𝐿∞ = 0.033 (red), 0.027 (blue) and 0.028 (black); (b) relatively high-error
cases with 𝐿∞ = 0.075 (red), 0.14 (blue) and 0.09 (black).
Figure 4.10 shows the performance of 𝑈 (𝑦) predictions at different 𝐿∞ errors. The DNS-extracted
𝑈∗(𝑦∗) and the GPR predictions of the three test cases are compared in the figure, for three
representative runs that yielded different maximum 𝐿∞ values of 0.03, 0.08 and 0.14, respectively.
Figure 4.10(a) shows three cases where the predictions of GPR are very close to the DNS-extracted
profiles, while Figure 4.10(b) shows another three cases where the errors are relatively larger.
The errors appear to be significantly attributed to the velocity profile in the lower portion of the
roughness sublayer, especially if a mean-flow reversal (i.e. 𝑈 (𝑦) < 0) layer is present. To analyze
whether some surfaces consistently lead to higher prediction errors than other surfaces among the
1300 runs, the number of times a particular case having 𝐿∞ > 0.1 (while being a test case) is shown
in Figure 4.11(a). Note that all cases are used as test cases for similar total occurrences for this
analysis, as almost all possible combinations of test cases are analyzed. It shows that C04, C06
and C14 have a higher frequency of 𝐿∞ > 0.1. The corresponding DNS-extracted 𝑈 (𝑦) profiles of
these three cases are highlighted in Figure 4.11(b) in black and compared to other profiles.

68

Figure 4.11 (a) Occurrences of each individual flow cases with 𝐿∞ > 0.1, when the set of optimal
input features are used. (b) DNS-extracted mean velocity distributions of all cases, with the three
),
cases with the highest large-error occurrences in (a) highlighted in black: C04 (
C14 (

), C06 (

).

Cases C04 and C06 are the only two surfaces with uniform streamwise inclination for all grains in
the datasets, while in Case C14 rough elements are rather sparsely located.The under-representation
of inclined and sparse surfaces in the database may explain the consistently higher prediction errors
for these three surfaces.

4.2.4 Explanation of the importance of inclination and porosity

The spatial distribution of rough peaks is shown by the literature for equilibrium flows [40] and
𝑢. Hence, 𝑃𝑜 (or 𝑅𝑙𝑥) that
non-equilibrium flows in Chapter 3 to impact the formation and growth of (cid:101)
characterizes the overall distance between rough peaks is studied in this section. Additionally, the
observed importance of 𝐼𝑥 on 𝑈 (𝑦) is not well understood and is also studied here. Figure 4.12 and
Figure 4.13 examine the iso-surfaces of 𝑢(𝑥, 𝑦, 𝑧) with small magnitudes, which shows the wake
region with a small positive value (𝑢+ = 1.0) and the separated-flow region with a negative value
(𝑢+ = −0.05) around each roughness element, in a few cases with different element inclinations
and distribution densities (i.e. porosity). Comparison between Figure 4.12(a) and Figure 4.12(b)

69

Figure 4.12 Effect of roughness-element inclination on 𝑢 distribution in the RSL: (a) Case C04
with elements inclined toward the direction of flow and (b) Case C06 with elements inclined in
the opposite direction of flow. Iso-surfaces of 𝑢+ = 1.0 and 𝑢+ = −0.05 are shown in yellow and
magenta, respectively, overlaid on visualized rough surfaces.

shows that an inclination of roughness elements toward the direction of flow (C04) led to smaller re-
circulation regions in the wake of roughness (shown by the magenta surfaces) compared to case C06
where the inclination is toward the opposite direction of the flow. The larger re-circulation regions
in case C06 are because, in this case, the 𝑢 = −0.05 iso-surfaces grow wider (in 𝑧) downstream of
each roughness element, as compared to case C04. Therefore, the RSL 𝑈 (𝑦) profile takes lower
values in C06 than in C04, which is shown in Figure 4.14(a) where case C06 shows a stronger
re-circulation region closer to the wall. In addition, channeling regions of locally high-speed flows
between streamwise ridges of roughness elements are seen here for all cases, similar to other studies
in the literature for roughness of various geometries, such as irregular surfaces with Gaussian height
distribution [111], sand grain roughness [112], pyramid roughness [113], among others. Around
70% of surfaces in the database are randomly distributed and/or rotated. Figure 4.13 compares two
cases with the same roughness elements but different values of distribution density. The higher
porosity in C18 (Figure 4.13(b)) leads to more significant fraction of the in-plane area

70

Figure 4.13 Effect of roughness porosity on 𝑢: (a) C05 case with densely distributed elements,
and (b) C18 case with the same 𝑘𝑐 as C05 but a sparser distribution. Iso-surfaces of 𝑢+ = 1.0 and
𝑢+ = −0.05 are shown in yellow and magenta, respectively, overlaid on visualized rough surfaces.

Figure 4.14 Mean velocity profiles for (a) the two cases with different 𝐼𝑥 (shown in Figure 4.12)
and (b) the two cases with different porosity values (shown in Figure 4.13).

In addition, the larger streamwise separation between
covered by high-speed channeling flows.
roughness elements reduces the fraction of the surface covered by the separation regions.
In
contrast, a more densely distributed roughness shown in Figure 4.13(a) displays separation regions
growing along 𝑥 and filling the streamwise distance between neighboring elements. These effects
of porosity on 𝑈 (𝑦) is shown in Figure 4.14(b), where near-wall reversal of the plane-average
velocity occurs in C05 but not in C18.

71

For surfaces mentioned above, both 𝐼𝑥 and parameters measuring the roughness peak distribution
(such as 𝑅𝑙𝑥 and 𝑃𝑜) are shown to affect the these separation regions. This suggests the importance
of 𝐼𝑥 and 𝑃𝑜 (as well as sparsity) on the prediction of 𝑈 (𝑦). Systematic analyses with controlled
surface variations is needed to understand the interplay between the effects of these features.
4.3 Prediction of 𝑘 𝑠
4.3.1 Datasets and ML technique

For the prediction of 𝑘 𝑠/𝑅𝑎, additional data from the Southampton Roughness Database [114;
67; 103; 115; 116] are added to the database used by [4]. A total of 90 rough surfaces are included.
A large portion of the additional datasets represent cuboid or bar roughness with constant or varying
heights which are not well represented in the datasets of [4]. The uncertainty in the estimation of
friction velocity is at most 1% for the present DNS cases. The uncertainty (for available data) in
𝑠 ) + 𝐶
experimental datasets is ≤ 3% [102]. A rough estimation based on Δ𝑈+ = (1/𝜅) log10(𝑘 +
(where 𝐶 is a constant), assuming that Δ𝑈+ is measured at a location with 𝑈+ = 10, shows that a
1% higher 𝑢𝜏 estimate is equivalent to a 10% higher 𝑘 𝑠. Though the calculation of 𝑘 𝑠 is not linearly
correlated to friction velocity, the uncertainty is hypothesized to be in similar range.

The roughness features used as inputs consist of single-point height statistics: 𝐼𝑥, 𝑃𝑜, 𝑠𝑘 , 𝑘𝑢,

𝑘𝑟𝑚𝑠/𝑅𝑎, 𝐸 𝑆𝑥 and 𝐸 𝑆𝑧, where

𝐸 𝑆𝑧 =

1
𝐴1

∫

𝑥,𝑧

(cid:12)
(cid:12)
(cid:12)

𝜕𝑘
𝜕𝑧

(cid:12)
(cid:12)
(cid:12)

𝑑𝐴,

(4.13)

𝐼𝑧 = tan−1(cid:110) 1
2
and two-point correlations of height fluctuations: 𝑅𝑙𝑥 defined in Section 4.2.1 and 𝑅𝑙𝑧 defined as

𝑠𝑘

,

(4.14)

(cid:19) (cid:111)

(cid:18) 𝜕𝑘
𝜕𝑧

𝑅𝑙𝑧 = 𝑟𝑧 𝑅𝑧=0.37,

(4.15)

𝑥,𝑧

𝑘′(𝑥, 𝑧)𝑘′(𝑥, 𝑧 + 𝑟𝑧) 𝑑𝐴/(𝑘 2

where 𝑅𝑧 (𝑟𝑧) = ∫
𝑟𝑚𝑠 𝐴1). The two coherence lengths are both nor-
|𝐼𝑧 |) is added, as the effect of spanwise inclination of
malized by 𝑅𝑎. The magnitude of 𝐼𝑧 (i.e.
roughness elements on the total drag is expected to be identical for inclinations toward positive or
negative 𝑧 directions. In addition, products of primary features and their second-degree polynomi-
als: (𝐸 𝑆𝑥)(𝐸 𝑆𝑧), (𝐸 𝑆𝑥)(𝑠𝑘 ), (𝐸 𝑆𝑧)(𝑠𝑘 ), (𝐸 𝑆𝑥) (𝑘𝑢), (𝑠𝑘 ) (𝑘𝑢), (𝐸 𝑆𝑧) (𝑘𝑢), 𝐸 𝑆2
𝑘 ) are
also included based on the importance of their primary features in characterizing rough surfaces
in the dataset. Inclusion of these quantities are shown to improve the current algorithm. However,
inclusion of 𝑅2
𝑙𝑧,etc. did not improve the algorithm significantly and they are not included. In
total, 19 features are used.

𝑧 and 𝑠2

𝑥, 𝐸 𝑆2

𝑙𝑥, 𝑅2

Figure 4.15 shows the pair plots between all roughness features and the actual 𝑘 𝑠/𝑅𝑎 values

measured from physical or numerical experiments. Specifically, the diagonal plots show the

72

Figure 4.15 Pair plots of roughness features and 𝑘 𝑠/𝑅𝑎. Orange and blue circles denote roughness
types: bar roughness and non-bar roughness, respectively.
probability density function of each single feature or 𝑘 𝑠; the off-diagonal plots show correlations
between each pair of two different roughness features or 𝑘 𝑠. In general, the bar-type roughness
tends to give higher skewness and 𝑘𝑟𝑚𝑠/𝑅𝑎 values due to the ideal nature of the roughness-element
geometry. Almost all features are shown to vary for a rather wide range, suggesting that the database
is overall comprehensive for the features shown. Wider ranges may be explored in the future for
𝐼𝑧 and 𝑘𝑢. The PDF of 𝑠𝑘 shows that most surfaces are positively skewed (i.e. peaky), while a
significant number of pitty surfaces (with 𝑠𝑘 < 0) are also included.

It is shown in Figure 4.15 that some correlation exists within the feature pairs (𝑃𝑜, 𝑆𝑘 ),
(𝐸 𝑆𝑥, 𝐸 𝑆𝑧), (𝑘𝑟𝑚𝑠, 𝑘𝑢), and (𝑠𝑘 , 𝑘𝑢), as highlighted in black boxes. The correlation between 𝑃𝑜 and

73

𝑆𝑘 is possibly because most of the surfaces are peaky; in this case, a higher skewness indicates a
sparser surface marked with higher porosity. The correlation between 𝐸 𝑆𝑥 and 𝐸 𝑆𝑧 is due to most
surfaces being isotropic in 𝑥 and 𝑧 (i.e. with mostly mild inclinations in 𝑥 and 𝑧, as shown by the
PDFs of 𝐼𝑥 and 𝐼𝑧). The correlation between 𝑘𝑟𝑚𝑠 and 𝑘𝑢 is because both features characterize the
intensity of 𝑘′ fluctuations. Lastly, the correlation between 𝑠𝑘 and 𝑘𝑢 is an interesting observation,
possibly due to more dominant peaks in sparser surfaces (i.e. higher 𝑠𝑘 ) contributed to significant
𝑘𝑢. In addition, the last row of plots show that 𝑘 𝑠/𝑅𝑎 does not correlate perfectly with any single
roughness feature.

The re-scaling of the input and output data is done similarly to that described in Section 4.2.2.
The Neural Network (NN) used by [4] is re-calibrated and retrained for the two-fold larger database
and the new set of inputs include both single- and two-point statistics. The hyperparameters of the
algorithm that are optimized based on the datasets include the number of neural layers and their
branches. Three layers are chosen with the numbers of branches in each of the three layer varying
from 14 to 21, 10 to 16, and 4 to 10, respectively. The over/under-fitting penalty parameter is
set to 0.2. The algorithm is trained over 250 iterations with a random selection of training and
testing data and the one with the lowest 𝐿∞ error (across both training and testing sets) is chosen.
A larger number of iterations for algorithm optimization is tested. The maximum error (with the
error defined as |𝑘 𝑠 − 𝑘 𝑠,𝑝𝑟𝑒𝑑 |/𝑘 𝑠,𝑝𝑟𝑒𝑑, where 𝑘 𝑠 and 𝑘 𝑠,𝑝𝑟𝑒𝑑 are the actual and predicted equivalent
sandgrain height values, respectively) changes by less than around 1% only, when the number of
iterations increases to 500 (while the optimization time increases by 3.5 times, from 7 hours to 24
hours). In every iteration, the training-to-testing-data ratio is maintained at 70/30. Additionally,
PCA is used to perform dimensional reduction to identify the minimum number of independent
input features needed for 𝑘 𝑠 prediction.

4.3.2 Results

Figure 4.16 shows the performances of the NN trained with the set of single-point surface
statistical inputs used by [4] (denoted as 𝑁 𝑁1) over the current dataset. Figure 4.16(a) compares the
predicted and actual 𝑘 𝑠 values for all cases, with the percentage differences shown in Figure 4.16 (b).
The maximum of percentage difference for all cases is around 75% and the mean error percentage
is around 11%. In comparison, the maximum and mean values of such error reported by [4] for a
database of 45 rough cases is around 30% and 10%, respectively. The lower errors reported there
is attributed to the limitation of the datasets used therein, with a lack of regular cuboid or bar-type
roughness, which are significantly different from multiscale surfaces or those consisting of non-bar
elements, associated with the higher skewness, 𝑘𝑟𝑚𝑠 and porosity values shown in Figure 4.15. In
comparison, Figure 4.17 shows the performance of the re-configured NN (denoted as 𝑁 𝑁2) with
hyperparemeters in the same range as the previous algorithm, but with the two additional two-point
surface statistical inputs (𝑅𝑙𝑥 and 𝑅𝑙𝑧).

74

Initially, similar to the analysis of identifying whether the algorithm is trainable on the rough surface
datasets in Section 4.2.3, the 𝐿′

∞ error for training datasets with all 19 rough features yielded

Figure 4.16 Performance of the re-calibrated NN with the same inputs as [4]. (a) Predicted values
(𝑘 𝑠,𝑝𝑟𝑒𝑑) compared to actual values (𝑘 𝑠) of the equivalent sandgrain height. (b) Percentage error.
Blue and red symbols denote numerical and experimental measurements, respectively.

Figure 4.17 Performance of the new NN with additional inputs of correlation lengths. (a) Predicted
values (𝑘 𝑠,𝑝𝑟𝑒𝑑) compared to actual values (𝑘 𝑠) of the equivalent sandgrain height. (b) Percentage
error. Blue and red symbols denote numerical and experimental measurements, respectively.
values up to 0.1. Then, the maximum and mean percentage errors are calculated as shown in
Figure 4.17(b). They are around 33% and 5%, respectively, reducing by almost a half from
predictions based on single-point inputs alone. The comparison highlights the importance of
including surface features that contain structural information of the surface topography besides
single-point height characteristics.

To compare the error distribution of the predictions of the 𝑁 𝑁1 and 𝑁 𝑁2 networks, the error
probability density functions (PDF) are shown in Figure 4.18. The higher probabilities of lower
errors for 𝑁 𝑁2 (with the PDF skewed to the left) confirms the overall better performance of the
model with the additional 𝑅𝑙𝑥 and 𝑅𝑙𝑧 inputs. To compare the two PDFs quantitatively, the 𝐿1 and
𝐿2 norms are calculated for each algorithms as

𝐿1 =

1
𝑛

𝑛
∑︁

𝑖=1

|𝑘 𝑠,𝑝𝑟𝑒𝑑,𝑖 − 𝑘 𝑠,𝑖 |
𝑘 𝑠,𝑖

,

75

(4.16)

𝐿2 =

1
𝑛

(cid:118)(cid:116) 𝑛
∑︁

𝑖=1

|𝑘 𝑠,𝑝𝑟𝑒𝑑,𝑖 − 𝑘 𝑠,𝑖 |2
𝑘 2
𝑠,𝑖

,

(4.17)

where 𝑛 is the total number of rough-wall flow cases. (𝐿1, 𝐿2) for 𝑁 𝑁1 are (10.7, 2.2) and for 𝑁 𝑁2
are (5.3, 1.0), again confirming the better prediction of 𝑁 𝑁2.

Figure 4.18 Comparison of error distributions of the NN predictions with (red) and without (blue)
correlation lengths as inputs. 𝐿′

∞ ≡ |𝑘 𝑠,𝑝𝑟𝑒𝑑 − 𝑘 𝑠 |/𝑘 𝑠.

The importance of correlation lengths is due to their direct representation of the overall distance
between the roughness peaks in 𝑥 and 𝑧. In Section 4.2, the formation of separation region around a
roughness element and its downstream growth in the element wake is shown to be dependent on the
distribution of roughness peaks. Since the total pressure drag is determined by the flow separation
zones, parameters describing the roughness peak distribution such as 𝑅𝑙𝑥 and 𝑅𝑙𝑧 are expected to
be important in determining 𝑘 𝑠. Inclusion of additional flow cases in the datasets including varying
peak element density for truncated cones [95], cuboidal roughness [115; 116], among others further
necessitates the need to include the correlation lengths as inputs.

Next, the minimum set of important topographical features for accurate predictions of 𝑘 𝑠 is
studied using principal component analysis (PCA). Instead using the 19 surface features as ML
inputs, principal components of the roughness-feature space are identified and used directly as ML
inputs to train the NN network. The number of principal components are systematically varied
from 1 to 19 (i.e. up to the number of total surface features). Figure 4.19(a) shows the variation of

76

the 𝐿∞-norm of the 𝑘 𝑠 prediction error, defined as

𝐿∞ = max

|𝑘 𝑠,𝑝𝑟𝑒𝑑,𝑖 − 𝑘 𝑠,𝑖 |
𝑘 𝑠,𝑖

,

(4.18)

as a function of the number of principal components used. The error decreases with increasing
number of principal components, as expected due to the inclusions or progressively richer

Figure 4.19 (a) Variation of 𝑘 𝑠 error norm with the number of principal components used as NN
inputs. (b) Comparison between the PDFs of 𝑘 𝑠 percentage errors predicted with 3 and 10 principal
components.
topographical information, roughly plateauing at 3 principal components and reaching a minima at
around 10 components. This suggests that, while a large number of topographical features (such
as 𝑘𝑟𝑚𝑠, 𝑆𝑘 , etc.) are relevant to drag generation, there are roughly three ‘orthogonal features’
only. Further increasing the number of principal components only slightly improve prediction.
Figure 4.19(b) compares the PDFs of 𝑘 𝑠 percentage errors for the predictions with 3 and 10
principal components. PDF distributions of both errors are skewed to the left for both the PCA
cases, indicating that most flow cases are successfully predicted with low errors. Predictions with
10 PCs give low errors in more cases.

To analyze the mapping between individual topographical features and the principal compo-
nents, Tables 4.2 and 4.3 show the weighted contributions of the features to each of the identified
principal components, for the PCA analyses with 3 and 10 PCs, respectively. The weighted contri-
butions are calculated as the product between the fraction of variance (defined in Section 4.2.3) of
the principal components and the squared correlation (also defined in Section 4.2.3). The fractional
variance of the principal components for the analysis with 3 PCs are (0.39, 0.23, 0.12), respectively,
covering 74% of variance of the topographical features, whereas the first six principal components
for the analysis with 10 PCs are (0.39, 0.23, 0.12, 0.08, 0.05, 0.05), respectively, covering 92%
of variance of the features. This explains the better overall prediction when using 10 principal
components than that using 3.

77

In addition, Table 4.2 and Table 4.3 show that, regardless of the number of PCs, the significant
contributors to the PCs (based on a weighted contribution ≥ 0.01) include: (i) 𝑘𝑟𝑚𝑠/𝑅𝑎, 𝑠𝑘 and 𝑘𝑢,
(ii) porosity, (iii) slope in both 𝑥 and 𝑧, and (iv) correlation lengths. Secondary features (products
and polynomial functions of slope and 𝑠𝑘 , not shown) also contribute significantly. The fact that
the NN predictions with 3 PCs are almost as good as those with more PCs (Figure 4.19(a)) suggests
that these features form the minimal set of important features for 𝑘 𝑠 prediction. As the number
of PCs increases from 3 to 10, the inclination angles in 𝑥 and 𝑧 are also added to the contributing
features. As the PCA illustrates the dominant variations of the rough surfaces in the Southampton
Roughness Database, the present findings suggest that almost all features studied herein are well
varied in the database, although not much so for the element inclination angles. This is indeed the
case as most studies on 𝑘 𝑠-correlation development (e.g. [103; 19; 67; 69]) focused on surfaces
with systematically varying single-point roughness height moments, slopes and element distribution
density, while systematically varied inclination angles were less examined (the present work and
[4]. Yet, the inclination angles are shown to affect the flow development in the roughness wake in
Section 4.2.4 and consequently the total drag. This explains the better prediction results of the NN
network with 10 PCs (where 𝐼𝑥 and 𝐼𝑧 contribute to the PCs) than that with 3 PCs (where 𝐼𝑥 and 𝐼𝑧
do not contribute to PCs). Future numerical/experimental datasets with wider range of inclinations
are needed to improve the database for 𝑘 𝑠 prediction development.

78

Surface

𝑘𝑟𝑚𝑠/𝑅𝑎

𝐼𝑥

𝐼𝑧

𝑃𝑜

ES𝑥

ES𝑧

𝑠𝑘

𝐾𝑢

𝑅𝑙𝑥/𝑅𝑎

𝑅𝑙𝑧/𝑅𝑎

PC-1
PC-2
PC-3

0.02
0.015
0.014

0.003

0.0003
0.00002 0.00046
0.0024
9.6E-08

0.03
0.02
0.015
0.016
0.0085 0.0021 0.00036 0.0038

0.04
0.0058

0.04
0.01

0.02
0.019
0.011

0.002
0.011
0.00038

0.0003
0.011
1.2E-05

Table 4.2 Weighted correlations between individual features with individual principal components for the 𝑘 𝑠 datasets, when three principal
components are used.

79

Surface

𝑘𝑟𝑚𝑠/𝑅𝑎

𝐼𝑥

𝐼𝑧

𝑃𝑜

ES𝑥

ES𝑧

𝑠𝑘

𝐾𝑢

𝑅𝑙𝑥/𝑅𝑎

𝑅𝑙𝑧/𝑅𝑎

PC-1
PC-2
PC-3
PC- 4
PC-5
PC-6

0.02
0.0155
0.014
0.00049
0.00034
1.2E-05

0.016
0.0003
0.0033
0.016
0.00002 0.00046
0.0024
9.6E-08
0.0085
0.00012 8.96E-07
0.0028
0.00048
0.0042
6.9E-05
0.044

0.036
0.0037

0.039
0.0058
0.0021
0.0021
4.2E-05
4.1E-07

0.031
0.039
0.015
0.01
0.00036
0.0038
2.42E-05 0.00048
0.00021
3.6E-05
3.7E-09
2.2E-06

0.0208
0.019
0.011
0.00070
2.5E-05
4.6E-05

0.0019
0.011
0.00038
0.032
0.00023
0.0011

0.00031
0.012
1.2E-05
0.030
3.7E-06
0.00086

Table 4.3 Weighted correlations between individual features with individual principal components for the 𝑘𝑠 datasets, when ten principal 
components are used.

80

4.4 Conclusions and discussions

This chapter examines the dependence of factors important for roughness modeling and the
roughness-sublayer mean velocity profile 𝑈 (𝑦) and equivalent
dispersive-stress transport (i.e.
sandgrain height 𝑘 𝑠) on the roughness features, as well as ML modeling based on a state-of-the-art
rough-wall-flow datasets. First, for the prediction of 𝑈 (𝑦), 21 DNS datasets are used. A data
clustering analysis of the different rough surfaces shows two groups, which corresponds to two
distinctive types of 𝑈 (𝑦) with or without near-wall flow reversal. Values of the surface features
demonstrate strong correlations among the pairs: (𝑃𝑜, 𝑠𝑘 ), (𝑘𝑟𝑚𝑠, 𝑘𝑢) and (𝑠𝑘 , 𝑘𝑢), attributed to
dataset limitation including the lack of pitty roughness. The observed correlations were used to
guide the reduction of feature inputs, to avoid over-fitting due to the limited datasets. One of the
novelties of the current work is the feature sensitivity analysis based on systematic elimination of
surface features by GPR and dimensionality reduction based on PCA. They identified 5 independent
surface features: (𝑅𝑙𝑥, 𝑠𝑘 , 𝐸 𝑆𝑥, 𝑘𝑟𝑚𝑠, 𝐼𝑥) as the most important ones. Detailed studies of time-
averaged velocity (𝑈 (𝑥, 𝑦, 𝑧)) patterns shows that the formation and growth of local separation
zones in the wake of roughness elements depend on the inclination and distribution of roughness
elements, explaining the observed importance of (𝐼𝑥) and 𝑅𝑙𝑥. The current analyses also identified
surfaces C04 and C06 (differing in 𝐼𝑥) and C14 (sparsely distributed roughness) to be difficult
to model. Therefore, additional data with systematically varying inclination angles and sparsity
(while keeping other features constant) are needed to make to improve 𝑈 (𝑦) predictions.

Next, for the prediction of 𝑘 𝑠, around 90 rough-wall-flow datasets from the Southampton
Roughness Database are used. Compared to past ML 𝑘 𝑠-prediction work [4], the present work
uses significantly larger datasets with (with bar-type roughness and pitty surfaces that were not
well represented) In addition, additional two-point statistical features of roughness (i.e. spatial
correlation lengths) are used as inputs, due to their capability in describing roughness sparsity
which affects downstream growth of local separation zones. Results show that the inclusion of
correlation lengths significantly improves NN algorithm predictions, reducing 𝐿1 and 𝐿2 errors
by a half compared to NN without these inputs. A PCA analysis shows that most of the surface
variation is contained in three principal components only. Comparisons between PCA analyses
with different numbers of principal components, as well as NN predictions based on these prin-
cipal components as inputs, show that the minimal set of important features for 𝑘 𝑠 predictions
is (𝑘𝑟𝑚𝑠, 𝑠𝑘 , 𝑘𝑢, 𝐸 𝑆𝑥, 𝐸 𝑆𝑥, 𝑅𝑙𝑥, 𝑅𝑙𝑧) and some of their combinations, which are overall consistent
with those observed for 𝑈 (𝑦) prediction due to the common physics (e.g. local separation zones)
affects both double-averaged velocity and pressure drag. It is also shown that the inclination angles
affect 𝑘 𝑠 but are not widely varied in the database. Systematic variations of individual features
(e.g. varying inclination angles of roughness elements of a given shape, varying correlation lengths
while keeping porosity constant if possible) are needed to further understand effects of the features

81

on 𝑘 𝑠 and to improve the database.

The work in this chapter demonstrates the capability of ML algorithms in predicting quantities
that are important for non-equilibrium turbulence modeling. For example, the predicted 𝑈 (𝑦)
profile can be used to model roughness-sublayer velocity in the framework of [75]. Results
also demonstrate the limitation of current available datasets and offer insights for future efforts
in systematic generation of rough-wall datasets. Future work based on the PCA analyses may
explore explicit expressions of the few orthogonal features based on the large number of surface
features used by the community, for effective feature reduction that will benefit roughness modeling.
Additionally, the standard deviation feature of GPR that measures the regions of high prediction
errors will be explored for 𝑈 (𝑦) prediction. This will offer further insights into the sensitivity
of 𝑈 (𝑦) on various features and what surface features are not sufficiently explored in the present
roughness database.

82

CHAPTER 5
CONCLUSIONS AND FUTURE WORK

5.1 Overall conclusions

Understanding the effect of roughness topography on the turbulent flow responses to temporal
and/or spatial acceleration or deceleration is important in many engineering and environmental
applications, such as flows in turbomachinery, around vehicles or sea-going vessels, around city
canopies, etc. Rough surfaces seen in real-world applications are usually irregular and multiscale,
adding an additional complexity in studying their effects. The current work aims to understand
the effects of wall roughness in non-equilibrium transient wall turbulence and their modeling for a
wide range of different roughness topographies.

In Chapter 2, direct numerical simulations of the half-channel turbulent flows are carried out
to study the effects of roughness texture on turbulence response to a sudden increase of bulk
velocity. The motivation is to quantitatively and qualitatively study the combined effect of the
mean flow acceleration and roughness texture, especially for the multiscale roughness. Two rough
surfaces: a homogeneous, densely distributed sandgrain roughness and a multiscale turbine-blade
roughness are simulated, along with a smooth-wall baseline case. The rough-wall flows develop
from a transitionally rough regime to a fully rough one during the transient. Immersed boundary
(IB) method is used to impose non-slip/penetration boundary conditions at the rough surface,
which is well-resolved by the grid. As seen in the existing literature of temporally and spatially
accelerated flows, turbulence over a smooth wall undergoes a quasi-laminarization process, which
is not observed in the presence of roughness. For rough cases, the rapid increase of hydrodynamic
drag and roughness sublayer (RSL) turbulence intensities in their early stages is attributed to high
turbulence production due to form-induced velocity gradients. This leads to a rapid increase in
turbulence kinetic energy (TKE) and isotropization of the RSL Reynolds stresses. Moreover,
the roughness geometry is shown to affect the recovery rate of the RSL turbulence. A slower
recovery of TKE redistribution (due to a slower amplification of pressure fluctuations) is seen on
the multiscale roughness, compared to the homogeneous one. This is associated with sparsely
distributed roughness peaks, producing new turbulence which takes time to spread and cover the
whole domain. For both rough cases, the outer-layer Reynolds stresses showed a slower recovery.
Additionally, the friction velocity (𝑢𝜏) is shown to scale on the mean velocity at the upper edge of the
roughness sublayer (𝑈𝑅). Also, the wall-normal profiles of dispersive stresses (when normalized
using 𝑈𝑅) are shown to be self-similar for both rough cases.

Motivated by the importance of form-induced fields in the dynamics of rough-wall turbulent
flows, Chapter 3 examines the effects of form-induced quantities on turbulence structure and the
TKE-redistribution process, important in determining non-equilibrium turbulence response, using
the same DNS half-channel flow data of Chapter 2. Analyses of two-point turbulent velocity

83

correlations show a rapid decrease in streamwise coherence and increases of inclination angles in
response to the acceleration in both rough cases, while the opposite is observed on the smooth wall.
In the rough cases, coherent motions of turbulence undergo a reduction in their streamwise length (as
shown by velocity spectra) in the homogeneous roughness case, whereas they are initially elongated
in 𝑥 on the multiscale rough wall, a phenomenon similar to the smooth-wall flow. This is further
verified by the higher local time-scale ratio (i.e. ‘frozen’ turbulence, with respect to the time scale
of the form-induced velocity shear) and their occurrence in large areas between sparsely distributed
peaks of the multiscale roughness. Furthermore, in addition to the form-induced production of
TKE as often discussed in the literature, a new pathway in which the form-induced velocity affects
the Reynolds stress transport through augmenting pressure fluctuations is identified. Analyses of
𝑖 gradients are
the pressure Poisson equation shows that both sources due to (cid:101)
important, with the spatial distribution of the latter correlated with the former, probably due to the
𝑢𝑖 in
production of TKE as a result of the form-induced fields. Results highlight the importance of (cid:101)
the dynamics of non-equilibrium turbulent flows.

𝑢𝑖 gradients and 𝑢′

To better understand how the transport of form-induced stresses depend on the roughness
texture, Chapter 4 analyzes and models the surface dependence of two quantities that are important
for the production of form-induced stresses (i.e. WKE): 𝑈 (𝑦) in RSL and a parameter measuring
the pressure drag, 𝑘 𝑠. A second objective is to model 𝑈 (𝑦) for different roughness topography to
be used for non-equilibrium turbulence modeling. For the prediction of 𝑈 (𝑦), 21 DNS datasets
of rough surfaces made up with ellipsoid elements in the fully rough regime are used. Data
classification analysis using clustering technique based on rough features shows two patterns in
𝑈 (𝑦), with and without mean-flow reversal. Feature sensitivity analysis by manual elimination of
inputs of the GPR algorithm identifies 5 features as the minimum set of important features in the
prediction of 𝑈 (𝑦): 𝑅𝑙𝑥, 𝐼𝑥, 𝐸 𝑆𝑠, 𝑘𝑟𝑚𝑠 and 𝑠𝑘 . The optimized GPR algorithm with these features as
inputs gives overall very good predictions of velocity profiles in all cases, with a maximum error
of 8%. Upon analyzing the iso-contours of 3D time-averaged velocity patterns, the formation and
growth of local mean separation zones downstream of roughness elements is shown to be dependent
on roughness rotation (𝐼𝑥) and its peak distribution (measured by 𝑅𝑙𝑥 or porosity), explaining the
dependence of 𝑈 on these features.

For the analysis and prediction of 𝑘 𝑠, the datasets of [4] is expanded to include cuboidal-type
roughness, multiscale roughness and a significant amount of pitty surfaces (𝑠𝑘 < 0). Additional
two-point features (i.e.
correlation lengths) are added to the list of inputs for the first time,
to incorporate spatial information such as separation of local peaks. This addition significantly
improves the NN algorithm by reducing both 𝐿1 and 𝐿2 errors by a half. A sweep of NN predictions
with different numbers of principal components (in the surface feature space) as inputs shows the
minimal set of important features for 𝑘 𝑠 prediction to be (𝑘𝑟𝑚𝑠, 𝑠𝑘 , 𝑘𝑢, 𝐸 𝑆𝑥, 𝐸 𝑆𝑥, 𝑅𝑙𝑥, 𝑅𝑙𝑧) and some

84

of their combinations. These features are overall consistent with those important ones for 𝑈 (𝑦)
prediction, due to similar flow physics involved (i.e. local separation zones) in determining mean
velocity and drag. 𝐼𝑥, 𝐼𝑧 are only shown to be important with addition of new principal components
(that contained remaining 20% of the variance of input features), indicating insufficient variation
of inclination angles in the available datasets. Addition of more datasets with varying inclination
angles while keeping the shape constant and repeating this for different shape of rough surfaces
will aid better prediction of 𝑘 𝑠.

To summarize, this work studied the response of mean and turbulent fields in non-equilibrium
wall turbulence dependent on roughness texture. The major contributors to the different mean-
momentum and Reynolds stress transports on rough surfaces as compared to a smooth wall are
identified. Qualitative estimation and quantitative prediction of their dependencies on roughness
texture are carried out, which also identifies a set of surface features that are important for dynamics
in accelerated turbulent flows. The identification of the important features will aid future generation
of additional rough surfaces to complement the datasets for better understanding and modeling
rough-wall flows. In addition, the mean-flow scaling relations and self-similarity discovered for
non-equilibrium accelerated turbulence provide significant insights for physics-based modeling of
rough-wall non-equilibrium flows in general.

5.2 Future work

First, identification of the important features in Chapter 4 will (1) aid in future generation of
regression models that can be integrated into computationally inexpensive RANS solvers, and (2)
offer guidance in generation of new rough surfaces to deepen the understanding of the roughness
texture effects and increase the accuracy and universality of roughness models. As these new
datasets are added, data classification techniques that can identify patterns in the datasets, similar
to clustering in 𝑈 (𝑦), can then be explored. This will survey the datasets for comprehensiveness
and may provide better strategies in splitting datasets into training-testing sets in ML algorithm
development. In addition, a more systematic analysis is required to understand individual effects of
the surface features, by synthesizing rough surfaces with variation of one feature at a time. This is
important for white box modeling using the identified important rough features and understanding
their influence on the wake field. A larger representation of multiscale roughness is also needed
(e.g. power-law roughness) to expand the roughness database. Future work based on the PCA
analyses may explore explicit expressions of the few orthogonal features based on the large number
of surface features used by the community, for effective feature reduction that will benefit roughness
modeling.

Next, future work is needed to expand the roughness models to non-equilibrium flows in
general, including accelerating or decelerating flows (in time or in space) over rough surfaces.
As an extension of this work, a recent study [3] compared non-equilibrium flat-plate turbulent

85

boundary layers under spatial acceleration and deceleration (including both attached and separating
flows). The mean-flow (𝑢𝜏, 𝑈, ⟨(cid:101)
𝑢 𝑗 ⟩, etc.) scaling on 𝑈𝑅 was found to apply also in the spatially
𝑢𝑖 (cid:101)
varying flows therein. For example, 𝑈 (𝑥, 𝑦)/𝑈𝑅 (𝑥) inside the RSL was found self-similar in the
attached flow with sufficiently high 𝑘 +, regardless of acceleration. This means that the current
𝑈 (𝑦) model developed for canonical channel flows may be applicable also to non-equilibrium
flows. However, a knowledge gap exists for non-equilibrium rough-wall flows under low 𝑘 +
(e.g. in APG/decelerating flows close to separation point), where the above mean-flow quantities
may depend also on 𝑘 + and pressure gradient/acceleration. Future investigations are needed on
understanding the dependencies of the mean flow and turbulence in flows with low 𝑘 + values.
To do so, simulation/experimental measurements are needed in temporally/spatially-varying wall-
bounded turbulent flows under systematically varying roughness geometry, 𝑘 + and acceleration.
This is a future endeavor needed to fully understand effects of roughness in non-equilibrium wall
turbulence.

86

BIBLIOGRAPHY

[1]

[2]

[3]

A. Scotti, “Direct numerical simulation of turbulent channel flows with boundary roughened
with virtual sandpaper,” Phys. Fluids, vol. 18, pp. 031701–1—4, 2006.

S. C. Mangavelli, J. Yuan, and G. J. Brereton, “Effects of surface roughness topography in
transient channel flows,” J. Turbul., vol. 22, pp. 434–460, 2021.

J. Yuan, S. C. Mangavelli, and G. Brereton, “Roughness sublayer in non-equilibrium wall-
bounded turbulence,” Bulletin of the American Physical Society, 2023.

[4] M. Aghaei Jouybari, J. Yuan, G. J. Brereton, and M. S. Murillo, “Data-driven prediction
of the equivalent sand-grain height in rough-wall turbulent flows,” J. Fluid Mech., vol. 912,
pp. A8–1—23, 2021.

[5]

[6]

[7]

[8]

R. Grenier, P. Oryschuk, and P. Léveillé, “Components surface deterioration and effect on
performance,” Tech. Rep. T032700 0323A, CEATI, 2007.

E. A. Murphy, J. M. Barros, M. P. Schultz, K. A. Flack, C. N. Steppe, and M. A. Rei-
denbach, “Roughness effects of diatomaceous slime fouling on turbulent boundary layer
hydrodynamics,” Biofouling, vol. 34, no. 9, pp. 976–988, 2018.

J. Nikuradse, “Strömungsgesetze in rauhen rohren,” VDI-Forsch., 1933.

C. F. Colebrook, “Turbulent flow in pipes with particular reference to the transition region
between smooth- and rough-pipe laws,” J. Inst. Civ. Eng., vol. 11, pp. 133–156, 1939.

[9] M. R. Raupach, R. A. Antonia, and S. Rajagopalan, “Rough-wall boundary layers,” Appl.

Mech. Rev., vol. 44, pp. 1–25, 1991.

[10]

J. Jiménez, “Turbulent flows over rough walls,” Annu. Rev. Fluid Mech., vol. 36, pp. 173–196,
2004.

[11] M. R. Raupach and A. S. Thom, “Turbulence in and above plant canopies,” Annu. Rev. Fluid

Mech., vol. 13, pp. 97–129, 1981.

[12]

J. Finnigan, “Turbulence in plant canopies,” Annu. Rev. Fluid Mech., vol. 32, pp. 519–571,
2000.

[13] H. S. Shafi and R. A. Antonia, “Anisotropy of the reynolds stresses in a turbulent boundary

layer on a rough wall,” Exp. Fluids, vol. 18, pp. 213–215, 1995.

[14] P.-Å. Krogstad, R. A. Antonia, and L. W. B. Browne, “Comparison between rough-and
smooth-wall turbulent boundary layers,” J. Fluid Mech., vol. 245, pp. 599–617, 1992.

[15] M. P. Schultz and K. A. Flack, “The rough-wall turbulent boundary layer from the hydrauli-
cally smooth to the fully rough regime,” J. Fluid Mech., vol. 580, pp. 381–405, 2007.

[16] G. J. Kunkel, J. J. Allen, and A. J. Smits, “Further support for Townsend’s Reynolds number
similarity hypothesis in high Reynolds number rough-wall pipe flow,” Phys. Fluids, vol. 19,
pp. 055109–1—6, 2007.

[17] K. A. Flack, M. P. Schultz, and J. S. Connelly, “Examination of a critical roughness height

for outer layer similarity,” Phys. Fluids, vol. 19, p. 095104, 2007.

[18] R. J. Volino, M. P. Schultz, and K. A. Flack, “Turbulence structure in rough-and smooth-wall

87

boundary layers,” J. Fluid Mech., vol. 592, pp. 263–293, 2007.

[19] K. A. Flack and M. P. Schultz, “Review of hydraulic roughness scales in the fully rough

regime,” J. Fluids Eng., vol. 132, pp. 041203–1—10, 2010.

[20] P.-Å. Krogstad and R. A. Antonia, “Surface roughness effects in turbulent boundary layers,”

Exp. Fluids, vol. 27, no. 5, pp. 450–460, 1999.

[21] R. J. Volino, M. P. Schultz, and K. A. Flack, “Turbulence structure in a boundary layer with

two-dimensional roughness,” J. Fluid Mech., vol. 635, pp. 75–101, 2009.

[22] P. Passalacqua, F. Porté-Agel, E. Foufoula-Georgiou, and C. Paola, “Application of dynamic
subgrid-scale concepts from large-eddy simulation to modeling landscape evolution,” Water
Resour. Res., vol. 42, pp. W06D11–1–11, 2006.

[23] R. Mejia-Alvarez and K. T. Christensen, “Wall-parallel stereo particle-image velocimetry
measurements in the roughness sublayer of turbulent flow overlying highly irregular rough-
ness,” Phys. Fluids, vol. 25, pp. 115109–1–24, 2013.

[24] C. Vanderwel and B. Ganapathisubramani, “Effects of spanwise spacing on large-scale
secondary flows in rough-wall turbulent boundary layers,” J. Fluid Mech., vol. 774, p. R2,
2015.

[25]

J. Yang and W. Anderson, “Numerical study of turbulent channel flow over surfaces with
variable spanwise heterogeneities: Topographically-driven secondary flows affect outer-layer
similarity of turbulent length scales.,” Flow, Turb. Combust., vol. 100, pp. 1–17, 2018.

[26] H. Bai, Kevin, N. Hutchins, and J. Monty, “Turbulence modifications in a turbulent boundary
layer over a rough wall with spanwise-alternating roughness strips,” Physics of Fluids, vol. 30,
no. 5, p. 055105, 2018.

[27] R. Mittal and G. Iaccarino, “Immersed boundary methods,” Annu. Rev. Fluid Mech., vol. 37,

no. 1, pp. 239–261, 2005.

[28] S. Leonardi, P. Orlandi, R. J. Smalley, L. Djenidi, and R. A. Antonia, “Direct numerical
simulations of turbulent channel flow with transverse square bars on one wall,” J. Fluid
Mech., vol. 491, pp. 229–238, 2003.

[29] P. Orlandi and S. Leonardi, “Direct numerical simulation of three-dimensional turbulent
rough channels: parameterization and flow physics,” J. Fluid Mech., vol. 606, pp. 399–415,
2008.

[30] M. Kanda, “Large-eddy simulations on the effects of surface geometry fo building arrays on
turbulent organized structures,” Bound.-Lay. Meteorol., vol. 118, pp. 151–168, 2006.

[31] S.-H. Lee and H. J. Sung, “Direct numerical simulation of the turbulent boundary layer over

a rod-roughened wall,” J. Fluid Mech., vol. 584, pp. 125–146, 2007.

[32]

[33]

J. H. Lee, H. J. Sung, and P.-Å. Krogstad, “Direct numerical simulation of the turbulent
boundary layer over a cube-roughened wall,” J. Fluid Mech., vol. 669, pp. 397–431, 2011.

I. P. Castro, A. Segalini, and P. H. Alfredsson, “Outer-layer turbulence intensities in smooth-
and rough-wall boundary layers,” Journal of Fluid Mechanics, vol. 727, pp. 119–131, 2013.

[34] M. Aghaei Jouybari, G. J. Brereton, and J. Yuan, “Turbulence structures over realistic and

88

synthetic wall roughness in open channel flow at Re𝜏 = 1000,” J. Turbul., vol. 20, pp. 723–
749, 2019.

[35] R. Ma, K. Alamé, and K. Mahesh, “Direct numerical simulation of turbulent channel flow

over random rough surfaces,” Journal of Fluid Mechanics, vol. 908, 2021.

[36] M. R. Raupach and R. H. Shaw, “Averaging procedures for flow within vegetation canopies,”

Bound.-Lay. Meteorol., vol. 22, pp. 79–90, 1982.

[37] E. Mignot, E. Bartheleemy, and D. Hurther, “Double-averaging analysis and local flow
characterization of near-bed turbulence in gravel-bed channel flows,” J. Fluid Mech., vol. 618,
pp. 279–303, 2009.

[38]

J. Yuan and U. Piomelli, “Roughness effects on the Reynolds stress budgets in near-wall
turbulence,” J. Fluid Mech., vol. 760, p. R1, 2014.

[39] O. Coceal, T. G. Thomas, I. P. Castro, and S. E. Belcher, “Mean flow and turbulence statistics
over groups of urban-like cubical obstacles,” Bound.-Lay. Meteorol., vol. 121, pp. 491–519,
2006.

[40]

J. Yuan and M. Aghaei-Jouybari, “Topographical effects of roughness on turbulence statistics
in roughness sublayer,” Phys. Rev. Fluids, vol. 3, p. 114603, 2018.

[41] E. Mignot, E. Barthélemy, and D. Hurther, “Turbulent kinetic energy budget in a gravel-bed

channel flow,” Acta Geophysica, vol. 56, pp. 601–613, 2008.

[42] R. Narasimha and K. R. Sreenivasan, “Relaminarization in highly accelerated turbulent

boundary layers,” J. Fluid Mech., vol. 61, pp. 417–447, 1973.

[43] C. Bourassa and F. O. Thomas, “An experimental investigation of a highly accelerated

turbulent boundary layer,” J. Fluid Mech., vol. 634, pp. 359–404, 2009.

[44] U. Piomelli and J. Yuan, “Numerical simulations of spatially developing, accelerating bound-

ary layers,” Phys. Fluids, vol. 25, pp. 101304–1—21, 2013.

[45] W. Schoppa and F. Hussain, “Coherent structure generation in near-wall turbulence,” J. Fluid

Mech., vol. 453, pp. 57–108, 2002.

[46]

J. Jiménez and A. Pinelli, “The autonomous cycle of near-wall turbulence,” J. Fluid Mech.,
vol. 389, pp. 335–359, 1999.

[47] R. J. Volino, “Non-equilibrium development in turbulent boundary layers with changing

pressure gradients,” J. Fluid Mech., vol. 897, pp. A2–1–48, 2020.

[48] Y. Na and P. Moin, “The structure of wall-pressure fluctuations in turbulent boundary layers
with adverse pressure gradient and separation,” J. Fluid Mech., vol. 377, pp. 347–373, 1998.

[49] M. Skote and D. S. Henningson, “Direct numerical simulation of a separated turbulent

boundary layer,” J. Fluid Mech., vol. 471, pp. 107–136, 2002.

[50] W. Cheng, D. I. Pullin, and R. Samtaney, “Large-eddy simulation of separation and reat-
tachment of a flat plate turbulent boundary layer,” J. Fluid Mech., vol. 785, pp. 78–108,
2015.

[51] W. Wu and U. Piomelli, “Effects of surface roughness on a separating turbulent boundary

89

layer,” J. Fluid Mech., vol. 841, pp. 552–580, 2018.

[52] R. Vinuesa, R. Örlü, and P. Schlatter, “Characterisation of backflow events over a wing

section,” J. Turbul., vol. 18, no. 2, pp. 170–185, 2017.

[53] Z. Harun, J. P. Monty, R. Mathis, and I. Marusic, “Pressure gradient effects on the large-scale
structure of turbulent boundary layers,” J. Fluid Mech., vol. 715, pp. 477–498, 2013.

[54] C. S. Vila, R. Vinuesa, S. Discetti, A. Ianiro, P. Schlatter, and R. Örlü, “Separating adverse-
pressure-gradient and reynolds-number effects in turbulent boundary layers,” Physical Re-
view Fluids, vol. 5, no. 6, p. 064609, 2020.

[55] S. He and J. D. Jackson, “A study of turbulence under conditions of transient flow in a pipe,”

J. Fluid Mech., vol. 408, pp. 1–38, 2000.

[56] S. He and M. Seddighi, “Turbulence in transient channel flow,” J. Fluid Mech., vol. 715,

pp. 60–102, 2013.

[57] R. B. Cal, B. Brzek, T. Johansson, and L. Castillo, “The rough favourable pressure gradient

turbulent boundary layer,” J. Fluid Mech., vol. 641, pp. 129–155, 2009.

[58]

J. Yuan and U. Piomelli, “Numerical simulation of a spatially developing accelerating
boundary layer over roughness,” J. Fluid Mech., vol. 780, pp. 192–214, 2015.

[59] S. He and M. Seddighi, “Transition of transient channel flow after a change in Reynolds

number,” J. Fluid Mech., vol. 764, pp. 395–427, 2015.

[60] C. D. Ghodke and S. V. Apte, “DNS study of particle-bed–turbulence interactions in an

oscillatory wall-bounded flow,” J. Fluid Mech., vol. 792, pp. 232–251, 2016.

[61] C. D. Ghodke and S. V. Apte, “A numerical investigation to study roughness effects in
oscillatory flows,” in Proceedings of the ASME 2017 Fluids Engineering Division Summer
Meeting, 2017.

[62] L. Djenidi, M. Kamruzzaman, and L. Dostal, “Effects of wall suction on a 2d rough wall

turbulent boundary layer,” Experiments in Fluids, vol. 60, no. 3, pp. 1–11, 2019.

[63] E. Napoli, V. Armenio, and M. De Marchis, “The effect of the slope of irregularly distributed
roughness elements on turbulent wall-bounded flows,” J. Fluid Mech., vol. 613, pp. 385–394,
2008.

[64]

J. Yuan and U. Piomelli, “Estimation and prediction of the roughness function on realistic
surfaces,” J. Turbul., vol. 15, pp. 350–365, 2014.

[65] L. Chan, M. MacDonald, D. Chung, N. Hutchins, and A. Ooi, “A systematic investigation of
roughness height and wavelength in turbulent pipe flow in the transitionally rough regime,”
J. Fluid Mech., vol. 771, pp. 743–777, 2015.

[66]

J. A. Van Rij, B. Belnap, and P. Ligrani, “Analysis and experiments on three-dimensional,
irregular surface roughness,” J. Fluids Eng., vol. 124, no. 3, pp. 671–677, 2002.

[67] M. Thakkar, A. Busse, and N. D. Sandham, “Surface correlations of hydrodynamic drag for

transitionally rough engineering surfaces,” J. Turbul., vol. 18, pp. 138–169, 2017.

[68] D. Chung, N. Hutchins, M. P. Schultz, and K. A. Flack, “Predicting the drag of rough

90

surfaces,” Annu. Rev. Fluid Mech., vol. 53, pp. 439–471, 2021.

[69] K. A. Flack and D. Chung, “Important parameters for a predictive model of 𝑘 𝑠 for zero

pressure gradient flows,” in AIAA SCITECH 2022 Forum, 2022.

[70] R. D. Sanhueza, I. Akkerman, and J. W. Peeters, “Machine learning for the prediction of
the local skin friction factors and nusselt numbers in turbulent flows past rough surfaces,”
International Journal of Heat and Fluid Flow, vol. 103, p. 109204, 2023.

[71] P. A. Durbin, “Reflections on roughness modelling in turbulent flow,” Journal of Turbulence,

vol. 24, no. 1-2, pp. 3–13, 2023.

[72] A. Busse and N. D. Sandham, “Parametric forcing approach to rough-wall turbulent channel

flow,” Journal of Fluid Mechanics, vol. 712, pp. 169–202, 2012.

[73] S. Wu, K. T. Christensen, and C. Pantano, “Modelling smooth-and transitionally rough-
wall turbulent channel flow by leveraging inner–outer interactions and principal component
analysis,” Journal of Fluid Mechanics, vol. 863, pp. 407–453, 2019.

[74] X. Yang and C. Meneveau, “Large eddy simulations and parameterisation of roughness
element orientation and flow direction effects in rough wall boundary layers,” Journal of
Turbulence, vol. 17, no. 11, pp. 1072–1085, 2016.

[75] G. J. Brereton, M. Aghaei Jouybari, and J. Yuan, “Towards modeling of turbulent flow over

surfaces of arbitrary roughness,” Phys. Fluids, vol. 33, pp. 065121–1—13, 2021.

[76]

J. Yuan, J. Nicolle, U. Piomelli, and A.-M. Giroux, “Modelling roughness and acceleration
effects with application to the flow in a hydraulic turbine,” IOP Conf. Series: Earth and
Environmental Science, vol. 22, p. 022016, 2014.

[77] R. Dutta, J. Nicolle, A.-M. Giroux, and U. Piomelli, “Evaluation of turbulence models on
roughened turbine blades,” in IOP Conf. Series: Earth and Environmental Science, vol. 49,
p. 062007, 2016.

[78] R. J. Volino, W. J. Devenport, and U. Piomelli, “Questions on the effects of roughness and
its analysis in non-equilibrium flows,” Journal of Turbulence, vol. 23, no. 8, pp. 454–466,
2022.

[79]

J. Yuan and U. Piomelli, “Numerical simulations of sink-flow boundary layers over rough
surfaces,” Phys. Fluids, vol. 26, pp. 015113–1—015113–28, 2014.

[80] A. Keating, U. Piomelli, K. Bremhorst, and S. Nešić, “Large-eddy simulation of heat transfer

downstream of a backward-facing step,” J. Turbul., vol. 5, p. N20, 2004.

[81] R. D. Moser and P. Moin, “The effects of curvature in wall-bounded turbulent flows,” J.

Fluid Mech., vol. 175, pp. 479–510, 1987.

[82] P. S. Jackson, “On the displacement height in the logarithmic velocity profile,” J. Fluid

Mech., vol. 111, pp. 15–25, 1981.

[83] S. B. Pope, Turbulent Flows. Cambridge, U.K.: Cambridge Univ. Press, 2000.

[84] H. M. Nagib and K. A. Chauhan, “Variations of von Kármán coefficient in canonical flows,”

Phys. Fluids, vol. 20, pp. 15180–1–10, 2008.

91

[85] V. Nikora, K. Koll, I. McEwan, S. McLean, and A. Dittrich, “Velocity distribution in the
roughness layer of rough-bed flows,” J. Hydr. Engng, vol. 130, pp. 1036–1042, 2004.

[86] D. Pokrajac, L. J. Campbell, V. Nikora, and I. Manes, C. adn McEwan, “Quadrant analysis
of persistent spatial velocity perturbations over square-bar roughness,” Exp. Fluids, vol. 42,
pp. 413–423, 2007.

[87]

J. Yuan and G. J. Brereton, “Structure of turbulent flows on heterogeneous irregular rough
surfaces,” in Proceedings of 11𝑡ℎ International Symposium of Turbulence and Shear Flow
Phenomena (TSFP-11), 2019.

[88] M. J. Lee, J. Kim, and P. Moin, “Structure of turbulence at high shear rate,” J. Fluid Mech.,

vol. 216, pp. 561–583, 1990.

[89] T. B. Nickels, “Inner scaling for wall-bounded flows subject to large pressure gradients,” J.

Fluid Mech., vol. 521, pp. 217–239, 2004.

[90] E. R. van Driest and C. B. Blumer, “Boundary layer transition-freestream turbulence and

pressure gradient effects,” AIAA J., vol. 1, pp. 1303–1306, 1963.

[91] S. Mangavelli and J. Yuan, “Effects of form-induced velocity in rough-wall turbulent channel

flows,” Journal of Turbulence, vol. 24, no. 1-2, pp. 14–35, 2023.

[92]

[93]

[94]

J. Fröhlich, C. P. Mellen, W. Rodi, L. Temmerman, and M. A. Leschziner, “Highly resolved
large-eddy simulation of separated flow in a channel with streamwise periodic constrictions,”
J. Fluid Mech., vol. 526, pp. 19–66, 2005.

J. H. Ferziger, M. Perić, and R. L. Street, Computational methods for fluid dynamics. Berlin:
springer, 2002.

J. C. R. Hunt, “A review of the theory of rapidly distorted turbulent flows and its applications,”
Fluid Dyn. Trans, vol. 9, pp. 121–152, 1978.

[95] K. M. Womack, R. J. Volino, C. Meneveau, and M. P. Schultz, “Turbulent boundary layer flow
over regularly and irregularly arranged truncated cone surfaces,” J. Fluid Mech., vol. 933,
p. A38, 2022.

[96] K. Bhaganagar, G. N. Coleman, and J. Kim, “Effect of roughness on pressure fluctuations in

a turbulent channel flow,” Phys. Fluids, vol. 19, pp. 028103–1—4, 2007.

[97] T. Meyers, J. B. Forest, and W. J. Devenport, “The wall-pressure spectrum of high-Reynolds-
number turbulent boundary-layer flows over rough surfaces,” J. Fluid Mech., no. 768,
pp. 261–293, 2015.

[98] Q. Yang and M. Wang, “Boundary-layer noise induced by arrays of roughness elements,” J.

Fluid Mech., vol. 727, pp. 282–317, 2013.

[99] O. Coceal, A. Dobre, T. G. Thomas, and S. E. Belcher, “Structure of turbulent flow over

regular arrays of cubical roughness,” J. Fluid Mech., vol. 589, pp. 375–409, 2007.

[100] X. Yang, J. Sadique, R. Mittal, and C. Meneveau, “Exponential roughness layer and analytical
model for turbulent boundary layer flow over rectangular-prism roughness elements,” J. Fluid
Mech., vol. 789, pp. 127–165, 2016.

[101] O. Coceal, T. G. Thomas, and S. E. Belcher, “Spatial variability of flow statistics within

92

regular building arrays,” Bound.-Lay. Meteorol., vol. 125, pp. 537–552, 2007.

[102] K. A. Flack, M. P. Schultz, J. M. Barros, and Y. C. Kim, “Skin-friction behavior in the

transitionally-rough regime,” Int. J. Heat Fluid Flow, vol. 61, pp. 21–30, 2016.

[103] P. Forooghi, A. Stroh, F. Magagnato, S. Jakirlic, and B. Frohnapfel, “Toward a universal

roughness correlation,” J. Fluids Eng., vol. 139, pp. 121201–1–12, 2017.

[104] Y. Kuwata and Y. Kawaguchi, “Direct numerical simulation of turbulence over systematically
varied irregular rough surfaces,” Journal of Fluid Mechanics, vol. 862, pp. 781–815, 2019.

[105] A. J. Musker, “Universal roughenss functions for naturally-occuring surfaces,” Trans. Can.

Soc. Mech. Eng., no. 1, pp. 1–6, 1980-1981.

[106] H. Ma, Y. Li, X. Yang, and L. Ye, “Data-driven prediction of the equivalent sand-grain

roughness,” Scientific Reports, vol. 13, no. 1, p. 19108, 2023.

[107] X. Yang, W. Zhang, J. Yuan, and R. F. Kunz, “In search of a universal rough wall model,” J.

Fluids Eng., vol. to be submitted, 2023.

[108] E. Florens, O. Eiff, and F. Moulin, “Defining the roughness sublayer and its turbulence

statistics,” Experiments in fluids, vol. 54, pp. 1–15, 2013.

[109] D. Chung, L. Chan, M. MacDonald, N. Hutchins, and A. Ooi, “A fast direct numerical

simulation method for characterising hydraulic roughness,” J. Fluid Mech., 2015.

[110] M. MacDonald, D. Chung, N. Hutchins, L. Chan, A. Ooi, and R. García-Mayoral, “The
minimal-span channel for rough-wall turbulent flows,” J. Fluid Mech., vol. 816, pp. 5–42,
2017.

[111] A. Busse and T. O. Jelly, “Influence of surface anisotropy on turbulent flow over irregular

roughness,” Flow, Turbulence and Combustion, vol. 104, pp. 331–354, 2020.

[112] W. Wu, U. Piomelli, and J. Yuan, “Turbulence statistics in rotating channel flows with rough

walls,” Int. J. Heat Fluid Flow, vol. 80, pp. 108467–1–13, 2019.

[113] J. Hong, J. Katz, and M. P. Schultz, “Near-wall turbulence statistics and flow structures over
three-dimensional roughness in a turbulent channel flow,” J. Fluid Mech., vol. 667, pp. 1–37,
2011.

[114] J. M. Barros, M. P. Schultz, and K. A. Flack, “Measurements of skin-friction of systematically
generated surface roughness,” International Journal of Heat and Fluid Flow, vol. 72, pp. 1–7,
2018.

[115] T. Medjnoun, C. Vanderwel, and B. Ganapathisubramani, “Characteristics of turbulent
boundary layers over smooth surfaces with spanwise heterogeneities,” Journal of Fluid
Mechanics, vol. 838, pp. 516–543, 2018.

[116] T. Medjnoun, E. Rodriguez-Lopez, M. Ferreira, T. Griffiths, J. Meyers, and B. Ganapathisub-
ramani, “Turbulent boundary-layer flow over regular multiscale roughness,” Journal of Fluid
Mechanics, vol. 917, p. A1, 2021.

93

A.1 Comparison between different ML methods

APPENDIX

Figure A.1 Comparison of different ML methods in 𝑈 (𝑦) prediction with different numbers of
roughness features used as inputs, based on (a) 𝐿∞ and (b) 𝐿2.

Figure A.1 compares the performances of different machine learning techniques on 𝑈 (𝑦)
prediction, based on the 𝐿∞ and 𝐿2 values (quantified as the maximum value among all cases
and averaged throughout all 𝑁 = 1300 runs). Different numbers of features used for training
are evaluated, similar to the procedure used for Figure 4.7(a). Two other regression methods:
Kernel Ridge (KR) and Decision Tree Regressor (DTR) are compared with GPR used herein.
Both GPR and KR (with third-degree polynomial) give overall similar predictions, with the error
norms decreasing with increasing number of features, reaching a minima at 4 or 5 features. This
consistency across different methods provides confidence in the findings in Chapter 4. However,
DTR gives significantly higher errors. GPR is finally selected as the method used in Chapter 4,
while another method may also be used.

A.2 Predictions with other sizes of test datasets

In the 𝑈 (𝑦) prediction in Chapter 4, 15% of the datasets are used as test data. To analyze
whether a change of this percentage affects the GPR prediction, the errors of the optimal GPR
predictions obtained with larger percentages (up to 50%) of data used as test data are shown in
Figure A.2(a). It is shown that 𝐿∞ and 𝐿2 error norms (the maximum values across all flow cases)
both increase sharply when more than 24% of the data are used as test data. This is due to the
lack of sufficient representation of different roughness types in the training set as a consequence of
limited data availability.

To provide more details, Figure A.2(b) shows that the percentages of the total datasets with
relatively high errors (i.e. 𝐿∞ > 0.1 and 𝐿2 > 0.005) increase with the size of the test datasets. It
is therefore concluded that, for the present available data, the applicability of the model is up to test
data being up to 24% of total datasets.

94

Figure A.2 Variation of percentage errors (a) and percentage of cases with errors larger than a
threshold (b) with the test-cases (as a percentage of all datasets).

95