Non-hydrostatic modelling of shallow water flow around a circular array of emergent cylinders

Vegetation in fluvial systems plays a significant role in the sustainable development of aquatic systems because of its biological and hydrodynamic effects. It improves the water equality by the production of oxygen and uptake of nutrients. The biodiversity in rivers and oceans is promoted by providing food and creating different habitats. Also, by altering hydrodynamic conditions, such as flow velocity, turbulence structure and mass mixing, the exchanges of sediment, metals and other contaminants between terrestrial and aquatic systems can be well mediated. Aquatic plants can change the velocity field across different scales, ranging from individual branches and blades on a single plant to a patch of plants, called canopy or meadow (Nepf 2012; Aberle & Järvelä 2015). Generally, vegetation often exists in large-scale in coastal regions and rivers, and it forms a complex hydrodynamic environment (Yu et al. 2014). The drag caused by large-scale vegetation in rivers can stabilize sediment and change the patterns of erosion and deposition. In nearshore regions, vegetation dissipates wave energy, which protects shoreline and riverbed from flood and wave attack. Therefore, it is of great significance to investigate the hydrodynamic influence of vegetation on shallow water flows (Wang et al. 2018).

A variety of laboratory experiments on flow through vegetation have been conducted over the last few decades. In these studies, the vegetation elements were usually simplified as rigid bodies, and the vegetation zone was modelled by an array of solid cylinders (Poggi et al. 2004; Nezu & Sanjou 2008; Zhan et al. 2017). Ball et al. (1996) measured the velocity vector fields associated with flow through pile groups in shallow water using particle tracking velocimetry (PTV). A steady low velocity near wake followed by an unstable far wake was observed downstream of the structure. Chen & Jirka (1995) investigated the wake structure behind a porous plate further. Their observations show that the von Kármán vortex street originates at a distance downstream of the porous obstruction. The region from the obstruction to the point where the vortex street forms is called the steady wake. Nicolle & Eames (2011) found that the vortex street did not appear at all for the solid volume fraction (SVF) less than 0.05. Additionally, as the SVF decreases, the vortex street moves further downstream, resulting in a longer steady wake. In the experiments of Zong & Nepf (2012), the vortex development behind a patch of cylinders was studied, and a monotonous relationship between the dimensionless length of wake and the SVF was generalized.

Both small-scale and large-scale turbulence are generated when the water flows through a patch of cylinders. The small-scale turbulence is generated by the stem wake behind single cylinders, and the large-scale turbulence is produced in the wake behind the patch. Takemura & Tanaka (2007) investigated the flow structures and drag characteristics around a colony-type model, which comprises seven equally spaced emergent cylinders. They concluded that the change in the large-scale and small-scale vortex structures was related to the interaction between the water flow and the solid obstruction. The production of patch-scale turbulence is mainly caused by the separated shear layers (SSLs) formed at the shoulders of the patch (Chang & Constantinescu 2015). The patch-scale turbulence is often considered as horizontal 2D or quasi-2D because of the low ratio of water depth to the patch scale. For shallow water flow, the turbulent flow is considered to be constrained by the rough bottom. Chu & Babarutsi (1988) concluded that the bed-friction generates not only the small-scale turbulence but, at the same time, exerts a stabilizing force on the large-scale turbulent coherent structures in shallow turbulent mixing layers. The vortex street behind a solid obstruction may be suppressed by the strong bed friction, which is characterized by a critical stability parameter, S_c = 0.20 (Ingram & Chu 1987). Therefore, the development of patch-scale turbulence of the wake flow can also be impeded by the bed friction.

There is a periodic unsteadiness in the flow when the vortex shedding off bluff bodies occurs. It is one of the toughest tasks for turbulence modelling to accurately predict the unsteady turbulent flows (Durbin 2002; Wegner et al. 2004; Xu & Ma 2009). Durbin (2002) argued that Reynolds-averaged Navier–Stokes (RANS) cannot be applied to unsteady flow, unless there is a spectral gap between the unsteadiness and the turbulence. In the experiments of Zong & Nepf (2012), the vortex shedding frequencies were calculated to be around 2 Hz (cylinder-scale) and 0.1 Hz (patch-scale), and the gap between these shedding frequencies and the background turbulence was also observed. Thus, the RANS methods were still commonly used in practical simulations. Brito et al. (2016) conducted the RANS simulations of compound open-channel flows with vegetated floodplains. Neary (2003) used a RANS-k-ω model to simulate flow through submerged vegetation. Yu et al. (2014) simulated flow through a circular patch of emergent cylinders using a horizontal 2D multi-body model and porous model, respectively. Their achievements agreed well with the experimental data regarding the mean velocities, but failed to predict the magnitude of the turbulent fluctuation. Moreover, the influence of free surface and bed shear were not included in the 2D model. When the patch diameter D is much larger than the water depth H, the bed friction is large enough to suppress the vortex street and cannot be ignored (Zong & Nepf 2012). Under this condition, a robust 3D model should be considered.

In this paper, a 3D non-hydrostatic FVM model (HydroFlow^®) based on unstructured grids is developed and used to simulate shallow water flow through a circular array of emergent cylinders. The flow patterns, turbulence intensity and wake structures are analysed in detail referring to different solid volume fractions and water depths.

A two-step predictor–corrector scheme is used in the numerical method. The flow driven by the hydrostatic pressure is first calculated in the predictor step, and then the flow driven by the non-hydrostatic pressure is updated in the corrector step. In the numerical model, the grid system consists of unstructured meshes in the horizontal plane and multi-layers in the vertical direction. The governing equations are discretized based on the finite volume method (FVM). A second order total variation diminishing scheme (OSHER scheme) is adopted to discretize the convective terms, and the central difference method is used to discretize the diffusion terms. In the corrector step, the Poisson-type equation for the non-hydrostatic pressure is numerically solved by a pre-conditioned BI-CGSTAB approach. The in-house codes (HydroFlow^®) are paralleled by the OpenMP library (Zhang et al. 2014).

The same array geometry from experiments (Zong & Nepf 2012) was chosen in the present work. Their measurements were conducted in a channel with a length of 13 m and a width of 1.2 m, which is shown in Figure 2. A circular patch of cylinders was put in the middle of the channel in a staggered arrangement. The leading edge of the patch was 3.0 m away from the inlet. The diameter d of a single cylinder was 0.6 cm, and the diameter D of the patch was 22 cm. The density of the cylinders within the patch was defined as the solid volume fraction (SVF), and was determined using φ=nπd²/4, in which n (cm⁻²) is the number of cylinders per unit horizontal area.

A steady flow discharge was specified at the inlet of the numerical flume, and a steady water elevation was specified at the downstream outlet. The integration time step was 0.0005 s in all of the study cases. The simulations proceeded until the temporal water elevation probed at a fixed upstream observing point reached nearly steady, i.e., the calculated water elevation fluctuation decreased to ±0.001 m. The computation was continually carried out during the next 60 s; meanwhile, the flow variables were recorded in a frequency of 20 Hz. The saved database was statistically analysed to identify the wake flow patterns.

Ten cases were numerically investigated in this study, with the SVF φ being 0.03 and 0.10, and the water depth H varying from 0.03 m to 0.266 m (see Table 1). U = 0.098 m/s was the free stream velocity at the inlet. The Reynolds number Re_d = Ud/ν was thus calculated to be 588.

The simulated mean velocity profiles (at y = 0) and (at y = D/2) in the longitudinal direction were compared with experiments in Figure 3. To be consistent with the experiments of Zong & Nepf (2012), the velocity was measured at mid-depth. In Figure 3, the outer range of the cylinder patch is outlined from x/D = 0 to x/D = 1.0. In general, good agreements between the predicted and measured velocities were achieved, which validated the accuracy of the present model.

Flow around and through an array of cylinders exhibited a great difference to the flow around a solid body. At the centreline, the streamwise velocity decreased gradually from about 2D upstream of the patch. Within the patch, the streamwise velocity u decreased downstream because of the drag forces induced by cylinders. The mean velocity profiles (Figures 4(a) and 5(a)) downstream a lower SVF patch (φ = 0.03) were nearly the same at different water depths. However, the spatial distributions of the mean streamwise velocity behind a dense patch (φ = 0.10) were distinct with the water depth varying. At large water depths (i.e., P1H5, P1H6, P1H7), there were only small discrepancies between the mean flow in the wakes. As shown in Figure 4(b), the steady near wake extended further downstream and there was a longer distance for the wake velocity deficit recovery at small water depths. The recovery of the centreline velocity is known so that the cross-wake mixing is mainly driven by the patch-scale vortex street. The analysis revealed that the patch-scale turbulence was weakened with H decreasing. The water depth was critical to the evolution of the wake flow; however, the length of the steady wake was found to be nonmonotonic with the water depth H. It was deduced that the change of wake structures was affected by the strong bed friction force. The flow reversal in the wake, which is indicative of a recirculation zone ⁠, was observed at large water depths and high SVF (Figure 4(b)). The recirculation zone moved downstream as H decreased (i.e., P1H3, P1H4, P1H5), but it disappeared at very low water depths (i.e., P1H1, P1H2). Instead of changing the water depth, Zong & Nepf (2012) observed the same trend for the length of the recirculation zone by changing the SVF.

The influence of water depth on the turbulent structures was further illustrated by the profiles of turbulent kinetic energy (TKE) (Figure 6). In the simulations, the TKE comprised two parts, i.e., TKE = k_resolved + k_modelled. The first constituent was obtained using the directly simulated velocity, which was calculated as k_resolved = u_rms² + v_rms². The other constituent was obtained from the turbulence Equations (6) and (7), which were modelled by the turbulence model, i.e., k_modelled = k. For a RANS-type model, the resolved TKE overwhelms the modelled TKE in the wake, i.e., x/D > 1 region, but the modelled TKE overwhelms the resolved TKE within the cylinder patch (x/D < 1) (Yu et al. 2014). For convenience, the turbulent kinetic energy was normalized by the mean velocity and denoted as the turbulence intensity, ⁠.

In the wake of a patch, there were two distinct zones of elevated turbulence. The upstream peak of I_TKE (denoted as I_{TKE, max}) within the patch (0 < x/D < 1) indicated the small-scale turbulence generated by individual cylinders (Figure 6(a)). I_{TKE, max} in Figure 6(b) appeared at a distance away from the patch (x/D > 1), which corresponded to the patch-scale turbulence. As shown in Figure 6, the variation of the small-scale turbulence within the patch was relatively weak. That indicated the impact of bed friction on the cylinder-scale turbulence within the patch was negligible. To facilitate analysis, the cylinder-scale turbulence intensity at the end of the patch (x/D = 1) was characterized by Istem. Commonly, the turbulence intensity I_TKE in the sparse patch (φ = 0.03) was higher than that in the dense patch (φ = 0.1), which was because the wake-to-cylinder interactions were suppressed at high SVF. For high SVF, the stem-scale turbulence decayed rapidly through the patch. Behind the patch, the turbulence intensity decreased in the near wake, and then increased approaching the downstream edge of the near wake. When the SSLs met, the turbulence intensity reached maximum and this peak contributed to the formation of large-scale Kármán vortex street (LKV). The profiles of I_TKE for large water depths (i.e., P1H5, P1H6 and P1H7) differed slightly. However, when H decreased under 0.133 m (i.e., P1H1, P1H2, P1H3 and P1H4), the position of the second I_{TKE, max} moved further downstream, which indicated that the onset of patch-scale turbulence was delayed. The strength of the second I_{TKE, max} decreased with H decreasing for both of the SVFs. The simulations highlighted the turbulence intensity was weakened as H decreased under one critical value, which resulted in a lower rate of the velocity recovery. Because the movement of SSLs contributed to the TKE behind the patch (Chang & Constantinescu 2015), the spatial evolution of SSLs was suppressed by the bed friction, which was gradually distinguished when the water depth decreased.

It is known that the time-averaged vorticity may smear out important information about the flow structures. Thus, the instantaneous vertical vorticity, ω_z was used to investigate the temporal and spatial evolution of the flow through the cylinder patch. When the upstream flow approached the cylinder patch, some water was forced around the patch and the residual water flowed through the patch. Because of the drag force on individual cylinders, there was a velocity deficit between the flow through the patch and the flow around the patch. Hence, two separated shear layers (SSLs) formed at each shoulder of the patch. These shear layers were separated by the flow through the patch (similar to a bleed flow) and extended a distance downstream of the patch until the LKV formed induced by the instability of the shear layer flow. For a high SVF and a large water depth (Figure 7(a)–7(c)), the LKV behind the steady wake was clearly simulated. The LKV was also observed for a small water depth (Figure 7(d) and 7(e)), but the onset point was delayed compared with those cases shown in Figure 7(a)–7(c). When the water depth varied, the instantaneous vorticity revealed different spatial patterns. In Figure 7(f) and 7(g), although the swirling vortices were visible, the vortex shedding process hardly occurred. This wake pattern was similar to the near wake behind a solid cylinder, which was labelled as unsteady bubble (UB) wake by Chen & Jirka (1995). The bed friction is much stronger as the water depth decreases. The strong bed friction dissipated the flow energy and suppressed the spatial evolution of the separated shear layers originating at the patch shoulders. As a result, the formation of LKV was limited and the wake pattern changed. In contrast, the vorticity field for a sparse patch at different water depths did not change greatly. Since the patch-scale turbulence for a sparse patch was relatively weak, the vortex street was hardly generated.

To reveal the effect of the bed friction force on the patch-scale turbulence, the profiles of turbulence intensity in the vertical direction are displayed in Figure 8. In particular, the turbulent intensity was measured at the position where SSLs met (position of I_{TKE, max} in Figure 6). Figure 8(a) shows the turbulence intensity at low SVF, in which the cases without cylinders are taken as comparisons. The maximum turbulence intensities of the cases with bare bed were found to be lower than the cylinder-scale turbulence intensities I_stem. This explained that the turbulence within the cylinder array was dominated by wake generation, and it agreed with the observations of Nepf et al. (1997).

For low SVF, the bleed flow is strong, which weakens the formation and evolution of SSLs, and the turbulence intensity in the patch-scale wake is weak. Commonly, the turbulence intensity in the wake zone is weaker than cases with the bare bed. In free surface flows without lateral shear, the sizes of the turbulence structures are mainly governed by the balance of two forces. One force is the two-dimensional up-cascading of turbulent kinetic energy towards larger length scales at the free surface, and the other force comes from the interaction between large structures and the rough bottom (Uijttewaal & Tukker 1998). Since the bed friction is the only force dominating the turbulence evolution in cases of bare bed, the turbulence intensities monotonously decrease from the near-bed region to the water surface. The curve for P2H3 in Figure 8(a) also exhibited this trend, and it can explain why the overall level of TKE in the patch-scale wake was lower than the other cases in Figure 6(a). Because the bed shear effect is more and more strong as the water depth H decreases, the turbulence intensity was gradually weaker (shown in Figure 8(a) and 8(b)).

By contrast, the turbulence intensity behind a dense patch was much stronger. The suppression of patch-scale turbulence was not distinctively observed until the water depth reduced to under 0.04 m (see P1H2 in Figure 8(b)). As the turbulence generated by SSLs was highly suppressed by the bed friction at small water depths, the large-scale Kármán vortex street was gradually stabilized by the enhanced bed friction force. Additionally, the values of I_TKE measured at middle-depth z/H = 0.5 for cases P1H5, P1H6 and P1H7 were almost the same. The tiny differences between the results of these three cases were also observed in the mean flow and instantaneous vorticity field.

In the steady wake of the cylinder patch, the SSLs were formed and enhanced downstream along the two sides of the wake region until they developed to the centreline of the wake region. Beyond the interaction position, the patch-scale vortex street formed. The length of the steady wake region, L₁, was therefore determined by the scale of the patch (D) and the growth rate of the SSLs. As shown in Figure 9, the characteristic width of the shear layer is δ, and a new width, δ₁, instead of δ, was used to calculate the growth rate of δ₁, i.e., dδ₁/dx.

δ₁ was defined as the distance from the shoulder of the array, i.e., y/D = 0.5, to the position where the mean velocity began to increase from U₁ to U₂. The values of δ₁ were obtained by measuring the lateral distribution of the mean streamwise velocity. Because of the porous obstruction blocking the flow, the velocity U₂ is always larger than the velocity U₁. The shear layer flows generated by the deficit between U₁ and U₂ developed laterally downstream until the two SSLs merged at the centreline of the wake zone. When the two SSLs merged at the centreline, the steady wake region was closed and the width δ₁ reached maximum, i.e., δ₁ = D/2. Beyond the steady wake region, the velocity at the centreline increased because of the strong exchange of the flow momentum. As a result of the interaction of SSLs, the flow speed recovered gradually to the upstream free-stream velocity U.

Figure 10 shows the dependence of dδ₁/dx on the stability parameter S. The shear layer growth rate for low SVF (φ = 0.03) was less than that for high SVF (φ = 0.10), which explained a longer steady wake region downstream of a sparse cylinder patch. The shear layer growth rate decreased as S increased, which indicated the transverse spreading rate of the shear layer decreased when the bed friction became the primary dynamic force on the turbulent flow (Chu & Babarutsi 1988).

In this paper, shallow water flows past a circular array of emergent cylinders have been numerically investigated. Ten study cases with different solid volume fractions of the patch and water depths were simulated, and the investigations focused on the correlations of the wake flow pattern with the SVF, water depth, bottom friction and the stability dynamic factor. The topic can help researchers investigate the influence of vegetation on patterns of deposition and erosion in shallow water flows.

The wake structure behind a circular array of cylinders was not only determined by the SVF, but also dependent on the bed friction. Two typical types of wake, the vortex street (VS) wake and the unsteady bubble (UB) wake were both observed in the simulations. For a solid cylinder, the VS–UB transition was characterized by a wake stability parameter, S_c ≈ 0.2 (Chen & Jirka 1995). However, Chen & Jirka (1995) also observed the unsteady bubble wake behind a porous plate for S= 0.073. This means that the critical wake stability parameter, S_c, for a porous patch may be different from a solid body. Based on the present study cases, the critical S was preliminarily found to be 0.02 for φ = 0.10.

The bed friction played a different role in the cylinder-scale and patch-scale turbulent structures. Within the patch, even for the case of the sparse patch, the production of cylinder-scale turbulence was higher than the production by the bed shear. However, the patch-scale turbulence behind the patch can be greatly influenced by the bed friction. The peak of patch-scale turbulence occurred at the position where SSLs met and the shear layer instability took place. The von Kármán vortex street may also form at this position, but it depended on the level of the patch-scale turbulence intensity. The VS was suppressed by the bed friction force, which happened when S > 0.02 for φ = 0.10.

The lateral growth of SSLs was slower when the bed friction increased or the SVF decreased. The length of the steady wake, L₁, was related to the growth rate of SSLs and the patch diameter D. Therefore, the L₁ obtained from the simulations were compared with an empirical relationship. It was noted that the empirical relationship failed to predict L₁ at S > 0.03 with φ = 0.10 and S > 0.01 with φ = 0.03. The assumption of linear shear layer growth and a constant empirical parameter S_δ1 possibly contributed to the failure of the empirical formula.

A critical parameter, S_c, was estimated to be 0.02 for high SVF (φ = 0.10), and the S_c decreased as φ decreased. In order to extend a robust correlation between S_c and φ, a series of study cases with different SVFs must be investigated.

This work was sponsored by the National Natural Science Foundation of China (No. 11572196) and Shanghai Science and Technology Committee (No. 17230741200).