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

LARGE-EDDY SIMULATION OF FREE-SURFACE TURBULENT CHANNELFLOW OVER SQUARE BARS

Razieh Jalalabadi
Civil, Environmental and Geomatic Engineering
University College London
London, WC1E 6BT, UK
r.jalalabadi@ucl.ac.uk

Thorsten Stoesser
Civil, Environmental and Geomatic Engineering
University College London
London, WC1E 6BT, UK
emailt.stoesser@ucl.ac.uk

ABSTRACTeraging method (Raupach & Shaw, 1982). In double-averaged
Navier-Stokes equations new terms known as dispersive shear
stresses emerge which represent the spatial variation of timeaveraged quantities. All spatial variations of first and second
order statistics are induced by near-bed flow heterogeneity andsecondary currents when it is present. In flow over rough bed
wake kinetic energy (WKE), calculated using the spatial variations of mean velocities, constitute one component of total
kinetic energy along with turbulent and mean kinetic energy(TKE and MKE respectively). The spatial variations of mean
flow variables are generally small compared to time-averaged
variables; however, for some rough surfaces they can contribute to flow dynamics substantially.
Considering the significant effects of geometrical surface
parameters on the distribution of first and second order statistics and flow structures (Jalalabadi et al., 2017), the normal
and shear stresses and kinetics energies for free-surface turbulent channel flow over spanwise-aligned square bars are investigated based on results of large-eddy simulation of flow overa transitional and a k-type roughness at constant Reynolds andFroude numbers. The distributions of mean stresses and kinetic energies present the effects of water surface deformations
induced by bed roughness on these flow variables. Quadrant
analysis is applied to reveal the contribution of different flow
movements on dispersive and Reynolds shear stresses. Shearstress and kinetic energy budgets are scrutinized to explore
their dependency on bed roughness and water surface deformations. The remainder of this paper includes numerical simulations, result and discussion and conclusions sections.

Large-eddy simulations of free surface turbulent channel
flow over spanwise-aligned square bars with transitional and ktype spacing are performed. Differences in the sign of spatial
and temporal velocity fluctuations entail sweep and ejection
events which contribute to the dispersive and Reynolds shear
stress only in flow over k-type roughness. The pressure-strain
and transport terms are the dominant components in the budgetof the dispersive and double-averaged Reynolds shear stresswhile the other terms are negligible. Pressure has the largest
contribution in the transport of these stresses. Investigating thebudgets of wake and double-averaged turbulent kinetic energy
reveals that the dominant terms contributing to the generationand fate of these energies are production, transport and convection. Water surface deformations have non-negligible effects
in the transport and convection of these stresses and energies.

INTRODUCTION
Turbulent open channel flows over rough surfaces are offundamental importance as the interaction of bed, bulk flow
and water surface influences the flow structure and water surface configuration. Better understanding of the hydrodynamics
of this flow type leads to significant improvements in analysing
free-surface flows in nature, like river flows, or in predicting
floods, respectively. In shallow flows the bed roughness alters the bulk flow structure and water surface directly and theroughness topography plays a crucial role.
Square bars perpendicular to the main flow direction, acanonical and rather simple geometry (Djenidi et al., 1999;
Jiminez, 2004) have been applied in numerous researches withthe aim to study principle effects of surface roughness on flow
characteristics. The main two types of this roughness are dtype and k-type based on the bar spacing. In d-type roughness a stable vortex is formed between the bars while in the
k-type roughness the flow reattaches to the bed at downstreamof bar. The transition from d-type to k-type roughness occurs
at λ /k = 4 − 5 where λ is the bar spacing and k is the barheight (Jiminez, 2004). The effects of bed roughness on the
averaged parameters can be represented using the double av-

NUMERICAL SIMULATIONS
The in-house LES code Hydro3D is employed which has
been validated for a large number of flows with similar complexity (Kara et al., 2015; McSherry et al., 2018). The codesolves the spatially filtered Navier-Stokes equations∇·u = 01

(1)

12th International Symposium on Turbulence and Shear Flow Phenomena (TSFP12)
Osaka, Japan (Online), July 19-22, 2022stress and turbulent kinetic energy along with spatial fluctuations of velocities, dispersive shear stress and wake kinetic
energy to investigate the effects of bars and water surface deformations on these flow variables. The dispersive shear stressis calculated as (Jelly & Busse, 2018)ũw̃ = (u − ⟨u⟩)(w − ⟨w⟩)

Figure 1.

where each parentheses on the right hand side represents thespatial fluctuation of a velocity component. Under the bar
crest all rms velocity fluctuations, Reynolds shear stress andturbulent kinetic energy are larger in flow over λ /k = 10.4
than those in λ /k = 5.2 case. Above the bar crest, only wallnormal rms velocity is slightly larger for λ /k = 10.4. The peakfor ⟨urms ⟩+ under the water surface in flow over λ /k = 10.4
represents the effects of the large water surface deformation as
this peak is not noticeable for flow over λ /k = 5.2. Similarvalues for the peak of ⟨vrms ⟩+ under water surface in both geometries reveal that only the water surface presence contributeto the generation of this peak regardless of its specific deformations. The wall-normal velocity experiences a weak local
increase under the water surface only in λ /k = 10.4 case. The
Reynolds shear stress is similar in both geometries except near
the water surface. The strong positive peak of ⟨u′ w′ ⟩+ under
the water surface in flow over λ /k = 10.4 is the effect of thestanding wave which is absent in flow over λ /k = 5.2 (Jalalabadi & Stoesser, 2022b). The strong peak of TKE under thewater surface is similar to the peak in ⟨urms ⟩+ in that area.
⟨eu⟩+ in flow over both bar spacings has similar values except
at the bar crest and under the water surface. At the bar crest⟨eu⟩+ is larger in flow over λ /k = 5.2 as, due to the formation of stable vortex between the bars, a significant momentum
transfer occurs at the bar crest in this geometry while in flow
over λ /k = 10.4 the momentum transfer is large between theflow everywhere under and above the bar crest. Under the water surface the standing wave causes large spatial variations inthe streamwise velocity in flow over larger bar spacing. ⟨ev⟩++e are negligible everywhere except close to the waterand ⟨w⟩e + experiences a weak local increase insurface where only ⟨w⟩
flow over λ /k = 10.4. Dispersive shear stress is similar forboth roughness spacing under the bar crest but larger for flow
over λ /k = 10.4 above that height. The peak under the wae + has the opposite sign of ⟨u′ w′ ⟩+ . Unlike
ter surface for ⟨euw⟩
other variable in figure 2, WKE is larger at all wall-normal locations in flow over λ /k = 10.4 than in flow over λ /k = 5.2
and reaches its maximum value at the free surface similar to⟨eu⟩+ .
The contribution of velocity fluctuations to the generationof Reynolds and dispersive shear stresses are shown in figure
3. These stresses are divided into four events, Q1 to Q4, basedon the sign of their fluctuating components (Wallace, 2016).
The second and forth quadrants, Q2+ and Q4+ , correspond tovertical flow away from the wall (w′ > 0) at low speed (u′ < 0)and towards the wall (w′ < 0) at high speed (u′ > 0) and known
as ejection and sweep events respectively. As shown in figure
3, Q2+ and Q4+ are the dominant events in both geometries
and shear stresses. For the Reynolds shear stress in both geometries, Q4+ is the strongest under the bar crest height while
Q2+ is the strongest above that. This is reversed for the dispersive shear stress in flow over λ /k = 10.4. This reveals thatthe sign of spatial fluctuations of streamwise and wall-normal
velocities are opposite of the sign of temporal fluctuations ofthese components. The negligible value of all quadrants of thedispersive shear stress in the flow over λ /k = 5.2 is due tosmall spatial fluctuations of velocities in this case. The stand-

Schematic of the computational domain.

∂u1+ u · ∇u = −∇p + ∇2 u − ∇ · τ
∂t
Re

(2)

where u = (u, v, w) is the velocity vector with the components
in the streamwise (x), spanwise (y) and wall-normal (z) directions, p is the pressure and τ is the subgrid scale stress tensor. A fractional-step method is applied on a staggered Cartesian grid where in the predictor step a second order RungeKutta scheme is applied to predict the velocities. In the corrector step the pressure Poisson equation is solved using a
multi-grid method to achieve a divergent flow field. The diffusive terms are computed using second order finite difference
scheme and the convective terms are computed using a fifthorder weighted, essentially non-oscillatory (WENO) scheme(Kara et al., 2015). The Wall-Adapting Local Eddy-viscosity(WALE) model is used to calculate eddy viscosity (Nicoud &Ducros, 1999). The free surface is captured using the Level
Set Method (Sethian & Smereka, 2003). In this method theinterface is tracked by solving a pure advection equation for alevel set signed distance function which is zero at the phase interface, negative in air and positive in water. Figure 1 presents
the schematic of the computational domain. The domain sizefor both geometries is 10.4k×5k×5k in the x, y and z directions respectively. Table 1 provides hydraulic parameters of
the simulations. Reynolds and Froude number are calculatedusing bulk velocity Ub and mean flow depth H. The superscript + represents the inner-scaled quantities calculated using
the friction velocity. Periodic boundary condition is applied in
the streamwise and spanwise directions. The governing equations are solved via parallel computing using Message Passing
Interface (MPI) and the present simulations were validated rigorously in Jalalabadi et al. (2021).

Table 1.

Hydraulic conditions and grid resolution.

λ /k

H/k

Ub

Re

Fr

∆x+ ,∆y+ ,∆z+

5.2

2.5

0.24

7.2×103

0.44

4.0,4.7,2.5

10.4

2.5

0.23

8.0×103

0.39

3.9,4.6,2.5(3)

RESULTS AND DISCUSSION
Figure 2 shows vertical profiles of spanwise averaged
root-mean-square (rms) velocity fluctuations, Reynolds shear2

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

Figure 2. Spatially-averaged (a) streamwise (b) spanwise (c) wall-normal RMS and dispersive velocity fluctuations, (d) Reynolds anddispersive shear stress and (e) turbulent and wake kinetic energy.

λ /k = 10.4 leads to the larger interaction of flow above andunder the bar crest while in this case PSC is large over thebar due to the larger form drag and gradient of dispersive velocities induced by the bars. The standing wave induces large
gradient of dispersive velocities (figure 2(a,c)) which results ina local peak in PSC under the water surface in figure 4(b). The
transport term T shown in figure 2(c,d) are generated mainly
by pressure transport and other terms are negligible. In flow
over λ /k = 5.2, T is large only under the bar crest and negligible above that implying that only bed roughness contribute to
the transport of dispersive shear stress. In larger bar spacing,on the other hand, both bars and water surface contribute to thetransport of dispersive shear stress.
∂ < uei uek >∂ < uei uek >∂ ⟨ui ⟩
= − uj+ − < uek uej >
∂t
∂xj
∂xj
C

e
e
∂
u
∂
u
∂ ⟨uk ⟩
i
k
′ u′
′ u′
g
− < uei uej >+ < ug>+<u>i j ∂xk j ∂x
∂xj
j
j
P




1
∂ uek
∂ uei
∂ uei ∂ uek
> +
(< pe
> + < pe>)+ −2ν <
∂xj ∂xj
ρ
∂ xi
∂ xk
PSC
ε

′ u′ ue′ u′ uge∂ ug∂u∂ < uei uej uek >
i
i j k
k j
+ −
−<
>−<
>
∂xj
∂xj
∂xj

D
E
1 ∂ peuek ∂ peuei− (+) + ν ∇2 uei uek (4)ρ ∂ xi∂ xk
TFigure 3. Quadrant components of Reynolds (black) and dispersive (red) shear stress for (a) /k = 5.2 and (b) /k = 10.4.

ing wave at the water surface contributes to the increase in
Q1+ to Q4+ under the water surface for λ /k = 10.4 while inflow over λ /k = 5.2 this increase is mild. The budgets of dispersive and double-averaged Reynolds shear stress and wakeand double-averaged turbulent kinetic energy will be shown to
discuss the dominant terms contributing to these second order
statistics and explore how they depends on roughness spacing
and water surface deformations.

The budgets of the double-averaged Reynolds shear stressare calculated using equation (5) for i =1 and k =3 and shown infigure 5. Similar to the dispersive shear stress, production (P),dissipation (ε) and convection (C) terms are negligible and notshown here and the main contribution to T is provided by pressure transport. Unlike the dispersive shear stress, large PSC
and T terms are not confined in the area below the bar crest inflow over λ /k = 5.2. In this geometry, large contribution of
PSC and T in the generation of ⟨u′ w′ ⟩+ are focused near the
bed and bar crest and slight undulations of water surface has
small contribution to that. In flow over λ /k = 10.4, the distribution of PSC is different than that in figure 4(b). The local
large value of PSC at the bar crest extends to the bed in thiscase and a strong positive contribution of PSC is localized under the surface wave. The local maxima of T under the surface
wave and at the downstream of bar crest are connected andextended to bed in figure 5(d) representing the interaction of

Dispersive and Double-Averaged Reynolds
Shear Stress Budget
Figure 4 presents the budgets of the dispersive shear stresscalculated using equation (4) for i =1 and k =3 for both bar
spacings. All terms in this section are normalized by H/uτ 3 .
In both cases the production (P), dissipation (ε) and convection (C) terms are negligible and not shown while the pressurestrain correlation (PSC) term is counterbalanced by the transport (T) term. In both geometries PSC is large near the bed
consistent with relatively large variations of ⟨eu⟩+ under thebar crest (figure 2(a)). On the vertical face of the bars, PSC
is larger in flow over λ /k = 5.2. This is attributed to thelarge spatial fluctuations of the mean pressure (Jalalabadi &Stoesser, 2022a) due to presence of stable vortex in the cavity
between the bars. Absence of this stable vortex in flow over3

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

Figure 4.

Figure 5.

Dominant terms of dispersive shear stress budget for (a,c) λ /k = 5.2 and (b,d) λ /k = 10.4.

Dominant terms of double-averaged Reynolds shear stress budget for (a,c) λ /k = 5.2 and (b,d) λ /k = 10.4.
"

#
"
#
′
′
∂ u′i ∂ u′k∂u1′ ∂ uki+ −2ν <
> +
(< p
+ p′
>)
∂xj ∂xj
ρ
∂ xi
∂ xk
ε
PSC

′
′
g
∂ < u′j u′j u′k >∂ uej ui uk
+ −
−<
>
∂xj
∂xj

1 ∂ < p′ u′k > ∂ < p′ u′i >) + ν < ∇2 u′i u′k > (5)
− (
+
ρ
∂ xi
∂ xk
T

water surface, bar and bed in the transport of double-averaged
Reynolds shear stress. To sum, as the large changes in spatial fluctuations of mean pressure and velocity occur under thebar crest, there is clear distinction in the distribution of dome + under and above bar crest in both cases
inant terms of ⟨euw⟩while in flow over larger bar spacing the effects of water surface leads to larger transport of dispersive shear stress in thebulk flow. The dominant term in budget of ⟨u′ w′ ⟩+ are strong
near the bed and bar crest in flow over λ /k = 5.2 while inlarger bar spacing the water surface contributes to the increaseof PSC and T and the larger interaction of bed, bar and water
surface enhances pressure transport in bulk flow in this case.

Wake and Double-Averaged Turbulent Kinetic
Energy Budget
The dominant terms contributing to the generation of
wake and double-averaged turbulent kinetic energy (1/2<
uei uei > and 1/2< u′i u′i >) are production (P), transport (T) andconvection (C) terms (equations 6 and 7). Figure 6(a,c,e)shows that the main contributors to the generation of WKE
are negligible everywhere except under the bar in flow over
λ /k = 5.2 which is consistent with figure 2(e). The large WKE
for flow over λ /k = 10.4 is supplied by the convection term so"#∂ < u′i u′k >∂ < u′i u′k >∂ ⟨ui ⟩
= − uj+ − < u′k u′j >
∂t
∂xj
∂xj
C

ei
∂ ⟨uk ⟩′ u′ ∂ uek > − < ug
′ u′ ∂ u− < u′i u′j >− < ug>i j ∂xk j ∂x
∂xj
j
j
P

4

12th International Symposium on Turbulence and Shear Flow Phenomena (TSFP12)
Osaka, Japan (Online), July 19-22, 2022
WKE is mainly generated by the interaction of water surface,bar and bed which leads to the larger spatial fluctuations offlow variables. The dominant term contributing to transportterm in both geometries is the first term in T in equation 6 andis named here as dispersive transport. It is similar to turbulent
transport (first term in T in equation 7) and large over the bar
crest and under the water surface in flow over k-type roughness.
∂ < uei uei >1∂ < uei uei >= − uj
∂t
2
∂xj

C
∂ ⟨ui ⟩∂ uei′′g
+ − < uei uej >+ < ui u j
>
∂xj
∂xj
P



∂ uei ∂ uei
1 ∂ < uei uei uej >+ −ν <> + −∂xj ∂xj
2
∂xj
ε
′
′
g
E
∂ ui u j uei
1 ∂ peuei
1 D
−<
>− (
) + ν ∇2 uei uei∂xjρ ∂ xi2T

pressure-strain correlation and transport terms. The large contribution of the dominant terms to dispersive shear stress arevisible only under the bar crest in flow over transitional roughness while the standing wave at the water surface over larger
bar spacing enhances the contribution of these terms above thebar crest too. Pressure transport has the largest contribution
to the generation of double-averaged Reynolds shear stress inboth cases and it is affected by both bars and water surface.
The dominant terms generating wake and double-averaged turbulent kinetic energy are production, transport and convection
terms. The production term is mainly generated at the bar crest
height while the convection term is large almost everywhere in
the bulk flow.

ACKNOWLEDGEMENTS(6)

Financial support was provided by the EPSRC/UK project“Rapid Monitoring of River Hydrodynamics and Morphology
using Acoustic Holography”, Grant No. EP/R022135/1.

"#∂ < u′i u′i >∂ < u′i u′i >1= − uj
∂t
2
∂xj

C
e
⟨u
⟩
∂
∂
u
i
i
′ u′
+ − < u′i u′j >− < ugi j ∂x >
∂xj
j
P
#
"

′ u′ u′ >′′∂<u∂ ui ∂ ui
1
i i j
> + −
+ −ν <∂xj ∂xj
2
∂xj
ε

′
′
′
′
g
e
∂
u
u
u
1
1 ∂ < p ui >
j i i
−<
(7)
>− (
) + ν < ∇2 u′i u′i >
2 ∂xj
ρ
∂xj
T

REFERENCES
Djenidi, L., Elavarasan, R. & Antonia, R. A. 1999 The turbulent boundary layer over transverse square cavities. Journalof Fluid Mechanics 395, 271–294.
Jalalabadi, R., Hwang, J., Nadeem, M., Yoon, M. & J., Sung H.
2017 Turbulent boundary layer over a divergent convergent
superhydrophobic surface. Physics of Fluids 29, 085112.
Jalalabadi, R. & Stoesser, T 2022a Large eddy simulation offree surface flow over square bars in laminar, transitional
and turbulent flows. 39th IAHR World Congress .
Jalalabadi, R. & Stoesser, T. 2022b Reynolds and dispersive shear stress in free-surface turbulent channel flow oversquare bars. Physical Review E 105, 035102.
Jalalabadi, R., Stoesser, T., Ouro, P., Luo, Q. & Xie, Z. 2021
Level set methods for fluid interfaces. Journal of HydroEnvironment Research 36, 64.
Jelly, T. O. & Busse, A. 2018 Reynolds and dispersive shear
stress contributions above highly skewed roughness. Journal of Fluid Mechanics 852, 710.
Jiminez, J 2004 The turbulent boundary layer over transverse
square cavities. Annual Review of Fluid Mechanics 36(4),173.
Kara, S., Kara, M. C., Stoesser, T. & Sturm, T. W. 2015
Free-surface versus rigid-lid LES computations for bridgeabutment flow. Journal of Hydraulic Engineering 141(9),04015019.
McSherry, R., Chua, K., Stoesser, T. & Mulahasan, S. 2018
Free surface flow over square bars at intermediate relative
submergence. Journal of Hydraulic Research 56(6), 825.
Nicoud, F. & Ducros, F. 1999 Subgrid-scale stress modellingbased on the square of the velocity gradient tensor. Flow
Turbulence Combustion 62(3), 183.
Raupach, M.R. & Shaw, R.H. 1982 Averaging procedures forflow within vegetation canopies. Boundary-Layer Meteorology 22(1), 79–90.
Sethian, J. A. & Smereka, P. 2003 Level set methods for fluid
interfaces. Annual Review of Fluid Mechanics 35, 146.
Wallace, J. M. 2016 Quadrant analysis in turbulence research:history and evolution. Annual Review of Fluid Mechanics
131, 131.

The main source for production of the double-averaged turbulent kinetic energy are the bars (figure 7(a,b)) while the large
deformation of the water surface contribute to P in flow overlarger bar spacing too. The main contribution to the transportof 1/2< u′i u′i > is turbulent transport and as seen in figure 7(d)standing wave over the water surface has a large contribution
to its generation. This is consistent with figure 2(a-c). The
large interaction of bed and water surface leads to the increasein the convection of double-averaged kinetic energy in figure7(d). There is a large interaction between bar crest and standing wave in terms of negative convection and a large interaction between the reattachement area over the bed and upstreamof the standing wave in terms of positive convection.

CONCLUSIONS
The results of large-eddy simulations of open channel
flow over spanwise aligned square bars show the effects ofroughness (bars) and water surface deformation on first andsecond order statistics. The large dispersive shear stress andwake kinetic energy in flow over larger bar spacing represents large interaction of bed and water surface in this case.
Reynolds shear stress and turbulent kinetic energy are larger
above the bar crest in flow over smaller bar spacing due to thepresence of stable vortex between bars which skims the flow
and limit the interaction of the flow under and above the bar
crest. Quadrant analysis of Reynolds shear stress reveals thedownward motion of fluid under the bar crest and its upward
motion above this height. The non-negligible terms in dispersive and double-averaged Reynolds shear stress budgets are512th International Symposium on Turbulence and Shear Flow Phenomena (TSFP12)
Osaka, Japan (Online), July 19-22, 2022

Figure 6.

Figure 7.

Dominant terms of wake kinetic energy budget for (a,c,e) λ /k = 5.2 and (b,d,f) λ /k = 10.4.

Dominant terms of turbulent kinetic energy budget for (a,c,e) λ /k = 5.2 and (b,d,f) λ /k = 10.4.

6