Numerical simulations of turbulent flows within an infinite array of randomly placed cylinders

Journal of Fluids and Structures 80 (2018) 245–261

Contents lists available at ScienceDirect

Journal of Fluids and Structures journal homepage: www.elsevier.com/locate/jfs

Numerical simulations of turbulent flows within an infinite array of randomly placed cylinders Ana M. Ricardo a, *, Dimokratis G.E. Grigoriadis b , Rui M.L. Ferreira a a b

CERIS, Instituto Superior Técnico, Universidade de Lisboa, Portugal UCY-CompSci, Department of Mechanical and Manufacturing Engineering, University of Cyprus, Cyprus

article

info

Article history: Received 11 September 2017 Received in revised form 14 February 2018 Accepted 9 April 2018

Keywords: Infinite array Random distribution LES Domain size Drag-wake controlled stratum

a b s t r a c t We address the task of modelling numerically turbulent flows within random arrays of circular cylinders, relevant for several industrial and environmental applications. Numerical simulations are employed to model infinite domains of randomly placed emergent and rigid cylinders validated by a laboratory database acquired by a PIV system. The main goals are: (i) to discuss the effect of the numerical domain size and the grid resolution on the first and second order moments; and (ii) to characterise the spatial distribution of mean flow and turbulence variables in the drag-wake controlled stratum. Five domains of different sizes (9–44 cylinders) and four grid resolutions were independently tested. The results show that the time-averaged velocity field and the Reynolds stress tensor were not significantly affected by the size of the tested numerical domains. The analysis of the grid resolution influence shows how the results improve with mesh refinement, while none of the tested meshes produces un-physical results. The present work provides guidance on the acceptable compromises, in terms of mesh resolution and domain size, when predicting, with eddy resolving computational fluid mechanics tools, first and second-order moments of turbulent flows within infinite arrays or randomly placed cylinders. © 2018 Elsevier Ltd. All rights reserved.

1. Introduction Vegetation plays a key role in river systems as it affects flow resistance and fluxes of sediments, nutrients and contaminants and it provides a large range of ecosystem services, justifying the research growth on hydrodynamics of wetland vegetation (Nepf, 2012; Aberle and Järvelä, 2013; Marjoribanks et al., 2014; Vargas-Luna et al., 2015; Ricardo et al., 2016a; Etminan et al., 2017). Rigid and emergent vegetation is often simulated, in idealised conditions, by circular cylinders, extending the research topic to a classical problem of fluid mechanics: the flow within arrays of cylinders, which features a large range of applications. Several studies with two (e.g. Sumner, 2010 and the references herein), three (e.g. Zhang and Zhou, 2001) or four (e.g. Lam and Lo, 1992; Wang et al., 2013) cylinders have produced important advances on the understanding of vortex formation, shedding and interaction process. Sewatkar et al. (2012) examined the minimum number of cylinders required to simulate the flow around a large number of in-line cylinders, by varying the number of in-line cylinders from four to ten. Cases with many cylinders have mainly considered squared or staggered configurations (Kim and Stoesser, 2011; Nicolle and Eames, 2011; Chang and Constantinescu, 2015; Yan et al., 2017a). Configurations with a random distribution of cylinders, which constitute the typical distribution for riparian vegetation, are less studied (Ghisalberti and Nepf, 2005; Tanino and Nepf, 2008, 2009; Ricardo et al., 2014b, 2016b).

*

Corresponding author. E-mail address: [email protected] (A.M. Ricardo).

https://doi.org/10.1016/j.jfluidstructs.2018.04.004 0889-9746/© 2018 Elsevier Ltd. All rights reserved.

246

A.M. Ricardo et al. / Journal of Fluids and Structures 80 (2018) 245–261

In complex geometries, as in the case of arrays of cylinders, laboratory or field measurements of 3D velocity fields typically face constraints limiting the amount of data acquired. It is often required to remove several cylinders in order to provide access to velocity probes or to allow flow visualisation. Nowadays, the advances in computational power allow high-resolution computational fluid dynamics (CFD) simulations providing detailed turbulent flow field information within realistic boundary conditions. The main methods employed in turbulent hydraulic flow computations are: (i) methods solving the Reynolds-averaged Navier–Stokes (RANS) equations, (ii) large-eddy simulations (LES) and (iii) direct numerical simulations (DNS). Rodi (2017) presents a review of the various methods for treating and simulating turbulence and its effects in hydraulic flows, covering the main features of each method and presenting application examples. DNS solves the full Navier–Stokes equations for all scales of fluid motion, therefore it might become computationally excessively demanding for large Reynolds numbers. In the RANS approach, turbulent fluctuations are averaged out and the effect of turbulent motions are accounted for, statistically, by a turbulence model. In practice, RANS is largely the most employed method, but it is of limited application when handling complex phenomena or complex geometries, particularly for flows with relevant turbulent transport by large-scale structures and strong anisotropy (Rodi, 2017; Stoesser, 2014). LES provides a compromise of computational cost and detail in flow description. It resolves large-scale unsteadiness and asymmetries resulting from flow instabilities (large-scale velocity field) and models small-scale (subgrid-scale) structures as a function of the resolved velocity field (Stoesser, 2014; Rodi, 2017). LES has already proved to be a reliable tool for modelling complex flow around solid obstacles (Grigoriadis et al., 2003; Stoesser et al., 2010; Grigoriadis et al., 2013; Akbari and Montazerin, 2013; Chang and Constantinescu, 2015). The motivation of the present work arises from the importance and the challenge of modelling turbulent flows within random arrays of cylinders. To the authors knowledge, studies that address the conditions for sufficiently accurate LES modelling of flows within random arrays with resolved cylinders have not been presented yet. There are numerical works, on flows within this type of complex boundaries, in which the array of cylinders is represented as distributed drag (as, for example, in Yan et al., 2017b), however the limitations of that approach have not been deeply discussed. In this work, we aimed at modelling infinite domains of randomly placed emergent cylinders, as a proxy to large spatial domains of rigid vegetation. As premise, we assume that an infinite domain can be reproduced as a domain of finite size repeated periodically, as, computationally, infinity means a flow passing for a long time through enough spatial diversity. The main goal of this work is to address the following questions: (i) what is the appropriate grid resolution to guarantee the momentum and the turbulent kinetic energy (TKE) verisimilitude of the flow?, and (ii) how small can the numerical domain be to ensure that the spatial complexity of an infinite array is predicted? Exploring the completeness offered by a numerical database, we also aimed at providing a detailed spatial characterisation of the first and second order moments of the flow. We focus especially the layer controlled by the vertical cylinders, herein named stratum drag-wake layer as the momentum slope is governed or mainly influenced by the fluid interaction with cylinders. Laboratory studies on flows within random arrays of cylinders have highlighted the complexity of the flow-cylinders interaction. Ricardo et al. (2014b) showed that production and dissipation rates of TKE are not locally balanced creating complex spatial patterns of turbulence transport, leading to the generation of background turbulence. Ricardo et al. (2016b) showed that background turbulence impacts the detachment of vortices formed by unsteady separation of the cylinders viscous boundary layer and contributes to the faster loss of coherence of those vortices. Ricardo et al. (2016a) discussed the subsistence of time-averaged flow complexity generated by the upstream cylinders distribution. They concluded that Reynolds stresses are completely determined by the local number-density of cylinders showing that turbulent variables adjust rapidly to the local conditions. This observation points to the possibility of reproducing numerically a large region of randomly placed cylinders with a constant number-density by means of a finite domain periodically repeated. The methodology to tackle the stated goals employs the LES model developed by Grigoriadis et al. (2003, 2004) that has been validated for relevant flow cases (details in Section 3). The numerical data will be compared with laboratory data featuring instantaneous velocity fields were acquired with a two-dimensional, two-component (2D-2C) Particle Image Velocimetry system (PIV). The tested domain is characterised by a number density of cylinder equal to 980 cylinders/m2 randomly distributed but obeying a spatially constant areal density. The effect of the size of the domain and the grid resolution on the modelled flow variables are, herein, discussed independently. An intermediate grid resolution is employed to simulate four domains with different sizes, varying the number of cylinders modelled. The paper is organised as follows. In Section 2 the laboratory experiments are presented. The numerical simulations, with the related computational parameters are presented in Section 3. The impact of the domain size and the grid resolution are discussed in Sections 4 and 5. Section 6 characterises the spatial distribution of first and second order moments in the drag-wake controlled stratum. At the end of the document, the main conclusions are summarised. Throughout this work, a Cartesian reference frame is considered, where x, y, and z correspond to the longitudinal (streamwise), lateral (spanwise), and vertical (bed-normal) directions (see Fig. 1), u, v , and w represent the corresponding velocity components and ωx , ωy and ωz the vorticity components, respectively. The coordinate system refers to the system defined by the flow domain used by the numerical simulations and the distances are normalised ( )by the cylinders diameter. Concerning the notation, the overbar (¯ui ) represents the time-average operator and primes u′i the fluctuation relatively to the time-averaged variable. The flow variables are presented dimensionless, normalised with the bulk velocity (defined as U0 = Q /(bh) where Q is the discharge, b is the channel width and h is the flow depth). The shear stresses are considered as stresses per unit mass. The limits of the drag-wake controlled layer are defined by the inflexions of the time and spaceaveraged profile of the longitudinal turbulence intensity, that, for the flow tested herein, are approximately at z /h = 0.1 and z /h = 0.9 (Ricardo et al., 2016a), corresponding to z /d = 0.5 and z /d = 4.2.

A.M. Ricardo et al. / Journal of Fluids and Structures 80 (2018) 245–261

247

Fig. 1. Plan view of part of the laboratory cylinder array comprising the reach where the measurements were carried out. The measuring gap is between 2.87 and 2.91 m. The solid lines aligned with flow direction identify the location of the vertical planes where two-component velocities were acquired. The top and right hand side axis (grey) correspond to the non-dimensional coordinate system employed in the numerical simulations. The rectangles represent the boundaries of the numerical domains tested. The grey shaded area identifies the part of the numerical domain where the statistics were collected. The flow direction is from left to right.

2. Laboratory experiments The laboratory test was carried out in a recirculating tilting flume of the Laboratory of Hydraulics and Environment of Instituto Superior Técnico, Universidade de Lisboa. This is a 12.5 m long, 41 cm wide and 50 cm deep flume with adjustable slope between −0.5% and +2.5%. The flume has glass side walls, enabling flow visualisation and laser illumination. In the present test, the flume bottom was horizontal and covered with a horizontal layer of gravel of mean diameter D50 = 7.5 mm and sand of D50 = 0.8 mm. Rigid, slightly rough, vertical and circular cylinders were randomly placed along of a 3.5 m long reach with a cylinder number-density of 980 cylinders/m2 . The diameter of the cylinders was 1.1 cm. Downstream the reach populated with cylinders, a coarse gravel weir controlled the flow, which was subcritical both downstream and upstream of the cylinder covered reach. Fig. 1 illustrates part of the cylinder-covered reach in the area where the velocity measurements were carried out. The measuring gap (narrow region without cylinders in the spanwise direction) was enforced with a width approximately equal to the mean inter-cylinder distance to enable the measurements. Measurements consisted in 2D-2C instantaneous velocity maps, in horizontal and vertical planes, with a Particle Image Velocimetry system (PIV). The PIV system is composed of an 8-bit 1600 × 1200 px2 CCD camera and a double-cavity Nd-YAG laser with pulse energy of 30 mJ at wavelength of 532 nm. PIV image pairs were acquired at a frequency of 15 Hz with a time delay of 1500 µs between frames. The instantaneous velocity maps were acquired at eight vertical planes, shown as parallel lines in Fig. 1, and one horizontal plane (82% of the flow depth). Each acquisition recorded 5000 image pairs corresponding to 5′ 33′′ of consecutive data and the field of view was approximately 12 × 10 cm2 , corresponding to interrogation areas of about (1.2 × 1.2) mm2 . Polyurethane particles with a mean diameter of 60 µm in a range from 50 µm to 70 µm and a density of 1.31 g/cm3 were employed as seeding material. Using this seeding, the cut-off frequency of the turbulent signal, calculated with the theory revised by Ferreira and Aleixo (2017), is about 100 Hz. The Nyquist frequency of the PIV measurements is 7.5 Hz, therefore the seeding particles ensure the quality of the data acquired in the time domain. Regarding the space domain, assuming a mean velocity of 0.1 m s−1 and applying Taylor frozen turbulence hypothesis, the cut-off frequency of 100 Hz implies that the smallest resolved turbulence scale is lc = 1.0 mm. This length scale is smaller than the size of the interrogation windows ensuring also the quality of the data in the space domain. For further details on the PIV system and seeding particles employed see Ricardo (2014). The laboratory test featured a flow with a mean longitudinal velocity of U0 = 0.099 m s−1 and flow depth h = 5.1 cm in the measuring gap. The relevant Reynolds number based on the cylinders diameter, defined as Red = U0 d/ν where d is the cylinder diameter and ν is the water kinematic viscosity, was equal to 1064. The flow was steady and gradually√ accelerating with a streamwise depth gradient of dh/dx = −0.014. The flow was sub-critical with Froude number Fr = U0 / gh = 0.14. The free surface exhibited an oscillating behaviour with amplitude approximately 0.5 cm. The most energetic mode of these oscillations occur at f ≈ 2 Hz, which is close to the vortex shedding at the vertical cylinders. During PIV measurements, a 0.2 mm acetate sheet was placed at the free surface, to avoid laser sheet scattering and to promote optical stability. The introduction of this acetate sheet creates a thin boundary layer and acts as a spatial and temporal filter, suppressing low amplitude and low wavelength oscillations. The combination of free surface oscillations and the acetate sheet do not have a measurable impact on the flow dynamics on the drag-wake layer. Only the mean vertical velocity component is expected to be affected by the fluctuations of the free surface.

248

A.M. Ricardo et al. / Journal of Fluids and Structures 80 (2018) 245–261

3. Numerical approach 3.1. Governing equations Under the LES approach, the flow structures are separated into large energy-containing eddies and subgrid scale (SGS) motions. The former are directly resolved by the computational grid while the subgrid scales are parametrised using appropriate turbulence models. Using d and U0 as the characteristic scales for length and velocity, the resulting nondimensional equations of motion for incompressible flow have the following form:

∂ uˆ i =0 ∂ xi

(1)

∂ uˆ i ∂ uˆ i uˆ j ∂p ∂τij 1 ∂ 2 uˆ i + =− − + + fi ∂t ∂ xj ∂ xi ∂ xj Red ∂ xj ∂ xj

(2)

where xi are the Cartesian coordinates, uˆ i are the resolved velocity components in the corresponding directions, p is the ˆ i uˆ j are the SGS stresses. The term fi , is an external forcing term which is used to impose the ˆ pressure field, τij = u i uj − u proper boundary conditions, using the immersed boundary (IB) method. ∂τ In the present work the numerical simulations were performed without SGS model, i.e., the term ∂ xij in Eq. (2) was set j equal to zero. The range of scales resolved by the computational grid proved to cover sufficiently small scales allowing the absence of a dissipative model, corroborating the indication that, in hydraulic flows, the effect of the SGS-model is relatively small in the presence of appropriate grid resolution (Stoesser, 2014). The effect of the SGS model parametrisation was examined carefully during a set of preliminary LES simulations using two different turbulence models, namely the classical Smagorinsky and the Filtered Structure function model (as in Grigoriadis et al., 2012). Employing the smallest domain with 270 × 384 × 97 cells with and without SGS model, we observed differences of about 3%. The averaged square difference was smaller than 6 × 10−4 in the time-averaged velocities, 10−4 in turbulence intensities and 4 × 10−5 in Reynolds shear stresses. Thus, for the Reynolds number considered and the resolutions used in the present study, the effect of the SGS model is not significant. 3.2. Numerical method Spatial discretisation is performed with a second order central finite difference scheme on Cartesian grids. A staggered grid arrangement is followed with pressure nodes located at the cell centres and velocities defined along the corresponding cell face for each component. Time advancement is achieved with a pressure correction algorithm using a fully explicit second order Adams–Bashforth scheme and the fractional step approach. The simulation time step is adjusted dynamically according to the convective (CFL) and viscous time scale (VSL) criteria: CFL < 0.2 and VSL < 0.05. The efficiency of the numerical simulations is mainly because of the direct solution of the derived Poisson’s equation and the adoption of the IB method to describe the cylinders. A fast parallel direct Poisson’s solver is used for the pressure correction which is optimised for shared memory architectures with an operation-oriented parallelisation. The description of solid boundaries with the IB method allows the simulation of geometrically complicated domains with Cartesian solvers (Lai and Peskin, 2000). In IB method, the presence of solid bodies is mimicked by means of suitably defined body forces (term fi in Eq. (2)) applied to the discretised set of the momentum equations. In doing so, one can use non-boundary conforming grids, without having to modify the coefficients of the Poisson equation locally. The body forces are computed at every time step for each velocity component, so that the velocity field on an arbitrary, immersed surface is driven to a specified velocity value. This is achieved with appropriate forcing schemes applied on cells adjacent to the immersed boundary. In case forces on the immersed boundary are needed, a second order local velocity reconstruction scheme is applied along the immersed body, as proposed by Yang and Balaras (2006). This local reconstruction is applied both for the pressure and the velocity fields in the vicinity of the cylinders. The resulting code has an excellent parallel efficiency, reaching performances of 200 MDof (millions degrees of freedom) per second. This means generating approximately 45 consecutive flow-fields (snapshots, or time steps) within a minute when executed on a 16 core personal computer or on a single computational node (equipped with 32 cores and two AMD8350 CPU processors) on a grid with 1024 × 512 × 129 cells, or 67.3 million cells. The physical (RAM) memory required by the algorithm is also low, at the order of 90 MB per million nodes computed. Further details of the numerical implementation can be found in Grigoriadis et al. (2003, 2004). The numerical code has been validated for several flow cases including the unidirectional flow over dunes (Grigoriadis et al., 2012) or the oscillatory flow over ripples (Grigoriadis et al., 2013). 3.3. Boundary conditions and simulation setup Five numerical setups were tested, changing the length (Lx ) or the width Ly of the domain. The vertical size of the domain (Lz ) was kept constant and equal to the flow depth of the laboratory flow, Lz = 4.636d. The relative location of the cylinders, in all numerical simulations, corresponds exactly to the location in the laboratory experiments (Fig. 1).

( )

A.M. Ricardo et al. / Journal of Fluids and Structures 80 (2018) 245–261

249

Table 1 Computational grids and numerical resolution in space and time using the domain DM . Grid

Resolution

Min cell from bed (mm)

d×d (cells)

Time step (ms)

G1 G2 G3 G4

512 × 256 × 65 768 × 384 × 97 1024 × 512 × 129 1536 × 768 × 192

0.418 0.220 0.110 0.077

22 × 25 32 × 38 43 × 51 64 × 76

2.3 1.7 1.2 0.8

Periodic boundary conditions are used along x–y directions and a dynamically computed pressure gradient is adjusted at each time step, so that a constant volumetric flow rate is maintained and the bulk Reynolds number is imposed. The surface of each cylinder was considered as hydrodynamically smooth throughout their height. The numerical domain at the bottom was defined as a horizontal smooth bed with non-slip conditions. We must note that the simulated conditions are not exactly equal to those used in the laboratory test, in terms of roughness. However this work mainly targets on the discussion of flow variables within the drag-wake controlled stratum, which are not significantly affected for the range of parameters examined. This observation is supported by the results presented by Ferreira et al. (2009) and Ricardo et al. (2016a), that considered a sand and gravel bed, respectively. Along the top of the domain, a rigid lid was employed to fix a uniform flow depth equal to the flow depth measured in the laboratory experiments. Thus, a slip boundary condition was employed with du/dz = 0, dv/dz = 0 and w = 0. All numerical simulations were performed in a non-dimensional fashion where distances are scaled with the cylinder diameter d = 11 mm, velocities with U0 = 0.099 m s−1 , leading to a time scale defined as t0 = d/U0 = 0.111 s. Starting from a plug flow as an initial field, the flow developed quickly to a statistically stationary state within 100–150 time units t0 , because of the strong interaction of the bed boundary layer with the shear layers of the cylinders. After achieving the statistically steady state, statistics were collected for a period of 600 time units, i.e. over a time interval of approximately 66.7 s. The temporal evolution of the mean bed shear stress was employed to ensure that this interval was sufficiently large to estimate second order moment statistics. 4. Effect of grid resolution To assess the effect of the grid size on the predicted flow fields, four different grids were employed to model the medium size domain DM , which corresponds to a box with dimensions 24d × 10d × 4.64d (longitudinal × lateral × vertical) containing a total of 26 cylinders (Fig. 1). The numerical tests with different grids were performed in non-dimensional fashion and with the boundary conditions mentioned above. Herein the statistics were collected for a period of 500 time units. Table 1 shows the parameters of four numerical grids tested for this setup. The comparison, shown in Fig. 2, of flow variables computed with different grids reveals that the effect of the grid resolution on the magnitude of the simulated variables is relatively small for the meshes tested herein. It should be noted that the tested meshes are all of relatively high resolution, for the Reynolds number considered, what explains the small difference in the results. The impact of the mesh on the turbulence intensity is larger than in the time-averaged velocity field, with exception of the vertical velocity component. To have a quantitative analysis of the grid effect, the first and second order moments were compared for any combination of two meshes (G1–G2, G1–G3, G1–G4, G2–G3, G2–G4, G3–G4) by computing the square difference, (φGi − φGk )2 , where φGi stands for any flow quantity computed with grid Gi. Herein, this difference was computed, for time-averaged flow variables (first and second order moments), for all points of a common mesh (coarser than G1) defined in the volume where the statistics were collected (grey shaded zone of Fig. 1). The values φGi in the common mesh were calculated with interpolation. This analysis demonstrated that the differences are bounded and have values within order of magnitude 10−3 or smaller. In case of removing the coarsest grid (G1) from the comparisons, the mean square differences become bounded by 5 × 10−4 . Fig. 3 exemplifies the mean square difference (space averaging (φGi − φGk )2 for all points) corresponding to the longitudinal turbulence intensity. Although the differences present the same order of magnitude for all grid combinations, G1 stands out as the mesh resolution with worst results. The analysis of the flow quantities modelled with different grid, as exemplified in Fig. 2 indicates that any of the tested grids would lead to fair estimates of the first and second order moments of the tested flow. The quantitative analysis stressed that the coarsest grid employed, G1, is associated to the largest errors while the other meshes lead to similar differences. This analysis indicates that, for the tested Red , grid resolutions with 30 or more cells are advisable. 5. Effect of domain size To assess the impact of the number of cylinders modelled on the flow variables, five different domains were simulated, DXS , DS , DM , DL and DLy , as shown in Fig. 1. Table 2 summarises the characteristics of the tested setups. The intermediate grid resolution G2 was employed. The numerical tests were performed with the conditions mentioned in Section 3.3, collecting

250

A.M. Ricardo et al. / Journal of Fluids and Structures 80 (2018) 245–261

Fig. 2. Time-averaged profiles of: (a) longitudinal velocity, (b) lateral velocity, (c) vertical velocity, (d) longitudinal turbulence intensity, (e) lateral turbulence intensity and (f) vertical turbulence intensity computed with different grid resolutions. All profiles are longitudinally placed at x/d = 16.5. Lateral profiles (a, b, d and e) were sampled at z /d = 3.8 while vertical profiles (c and f) are located at y/d = 3.4.

Fig. 3. Mean square differences of the longitudinal turbulence intensity computed with any two meshes.

statistics for a period of 600 time units. In the simulation with the smallest domain, DXS , the application of periodic boundary conditions results in a distribution of cylinders in the close vicinity of the measuring gap that is different of that in the other domains and the laboratory. Thus, direct comparisons with results from DXS require particular caution. Figs. 4 and 5 display time-averaged lateral profiles of the first and the second order moments for each tested domain, sampled at z /d = 3.8. Profiles at two longitudinal positions are presented; profiles at x/d = 16.5 (Fig. 4), close to the upstream end of the measuring gap, and further downstream at x/d = 18.5 (Fig. 5). When possible, comparison with the laboratory data is also presented (asterisks markers). Fig. 6 shows the time-averaged vertical profiles of longitudinal and vertical velocities and turbulence intensities, illustrating that vertically the flow variables are approximately constant

A.M. Ricardo et al. / Journal of Fluids and Structures 80 (2018) 245–261

251

Table 2 Parameters for the computational domains shown in Fig. 1 and computational time to compute one second of the flow. Domain

Resolution

Size Lx × Ly × Lz

No. of cylinders

Time step (ms)

Simulation time per second computed

DXS DS DM DL DLy

270 × 384 × 97 480 × 384 × 97 768 × 384 × 97 1080 × 384 × 97 480 × 1024 × 97

8.44 × 10 × 4.636 15 × 10 × 4.636 24 × 10 × 4.636 33.75 × 10 × 4.636 15 × 28 × 4.636

9 16 26 36 44

0.31 0.31 0.21 0.15 0.31

16 min 30 min 50 min 1h 20 min 1h 45 min

throughout the drag-wake controlled layer in all the scenarios. The comparison of time-averaged profiles of numerical and laboratory data is a demanding comparison, however we considered it to be more fruitful for the numerical modelling community than the comparison of spatial averages that would encompass all the flow complexity and uncertainty. Figs. 4 and 5 allow the conclusion that the impact of the effect of the domain length is stronger than the domain width effect. The profiles obtained for DS and DLy which have the same streamwise extent, Lx , and different width, Ly , show, generally, quasi-coincident values. Regarding the mean flow field, Figs. 4 and 5(a) and (b) show that longitudinal and lateral velocity components are well predicted by the numerical model independently of the simulated domain size. The largest differences among the numerical tests are observed in the profiles sampled at x/d = 18.5 (downstream end of the measuring gap), where also the differences to laboratory data are larger. The mean vertical velocity, w ¯ (Figs. 4c and 5c), seems more sensitive to the domain size that u¯ and v¯ . As lateral profiles of w ¯ are not available in the laboratory database a definitive conclusion is not possible. Nevertheless, it should be noted that this flow variable is, in the present study, particularly sensitive to the simplifications considered in the numerical simulations. Furthermore, w ¯ is very small, presenting magnitudes at least one order magnitude smaller than u¯ and v¯ . For the second order moments, although the differences are small, the sorting according to the domain length is more visible than for the velocities, with the exception of DXS for which the results do not reveal a consistent pattern. The profiles for DS and DLy , that almost overlap, show generally smaller values than profiles for DM and DL , which are closer to each other than to the formers. According to Ricardo et al. (2014b), in flows within arrays of cylinders, the cumulative effect of convection and turbulent transport of turbulent kinetic energy leads to the generation of background turbulence, that subsists further downstream and interacts with the locally generated turbulence. Therefore, it is expected to have differences on the second order moments when modelling one or two rows of cylinders or when a large number of cylinder is considered. √

Fig. 7 illustrates the highest values of u′ u′ /U0 computed with domains DXS , DS , DM and DL (the same Ly and increasing Lx ). This figure supports the conclusion that the magnitudes are similar and that the spatial distribution exhibit differences. Moreover, it shows that at a given stage increasing the number of cylinders do not change the turbulence level on the wakes. This is visible from DM to DL . When a large degree of detail on the flow complexity is required, the domains DXS and DS might not capture correctly that complexity. The analysis of the effect of the numerical domain size leads to the conclusion that good predictions of first order moments of a flow within an infinite region populated with cylinders randomly placed are achievable by modelling a small number of cylinders with periodic boundary conditions. Fair estimates of the second order moments are also possible, however losses in quality are expected due to the background turbulence and its interaction with the locally generated turbulence. Although, the minimum number of cylinders required may depend on the number-density and on the Reynolds number, the present work brings evidence that it is possible to forecast the mean flow velocities and Reynolds stresses with sufficient accuracy setting the domain to the largest size ensuring a reasonable simulation time. The stated conclusion is of practical importance for the community of CFD modellers, as they can henceforth minimise the sensitivity analysis of domain size without decreasing the level of results confidence. The flow characterisation that follows is based in data produced with the intermediate domain size, DM . We chose domain DM not only on the basis of results performance but on the basis of computational time, i.e., the option is for the largest domain that guarantees affordable computational time. 6. Characterisation of the time-averaged flow Comparisons of time-averaged profiles of numerically produced data and the laboratory data at selected locations were presented above leading to a positive evaluation of the numerical model ability for modelling flows within infinite arrays of randomly placed cylinders. In the present section, the comparison of numerical and laboratory data will be further commented at the light of the spatial distribution of first and second order moments. Employing the numerical database as valuable extension to the data that was possible to measure in the laboratory, a detailed characterisation of the flow is also presented. The numerical database employed herein corresponds to the simulation with domain DM and grid G4 (1536 × 768 × 192) and for this setup the statistics were collected in a box larger than in the previous simulations, so that the flow characterisation outside the laboratory measuring gap is also available.

252

A.M. Ricardo et al. / Journal of Fluids and Structures 80 (2018) 245–261

Fig. 4. Time-averaged lateral profiles of: (a) longitudinal velocity, (b) lateral velocity, (c) vertical velocity, (d) longitudinal turbulence intensity, (e) lateral turbulence intensity, (f) vertical turbulence intensity, (g) longitudinal-lateral shear stress, (h) longitudinal-vertical shear stress and (i) lateral-vertical shear stress. These profiles were sampled at z /d = 3.8 and x/d = 16.5.

Fig. 8 shows vertical distributions of velocity and turbulence intensities at y/d = 1.4 computed with laboratory and numerical data. The vertical distribution of the flow variables within arrays of cylinders is typically characterised by three layers: the layer close to the bottom, where the flow is 3D due to interaction with the bed; a thin layer affected by the free surface oscillations and the intermediate layer, drag-wake controlled stratum, where the flow is controlled by the vertical elements and the flow properties are approximately constant in the vertical direction (Ricardo et al., 2016a). Fig. 8 shows that in the drag-wake controlled stratum the flow is well predicted by the numerical model. Figs. 9 and 10 present the spatial distributions of flow variables computed with numerical and laboratory databases in the horizontal plane placed at z /d = 3.8. Fig. 11 shows the spatial distribution, in the drag-wake controlled stratum, of flow variables that are, generally, not available in laboratory databases. The time-averaged velocity field is characterised by low longitudinal velocity in the wake of the cylinders alternated by high longitudinal velocity, u¯ , zones between cylinders (Fig. 9a and b), evidencing heterogeneity at large scales. The lateral velocity, v¯ , is at least one order of magnitude smaller than u¯ and its largest values occur in the close vicinity of the cylinders,

A.M. Ricardo et al. / Journal of Fluids and Structures 80 (2018) 245–261

253

Fig. 5. Time-averaged lateral profiles of: (a) longitudinal velocity, (b) lateral velocity, (c) vertical velocity, (d) longitudinal turbulence intensity, (e) lateral turbulence intensity, (f) vertical turbulence intensity, (g) longitudinal-lateral shear stress, (h) longitudinal-vertical shear stress and (i) lateral-vertical shear stress. These profiles were sampled at z /d = 3.8 and x/d = 18.5.

predominantly upstream of each cylinder (Fig. 9c and d), in patterns of negative and positive values. The lateral velocity presents larger magnitudes between x/d = 10 and x/d = 15 than within the measuring gap, as due to the proximity between cylinders the flow is forced to deviated finding preferential paths which are not parallel to the streamwise direction. The vertical velocity, illustrated in Fig. 11(a), is vanishing small, the largest magnitudes being found in the wake regions. Fig. 9(e) and (f) show the bed-normal vorticity distribution, revealing a quasi-antisymmetric pattern on the wake of the cylinders identifying von Kármán vortex streets of the paired vortices caused by the unsteady separation of the flow on the cylinder. As previously reported by Ricardo and Ferreira (2016) and Ricardo et al. (2014a), the vorticity distribution in the wakes of cylinders within arrays is similar to that on wakes of isolated cylinders. However, within arrays of cylinders, the vortex streets appear suppressed when the vortex streets have less space to develop due to the local distribution of the neighbouring cylinders. Fig. 9 show that the flow patterns in the wake of cylinders in random arrays do not exhibit the quasi-perfect symmetry characteristic of the wakes behind isolated cylinders. The presence of the neighbouring cylinders conditions the flow circulation behind a given cylinder.

254

A.M. Ricardo et al. / Journal of Fluids and Structures 80 (2018) 245–261

Fig. 6. Time-averaged vertical profiles of: (a) longitudinal velocity, (b) vertical velocity, (c) longitudinal turbulence intensity and (d) vertical turbulence intensity. These profiles were sampled at x/d = 17.5 and y/d = 3.4.

Fig. 7. Contours of the highest magnitudes of longitudinal turbulence intensity in the numerical domains DXS , DS , DM and DL , at z /d = 3.8.

The comparison between laboratory and numerical results corroborates the above conclusion about the suitability of the employed LES approach to model such a complex flow. The magnitude and the spatial distribution of the first order moments are well predicted, the largest differences being found at the domain lateral boundaries where the local distribution of cylinders is not equal. The direct comparison between numerical and laboratory data, against each other, indicated that only for the vertical velocity component the agreement is poor. The numerical values of w/ ¯ U0 are nearly zero everywhere while in the laboratory database the range of variation is from w/ ¯ U0 = −0.2 to w/ ¯ U0 = 0.2. The mean square difference

A.M. Ricardo et al. / Journal of Fluids and Structures 80 (2018) 245–261

√

(a) Lab: u/U0 .

(b) Lab:

(d) Num: u/U0 .

(e) Num:

u′ u′ /U0 .

√

u′ u′ /U0 .

255

(c) Lab:

√

(f) Num:

w′ w ′ /U0 .

√

w′ w ′ /U0 .

Fig. 8. Two-dimensional maps of longitudinal velocity (left), longitudinal turbulence intensity (centre) and vertical turbulence intensity (right) at y/d = 1.4 of the laboratory (top) and numerical (bottom) databases. The flow direction is from left to right.

between numerical and laboratory data, considering the data in the drag-wake stratum from vertical and horizontal planes where laboratory data is available, is 2.7 × 10−2 , 7 × 10−3 and 3 × 10−3 for u¯ /U0 , v¯ /U0 and w/ ¯ U0 , respectively. The pattern of the spatial distribution and the magnitude of the second order moments presented in Fig. 10 are also fairly −3 well predicted by the numerical √ difference between numerical and laboratory data is 4 × 10 , √ model.√The mean square

6 × 10−3 and 3 × 10−3 for u′ u′ /U0 , v ′ v ′ /U0 and w ′ w ′ /U0 , respectively. The lateral variation is accurately forecast, i.e., the peaks and troughs are found at the same y−location both in the experimental and numerical domains. Differences are found in the location of the variables maximum relatively to each cylinder: in the laboratory database, the highest values were found closer to the cylinder and seem to decay faster than in the numerical database. The regions of highest magnitudes for the laboratory turbulence intensities and shear stress are slightly shorter, in the streamwise direction, than the corresponding numerical distributions. This difference might be explained by the absence of a turbulence SGS model or by the summing up effect of the small differences in the laboratory and modelled flow. Nevertheless, considering the simplifications included in the simulation, differences between numerical and laboratory results are small, as Figs. 9 and 10 show. The spatial distribution of the longitudinal turbulence intensity, presented in Fig. 10 (a and b), shows the largest magnitudes on the vortex street downstream of each cylinder and lower values in the corridors between cylinders. While the distribution of the lateral (Fig. 10c and d) and vertical (Fig. 11b) turbulence intensities present larger magnitudes on the wake of each cylinder. Therefore, these patterns are similar to the symmetric distribution of turbulence intensities of a flow around an isolated cylinder. Within the array the symmetry is disrupted by the presence of the neighbouring cylinders, what can be observed in Figs. 10(a–d) and 11(b) particularly in the region between x/d = 10 and x/d = 15. Regions with almost zero turbulence intensities that can be observed in the regular distributions of cylinders (Stoesser et al., 2010) practically do not exist in random distributions. Figs. 10(a–d) and 11(b), corroborating with Ricardo et al. (2014b), indicate that the randomness of the cylinders placement promotes an efficient generation of background turbulence, relatively to the case of squared and staggered distributions. Fig. 10(e) and (f) present the non-dimensional longitudinal-lateral component of the Reynolds shear stresses, u′ v ′ /U02 . It shows an almost anti-symmetric pattern relatively to the axis of each cylinder of high positive and negative values on the wake of each cylinder, and almost zero values between cylinders. The longitudinal-vertical and lateral-vertical components of the turbulent shear stress, Fig. 11(c) and (d), are much smaller than u′ v ′ /U02 , exhibiting the largest magnitudes within the cylinders wake. The spatial distribution of the longitudinal and lateral components of the time-averaged vorticity field is, generally, difficult to characterise with laboratory data. Fig. 11(e) and (f) represent these flow variables within the drag-wake controlled

256

A.M. Ricardo et al. / Journal of Fluids and Structures 80 (2018) 245–261

(a) Num: u/U0 .

(b) Lab: u/U0 .

(c) Num: v/U0 .

(d) Lab: v/U0 .

(e) Num: ωz /dU0 .

(f) Lab: ωz /dU0 .

Fig. 9. Two-dimensional maps of time-averaged longitudinal velocity (top), lateral velocity (centre) and vertical vorticity (bottom) at z /d = 3.8 of the numerical (left) and laboratory (right) databases. The dotted line identifies the limits of the laboratory measuring gap, easing numerical/laboratory comparisons. The flow is from left to right.

stratum, where their magnitude is nearly zero. Within this layer, the most important coherent structures are the vortices of vertical axis generated by the shedding on the cylinders. To further highlight the capabilities of the numerical model, the database produced is employed to the visualisation of the time-averaged flow structures around the circular cylinders, with particular interest on the vortex system associated to

A.M. Ricardo et al. / Journal of Fluids and Structures 80 (2018) 245–261

257

(a) Num:

√

u′ u′ /U0 .

(b) Lab:

√

(c) Num:

√

v ′ v ′ /U0 .

(d) Lab:

√

(e) Num: u′ v ′ /U02 .

u′ u′ /U0 .

v ′ v ′ /U0 .

(f) Lab: u′ v ′ /U02 .

Fig. 10. Two-dimensional maps of longitudinal turbulence intensity (top), lateral turbulence intensity (centre) and longitudinal-lateral turbulent shear stress (bottom) at z /d = 3.8 of the numerical (left) and laboratory (right) databases. The dotted line identifies the limits of the laboratory measuring gap, easing numerical/laboratory comparisons. The flow is from left to right.

the junction flow, around the cylinders and close to the bed, which require high-resolution 3D databases. Figs. 12 and 13 illustrate the flow structures close to the bed employing the Q -criterion method. It is possible to identify the horse shoe vortex due to separation and roll-up of boundary layer vorticity in front of the cylinder and the counter-rotating vortices (base vortices) featured by the base flow around each cylinder (Porteous et al., 2014).

258

A.M. Ricardo et al. / Journal of Fluids and Structures 80 (2018) 245–261

(a) Num: w/U0 .

(b) Num:

√

w ′ w′ /U0 .

(c) Num: u′ w ′ /U02 .

(d) Num: v ′ w ′ /U02 .

(e) Num: ωx /dU0 .

(f) Num: ωy /dU0 .

Fig. 11. Numerically computed two-dimensional maps of time-averaged: (a) vertical velocity, (b) vertical turbulence intensity, (c) longitudinal-vertical turbulent shear stress, (d) lateral-vertical turbulent shear stress, (e) longitudinal vorticity and (f) lateral vorticity at z /d = 3.8. The flow is from left to right.

7. Conclusions Numerical simulations were performed and compared with PIV experiments of a fully developed flow around randomly placed emergent circular cylinders. The numerically computed time-averaged velocity field and Reynolds stress tensor, in the drag-wake controlled stratum, were found in fair agreement with the respective laboratory variables.

A.M. Ricardo et al. / Journal of Fluids and Structures 80 (2018) 245–261

259

Fig. 12. Iso-surface of the Q -criterion (Q = 4) for identification of three-dimensional time-averaged flow structures. The vertical axis is scaled by a factor of 2.

(a) z /d = 0.05.

(b) z /d = 0.10.

(c) z /d = 0.15.

(d) z /d = 0.20.

(e) z /d = 0.50.

(f) z /d = 2.50.

Fig. 13. Bed parallel planes of positive Q (Q -criterion method) for identification of three-dimensional time-averaged flow structures.

The numerical database configured a valuable extension of the laboratory database allowing for data in a spatial domain wider than the experimental measuring gap, including the near vicinity of cylinder surfaces, 3D data and three components of the velocity field and derived turbulence quantities. The spatial characterisation of the flow at the drag-wake controlled stratum revealed that: – the three turbulence intensity components exhibit values within the same order of magnitude, but the vertical component is smaller than the longitudinal and vertical components; – the longitudinal-lateral is the most important shear component of the Reynolds stress tensor, with values of the same order of magnitude of the normal stresses; the longitudinal-vertical and lateral-vertical shear components are at least one order of magnitude smaller; – the vertical vorticity is the most relevant component of the time-averaged vorticity field as the longitudinal and the lateral components are vanishing small; – the proximity to the neighbouring cylinders impacts patterns of the spatial distribution of the flow variables causing the lack of symmetry (or anti-symmetry) characteristic of wakes of isolated cylinders. Two questions related to the grid resolution and the domain size directed this work. To answer the question ‘‘what is the appropriate grid resolution to guarantee the momentum and the TKE verisimilitude of the flow?’’ four meshes were tested resolving the cylinder diameter (in the lateral direction) with 25, 38, 51 and 76 grid points. The computed flow variables presented similar values for all the tested meshes with the coarsest grid, G1, being the grid with the largest differences.

260

A.M. Ricardo et al. / Journal of Fluids and Structures 80 (2018) 245–261

All the mean square differences were bounded by 3 × 10−3 and if grid G1 is excluded those are bounded by 5 × 10−4 . For the considered Reynolds number and boundary conditions, a grid with 30 points per cylinder diameter provide appropriate resolution to compute the first and the second order moments of a turbulent flow within random arrays of cylinders. The effect of the size of the numerical domain (or the number of cylinders modelled) was discussed by comparing results of five different numerical domains containing 9, 16, 26, 36 and 44 cylinders. The first order moments were well predicted independently of the simulated domain size. Regarding the second order moments, although the magnitudes are similar in all the domains, the smallest domains lead to losses in quality due to the background turbulence effects. Thus, to the question ‘‘how small can the numerical domain be to ensure that the spatial complexity of an infinite array is predicted?’’, the present work indicates that domain length should encompass the equivalent to five rows of cylinders in the streamwise direction. However our recommendation is to choose the largest domain that guarantees affordable computational time. The general conclusion points out to the possibility of predicting the first and second order moments of flow in very large (infinite) domains of randomly placed cylinders by modelling a relatively small domain with periodic boundary conditions. This is an important conclusion for the CFD community, as it allows for advances in the understanding of these complex flow with minimal computational resources. The present work configured a step forward on the modelling of turbulent flows within arrays of cylinders targeting a random distribution, nevertheless, there are still many challenges to address. As future steps we consider testing distributions with larger number-density of cylinders and different Reynolds numbers. Also improvements on the modelling of the free surface will be included in future works. Overall, the findings presented in this work confirm the capability of efficient numerical modelling techniques to study the dynamics of large-scale flows through vegetation canopies within affordable computational times. Acknowledgements This work was partially funded by FEDER, program COMPETE, and by national funds through Portuguese Foundation for Science and Technology (FCT) project PTDC/ECM-HID/6387/2014 - POCI-01-0145-FEDER-016825 and the grant SFRH/BPD/93903/2013. References Aberle, J., Järvelä, J., 2013. Flow resistance of emergent rigid and flexible floodplain vegetation. J. Hydraul. Res. 51, 33–45. http://dx.doi.org/10.1080/ 00221686.2012.754795. Akbari, G., Montazerin, N., 2013. On the role of anisotropic turbomachinery flow structures in inter-scale turbulence energy flux as deduced from spiv measurements. J. Turbul. 14, 44–70. http://dx.doi.org/10.1080/14685248.2013.861073. Chang, K., Constantinescu, G., 2015. Numerical investigation of flow and turbulence structure through and around a circular array of rigid cylinders. J. Fluid Mech. 776, 161–199. http://dx.doi.org/10.1017/jfm.2015.321. Etminan, V., Lowe, R.J., Ghisalberti, M., 2017. A new model for predicting the drag exerted by vegetation canopies. Water Resour. Res. 53, 3179–3196. http://dx.doi.org/10.1002/2016WR020090. Ferreira, R.M., Aleixo, R., 2017. Experimental Hydraulics: Methods, Instrumentation, Data Processing and Management: Vol. II: Instrumentation and Measurement Techniques. In: Velocity, CRC Press, pp. 35–209 (Chapter 3). Ferreira, R.M.L., Ricardo, A.M., Franca, M.J., 2009. Discussion of ‘‘Laboratory investigation of mean drag in a random array of rigid, emergent cylinders’’ by Heydi M. Nepf and Yukie Tanino, J. Hydraul. Eng., vol. 134, no 1, 2008. J. Hydraul. Eng. 135, 690–693. http://dx.doi.org/10.1061/(ASCE)HY.19437900.0000021. Ghisalberti, M., Nepf, H., 2005. Mass transport in vegetated shear flows. Environ. Fluid Mech. 5, 527–551. http://dx.doi.org/10.1007/s10652-005-0419-1. Grigoriadis, D.G.E., Balaras, E., Dimas, A.A., 2013. Coherent structures in oscillating turbulent boundary layers over a fixed rippled bed. Flow Turbul. Combust. 91, 565–585. http://dx.doi.org/10.1007/s10494-013-9489-1. Grigoriadis, D.G.E., Bartzis, J.G., Goulas, A., 2003. Les of the flow past a rectangular cylinder using the immersed boundary concept. Internat. J. Numer. Methods Fluids 41, 615–632. http://dx.doi.org/10.1002/fld.458. Grigoriadis, D., Bartzis, J., Goulas, A., 2004. Efficient treatment of complex geometries for large eddy simulations of turbulent flows. Comput. & Fluids 33, 201–222. http://dx.doi.org/10.1016/S0045-7930(03)00038-0. Grigoriadis, D.G., Dimas, A.A., Balaras, E., 2012. Large-eddy simulation of wave turbulent boundary layer over rippled bed. Coast. Eng. 60, 174–189. http://dx.doi.org/10.1016/j.coastaleng.2011.10.003. Kim, S.J., Stoesser, T., 2011. Closure modeling and direct simulation of vegetation drag in flow through emergent vegetation. Water Resour. Res. 47, W10511. http://dx.doi.org/10.1029/2011WR010561. Lai, M.-C., Peskin, C.S., 2000. An immersed boundary method with formal second-order accuracy and reduced numerical viscosity. J. Comput. Phys. 160, 705–719. http://dx.doi.org/10.1006/jcph.2000.6483. Lam, K., Lo, S., 1992. A visualization study of cross-flow around four cylinders in a square configuration. J. Fluids Struct. 6, 109–131. http://dx.doi.org/10. 1016/0889-9746(92)90058-B. Marjoribanks, T.I., Hardy, R.J., Lane, S.N., 2014. The hydraulic description of vegetated river channels: the weaknesses of existing formulations and emerging alternatives. Wiley Interdiscip. Rev.: Water 1, 549–560. http://dx.doi.org/10.1002/wat2.1044. Nepf, H., 2012. Flow and transport in regions with aquatic vegetation. Annu. Rev. Fluid Mech. 44, 123–142. http://dx.doi.org/10.1146/annurev-fluid-120710101048. Nicolle, A., Eames, I., 2011. Numerical study of flow through and around a circular array of cylinders. J. Fluid Mech. 679, 1–31. http://dx.doi.org/10.1017/ jfm.2011.77. Porteous, R., Moreau, D.J., Doolan, C.J., 2014. A review of flow-induced noise from finite wall-mounted cylinders. J. Fluids Struct. 51, 240–254. http: //dx.doi.org/10.1016/j.jfluidstructs.2014.08.012. Ricardo, A.M., 2014. Hydrodynamics of Turbulent Flows Within Arrays of Circular Cylinders. (Ph.D. thesis), Instituto Superior Ténico & Ecole Polytechnique Fédérale de Lausanne.

A.M. Ricardo et al. / Journal of Fluids and Structures 80 (2018) 245–261

261

Ricardo, A.M., Di Carlo, S., Franca, M.J., Schleiss, A.J., Sanches, P.M., Ferreira, R.M., 2014a. Vortex interaction in patches of randomly placed emergent cylinders. In: S., , et al. (Eds.), Proc. of River Flow 2014. pp. 63–69. Ricardo, A.M., Ferreira, R.M.L., 2016. Vorticity fluxes on the wake of cylinders within random arrays. In Proceedings of 4th IAHR Europe Congress. Ricardo, A.M., Franca, M.J., Ferreira, R.M., 2016a. Turbulent flows within random arrays of rigid and emergent cylinders with varying distribution. J. Hydraul. Eng. 142, 04016022. http://dx.doi.org/10.1061/(ASCE)HY.1943-7900.0001151. Ricardo, A.M., Koll, K., Franca, M.J., Schleiss, A.J., Ferreira, R.M., 2014b. The terms of turbulent kinetic energy budget within random arrays of emergent cylinders. Water Resour. Res. 50, 4131–4148. http://dx.doi.org/10.1002/2013WR014596. Ricardo, A.M., Sanches, P., Ferreira, R., 2016b. Vortex shedding and vorticity fluxes in the wake of cylinders within a random array. J. Turbul. 17, 999–1014. http://dx.doi.org/10.1080/14685248.2016.1212166. Rodi, W., 2017. Turbulence modeling and simulation in hydraulics: A historical review. J. Hydraul. Eng. 143, 03117001. http://dx.doi.org/10.1061/(ASCE) HY.1943-7900.0001288. Sewatkar, C.M., Patel, R., Sharma, A., Agrawal, A., 2012. Flow around six in-line square cylinders. J. Fluid Mech. 710, 195–233. http://dx.doi.org/10.1017/jfm. 2012.359. Stoesser, T., 2014. Large-eddy simulation in hydraulics: Quo vadis. J. Hydraul. Res. 52, 441–452. http://dx.doi.org/10.1080/00221686.2014.944227. Stoesser, T., Kim, S., Diplas, P., 2010. Turbulent flow through idealized emergent vegetation. J. Hydraul. Eng. 136, 1003–1017. http://dx.doi.org/10.1061/ ASCEHY.1943-7900.0000153. Sumner, D., 2010. Two circular cylinders in cross-flow: A review. J. Fluids Struct. 26, 849–899. http://dx.doi.org/10.1016/j.jfluidstructs.2010.07.001. Tanino, Y., Nepf, H.M., 2008. Lateral dispersion in random cylinder arrays at high reynolds number. J. Fluid Mech. 600, 339–371. http://dx.doi.org/10.1017/ S0022112008000505. Tanino, Y., Nepf, H.M., 2009. Laboratory investigation of lateral dispersion within dense arrays of randomly distributed cylinders at transitional reynolds number. Phys. Fluids 21, 046603. http://dx.doi.org/10.1063/1.3119862. Vargas-Luna, A., Crosato, A., Uijttewaal, W.S., 2015. Effects of vegetation on flow and sediment transport: comparative analyses and validation of predicting models. Earth Surf. Process. Landforms 40, 157–176. http://dx.doi.org/10.1002/esp.3633. Wang, X., Gong, K., Liu, H., Zhang, J.-X., Tan, S., 2013. Flow around four cylinders arranged in a square configuration. J. Fluids Struct. 43, 179–199. http://dx.doi.org/10.1016/j.jfluidstructs.2013.08.011. Yan, C., Huang, W.-X., Miao, S.-G., Cui, G.-X., Zhang, Z.-S., 2017a. Large-eddy simulation of flow over a vegetation-like canopy modelled as arrays of bluff-body elements. Bound. Layer Meteorol. http://dx.doi.org/10.1007/s10546-017-0274-x. Yan, C., Nepf, H.M., Huang, W.-X., Cui, G.-X., 2017b. Large eddy simulation of flow and scalar transport in a vegetated channel. Environ. Fluid Mech. 17, 497–519. http://dx.doi.org/10.1007/s10652-016-9503-y. Yang, J., Balaras, E., 2006. An embedded-boundary formulation for large-eddy simulation of turbulent flows interacting with moving boundaries. J. Comput. Phys. 215, 12–40. http://dx.doi.org/10.1016/j.jcp.2005.10.035. Zhang, H.J., Zhou, Y., 2001. Effect of unequal cylinder spacing on vortex streets behind three side-by-side cylinders. Phys. Fluids 13, 3675–3686. http: //dx.doi.org/10.1063/1.1412245.