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

SIMULATION OF THE TURBULENT FLOW OVER A CANOPY OF HIGHLYFLEXIBLE BLADES

Bastian Löhrer
Institute of Fluid Mechanics
Technische Universität Dresden
01062 Dresden, Germany
bastian.loehrer@tu-dresden.deLéo Guiot de la RochèreLMFA, CNRS – UMR 5509, Université de Lyon, École Centrale de Lyon, Université Lyon 1, INSA Lyon
Écully, France
leo.guiotdelarochere@engees.unistra.fr

Delphine Doppler
LMFA, CNRS – UMR 5509, Université de Lyon, École Centrale de Lyon, Université Lyon 1, INSA Lyon
Écully, Francedelphine.doppler@univ-lyon1.fr

Sara Puijalon
Univ Lyon, Université Claude Bernard Lyon 1, CNRS, ENTPE, UMR 5023 LEHNA
69622, Villeurbanne, France
sara.puijalon@univ-lyon1.fr

Jochen Fröhlich
Institute of Fluid Mechanics
Technische Universität Dresden
01062 Dresden, Germany
jochen.froehlich@tu-dresden.deABSTRACT

from controlled model experiments and simulations alike, employing relatively rigid elements to facilitate optical accessibility and to reduce the complexity of simulations. Long andhighly flexible elements pose challenges to experiments andsimulations, and are addressed here.

In the present contribution, a recently developed own numerical method is used, capable of representing individual
flexible vegetation elements as they reconfigure and closelyinteract with the surrounding flow. The method is applied toa model configuration composed of very long flexible rectangular blades whose length and high flexibility translates into a
Cauchy number around 25,000. This very detailed, turbulenceresolving and plant-resolving simulation at high Cauchy number is the first of its kind, enabled by a composition of efficient
numerical modeling techniques.

PHYSICAL MODEL
Subject of the presented work is the fully developed flow
with bulk velocity Ub in a laterally confined open channel ofwidth Lz . Flexible blades of length L = Lz are attached perpendicular to the bottom wall at staggered locations, such thatthey form a dense canopy while also avoiding any channelization of the flow in between the stems. The overall problem
can be characterized by the nominal Cauchy number Ca =ρ Ub2 W L3 /(2 EI), where I = T 3W /12 is the second moment ofarea, the buoyancy number B = (ρ − ρs ) gW T L3 /(EI), and byReynolds numbers ReH = Ub H/ν and ReW = UbW /ν, basedon the channel height and the width of the blades, respectively.

INTRODUCTION
Aquatic canopies formed by flexible submerged plants arean essential element in river hydraulics through the flow resistance they provide, their resilience from bed erosion, andtheir impact on the transport of nutrients and pollutants (Nepf,2012). Knowledge about these processes is frequently drawn

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

2.0 1.5 1.0 0.5 0.0view from front/above with center sliceyz

u/Ubx

y/H

top viewz

Figure 1. Visualization of an instantaneous solution, showing the streamwise velocity component u in the center plane z = 0.5Lz , andthe instantaneous shape of the canopy with each blade colored according to the vertical position y.

COMPUTATIONAL MODEL

The corresponding values for the present case are listed in Table 1, of which the huge value of Ca and the relatively large
value of B reflect the large flexibility of the blades, i. e. nominal drag and lift forces outweighing elastic forces.

The simulation was configured according to an experiment presently conducted at the university of Lyon. The
computational domain spans a length of 18.5 H ≈ 11.5 Lz in
streamwise direction, where H is the water depth and Lz is thespanwise extent. Periodic boundary conditions were appliedin x-direction, a free-slip condition was imposed at y = H, andno-slip conditions were enforced at the bottom (y = 0) and theside walls. With a spatial resolution of 24 fluid grid cells perblade width, ∆x ≈ W /24, this translates into a numerical domain containing 460.8 million cells. The flow was driven by aspatially constant volume force which was adjusted in time tomaintain the desired flow rate. Additional numerical parameters are provided in Table 2.

NUMERICAL METHOD
The problem was analyzed on the basis of numerical simulations using the in-house code PRIME developed at TU
Dresden (Kempe & Fröhlich, 2012; Tschisgale & Fröhlich,2020). This method solves the Navier-Stokes equations forincompressible flows employing a second-order finite volume
approach on a staggered Cartesian grid for the spatial discretization and a second-order scheme for time integration. A
Large Eddy Simulation (LES) approach with the Smagorinsky
model was employed with damping close to the blades. In thesimulation reported below, the horizontal average of the ratioof turbulent to molecular viscosity remains below 0.25. The
motion of the flexible structures was computed with a finitedifference approach according to Lang & Arnold (2012), using a geometrically exact Cosserat rod model (Antman, 1995)
which is suitable for large deflections. While accounting forbending, torsional, and shearing motion, the model is built
upon the assumption of rigid cross sections which prevents
warping. Since warping is, however, a central mechanismof deformation its effect was modeled by reducing torsionaland shearing stiffnesses as proposed by Freund & Karakoç
(2016). Collisions were treated based on the model proposedby Tschisgale et al. (2019) with several technical enhancements. Coupling between fluid and structures was establishedwith a specialized Immersed-Boundary Method (IBM) applicable for small and even vanishing thickness of the blades, imposing a no-slip condition on the surfaces of the rods (Tschisgale & Fröhlich, 2020).

Table 2.

Numerical parameters.

parameter

value

description

Lx /L

11.5

channel length per blade length

Nx

4608

cells in streamwise direction

Ny

250

cells in vertical direction

Nz400

cells in spanwise direction

CS0.15

Smagorinsky constant

Ns528

number of blades

Ne108

number of elements per blade

RESULTSTable 1.

The snapshot shown in Figure 1 above gives a first impression of the flow. Due to their remarkable flexibility andtheir buoyancy, the blades are oriented almost horizontally
along most of their length, forming a visually dense canopy
hull. This is underpinned by the roughness density λ , defined as the ratio of the mean frontal area of a plant divided
by its share of the the bed area (Wooding et al., 1973), listed
in Table 1. A value of λ well above 0.1 indicates a dense
canopy (Nepf, 2012). The observed reconfigured canopy
height L∗ ≈ 0.24H places the case in the regime of shallow
submergence, H > L∗ > H/5 where canopy scale turbulence
dominates and boundary layer turbulence is usually not de-Dimensionless numbers.

numbervalue

description

Ca

25,000Nominal Cauchy number

B

145

Buoyancy number

ReW

1922

Reynolds number

ReH

20,064

Reynolds number

λ0.41

Roughness density

2

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

amplitude / mm

20Ub /H

Ub /LzUb /Lx

15Ub /H

251510105510−210−1

100f / Hz

10−210−1100f / Hz

Figure 3. Amplitude spectra for the position of blade tips,
i.e. at l = L. (a) Spanwise component; (b) vertical component. Spectra were computed per blade, then averaged acrossthe blades to obtain the solid black line. The gray areas span
the range of all and of 50 % of the spectra, respectively. Frequencies associated with the channel dimensions are indicatedby the dotted lines.

are shown in Figure 5a, manifesting an inclined, elliptical
pattern, elongated in streamwise direction. Vertical velocity
fluctuations are only correlated over a significantly shorter distance (Figure 5b) resulting in the well-known spherical patternalso observed by Shaw et al. (1995).

CONCLUSIONS
Scale-resolving simulations of the flow over highly flexible, resolved vegetation elements, give access to a broad rangeof data, including field quantities related to the flow, and related to the motion of the blades. Despite being characterized
by a pronounced secondary flow due to the fixed side walls,
statistics are relatively constant in a core region of the channel geometry. Yet, future investigations are to feature wider
domains with spanwise periodicity. Although the blades form
a dense canopy surface which undulates under the flow, noneof these apparent frequencies are present in spectra of blade
tip motion. This indicates, that pointwise blade motion statistics provide only a limited window into the motion of thecanopy. It may thus be more insightful to inspect the continuous canopy hull formed by the union of all blades.

⟨u′ (xx,t) u′ (xx + r ,t)⟩x,t⟨u′ 2 (xx,t)⟩x,t ⟨u′ 2 (xx + r ,t)⟩x,t(a)

(b)

20

Ub /Lzρu′ u′ (y, z; r ) = q(a)Ub /Lxveloped (Nepf, 2012). The canopy hull is occasionally broken open by high-speed streaks, extending all the way to thechannel bed and making the blades form ridges to either side.
The blades, hence, undergo a pronounced meandering motion in spanwise direction, while also undulating vertically.
Their mean height, computed as the maximum elevation of
their time-averaged shape (i.e. the height of their tips), is indicated in Figures 4a-b,e-g (broken line). Away from the sidewalls, the resulting profile coincides with an iso-contour of thestreamwise mean-velocity component ⟨u⟩, and with the line ofmaximum ⟨u′ u′ ⟩ in Figure 4e. The averaged streamwise velocity is maximum in the upper quarter of the channel between
the two counter-rotating swirls in the secondary flow below thefree surface, recognizable in Figure 4b. The symmetry-plane
profile of mean streamwise velocity in figure 4c has a narrow
region with ⟨u⟩ ≈ const. within the canopy and two inflections
points related to the drag of the blades. In contrast to situations
with stiffer blades (Nepf, 2012), the second inflection point is
not located at the height of the canopy edge, but instead occurs
considerably lower.
The dynamic behavior of the blades is analyzed by studying the motion of their centerlines. The points along theskeleton line of a blade are c s (l,t), depending on time, andthe arc length distance l from the respective fixed end of theblade, c s0 = c s (0,t) = const. Here, s is the index of therespective blade. Defining the geometries relative to theirroots, (xs , ys , zs )T = c s (l,t)−ccs0 , Figure 2a shows the average
shapes of the blade centerlines, while the root-mean-square deviations from this mean geometry are visualized in Figures 2bd. Here, the pronounced spanwise fluctuations are striking.
They reflect the secondary flow structures in Figure 4b. In afirst attempt to identify frequencies associated with blade motion, spectra were computed for the spanwise and vertical coordinates of the tips of blades (Figure 3). They exhibit a broad
range of frequencies without any remarkable peak, other thanfor a frequency linked to the flow-through time.
Autocorrelations of streamwise velocity fluctuations⟨x ⟩/Lz

0.5Ll /L1.0(d)l /L

0.00.5l /L0.5p⟨z ′ z ′ ⟩ / L

Antman, S. S. 1995 Nonlinear Problems of Elasticity, Applied
Mathematical Sciences, vol. 107. New York, NY: Springer
New York.
Freund, J. & Karakoç, A. 2016 Warping displacement of Timoshenko beam model. International Journal of Solids andStructures 92–93, 9–16.
Kempe, T. & Fröhlich, J. 2012 An improved immersed boundary method with direct forcing for the simulation of particle laden flows. Journal of Computational Physics 231 (9),3663–3684.
Lang, H. & Arnold, M. 2012 Numerical aspects in the dynamic
simulation of geometrically exact rods. Applied Numerical
Mathematics 62 (10), 1411–1427.

102

0.5
1.0

z

0.5

500.0

REFERENCES

Lz

0.0

10

1.0z/

2.5
0.0
0.0

102

5.0p
⟨y ′ y ′ ⟩ / L(c)0.5This work was funded by the French-German ANR-DFG
project ESCaFlex (ANR-16-CE92-0020, DFG grant 634058).
Computation time was provided by ZIH Dresden.

102

0.0

0.5
1.0z/0.500.0

ACKNOWLEDGEMENTS

Lz

0.01z/0.00.02⟨y ⟩ / L0.1p⟨x ′ x ′ ⟩ / L(b)

Figure 2. One-point statistics of the geometry of the bladeskeleton lines depending on the z-location of their fixed roots.
(a) Mean geometry of the blades; (b-d) root-mean-square values of their deviation from the mean shapes, showing x, y andz-components separately.

3

12th International Symposium on Turbulence and Shear Flow Phenomena (TSFP12)
Osaka, Japan (Online), July 19-22, 2022(a)(b) ⟨v⟩ /U , ⟨w⟩ /Ut,xt,xbb⟨u⟩t,x / Ub

1.00

(c)

1.5

y/H

0.75

⟨u′ u′ ⟩t,x / Ub2

1.0

0.10

0.5

0.05

⟨v′ v′ ⟩t,x / Ub2⟨w′ w′ ⟩t,x / Ub2

0.50
0.25
0.00

0.0
0.0

0.5

1.0

1.5(e)⟨u′ u′ ⟩t,x / Ub2

1.00

0.5

1.00

1.5

0.75

15

10.0

0.1

⟨u⟩t,x / Ubz/H×10−2−⟨u′ v′ ⟩t,x / Ub2

0.00
0.0z/Hy/H(d)0.15(f ) ⟨w′ w′ ⟩2t,x / Ub×10−26(g) −⟨u′ v′ ⟩2t,x / Ub×10−2

3104

25210.50
0.25
0.00

00.0

0.5

1.0

1.5

00.0z/H

0.5

1.0

1.5

00.0z/H

0.5

1.0

1.5z/Hrz = 0

0.5
0.0

ry = 0rx = 0

0.5
0.0
−7.5

0.0 0.5

−5.0

−2.5

ry /H(b)0.0

2.5

5.0rz = 0

0.0
−0.5

0.5
0.0

rx = 00.5
0.0
0.0 0.5

0.5

7.5

ry = 0ry /Hρu′ u′

rx /H

−7.5

−5.0

−2.50.02.5

5.0

7.5

|rz |/H ry /H(a)|rz |/H ry /HFigure 4. One-point statistics of the fluid motion in the channel cross section and profiles located at the vertical center line z = 0.5Lz .
(a) mean streamwise velocity; (b) mean secondary flow with streamlines and magnitude of secondary velocity components; (c) mean
streamwise velocity in center plane; (d) Reynolds stresses in center plane; (e-g) Reynolds stresses, as indicated. In the contour plots,
the broken lines represent the average canopy height corresponding to the vertical position of the time-averaged blade geometries. ⟨·⟩t,xdenotes averaging in time, followed by averaging in the periodic x-direction, while u′ = u − ⟨u⟩t is the deviation from the local time
average ⟨·⟩t .

ρv′ v′
0.5
0.0−0.5rx /H

Figure 5. Two-point correlations of velocity fluctuations, computed for the center of the channel cross section, i.e. y = 0.5H, z = 0.5Lz .
(a) streamwise velocity component; (b) vertical velocity component.

Nepf, H. M. 2012 Flow and Transport in Regions with AquaticVegetation. Annual Review of Fluid Mechanics 44, 123–142.
Shaw, R. H., Brunet, Y., Finnigan, J. J. & Raupach, M. R.
1995 A wind tunnel study of air flow in waving wheat:Two-point velocity statistics. Boundary-Layer Meteorology
76 (4), 349–376.
Tschisgale, S. & Fröhlich, J. 2020 An immersed boundary
method for the fluid-structure interaction of slender flex-ible structures in viscous fluid. Journal of Computational
Physics 423, 109801.
Tschisgale, S., Thiry, L. & Fröhlich, J. 2019 A constraintbased collision model for Cosserat rods. Archive of Applied
Mechanics 89 (2), 167–193.
Wooding, R. A., Bradley, E. F. & Marshall, J. K. 1973 Drag
due to regular arrays of roughness elements of varying geometry. Boundary-Layer Meteorology 5 (3), 285–308.

4