Influence of submerged flexible vegetation on turbulence in an open-channel flow

Abstract

In this study, we present a three-dimensional numerical model for the interaction of flow with submerged flexible vegetation, based on a large-eddy simulation and the immersed boundary method. The model innovatively realises the interaction between the flow and highly flexible vegetation with clustered leaves. Besides being a three-dimensional model of motion with full degrees of freedom, this study improves the consideration of the motion of the vegetation in all directions, and in addition the energy and momentum transfer in the spanwise direction. Furthermore, we perform a flume experiment for the flow with submerged flexible vegetation, the results of which are used to validate the simulation effects of the numerical model. It is found that the numerical model can effectively simulate the velocity profiles and the movement of vegetation induced by the flow. Using the model to analyse the flow–vegetation interaction, we find that the movement of vegetation is closely related to the flow velocity. As the flow velocity increases, both the offset angle and the vegetation swaying amplitude increase. Compared to vertical rigid vegetation, the tilting of flexible vegetation does not significantly change the velocity difference and the magnitude of the turbulent kinetic energy between the inside and the outside of the vegetation canopy, but it does weaken the disturbance to flow, thus reducing the resistance to flow. However, the swaying of vegetation dose significantly increase the velocity difference between the inside and the outside of the canopy. It forms Kelvin–Helmholtz hairpin vortices intensifying the turbulence production, and enhancing the disturbance and resistance to flow.

1. Introduction

Natural rivers and lakes are complex three-dimensional (3-D) ecosystems with different types of vegetation at different levels. The existence of vegetation can not only slow down the flow velocity and stabilise the riverbed, but also play an important ecological role purifying the water, beautifying the environment and providing the habitat for invertebrates, fish and other aquatic organisms (Nepf 2012; Shan et al. 2020). Therefore, it is of great significance to study the flow–vegetation interaction for analysing the turbulence structure and material transport from the perspective of protecting the aquatic ecological environment. So far, a large number of studies have been carried out on the macroscopic characteristics of vegetation flow, both experimentally (Ikeda & Kanazawa 1996; Nepf 1999; Wilson et al. 2003; Wang & Wang 2010) and numerically (Li & Zhang 2010; Zhang et al. 2017; Xiang et al. 2020).

As compared to the numerical model, the plant models used in flume experiments can be closer to the real shape of plants. Besides, the physical models are easier to set up since there is no need to analyse the complex relationship between the stress on plants and their deformation. Some researchers have studied the effects of aquatic plants on the flow velocity distribution in flume experiments. They found that submerged plants form a mixing layer in the flow leading to the velocity forming a profile similar to a hyperbolic tangent curve (Ghisalberti & Nepf 2002, 2004). The flow–vegetation interaction forms coherent motions, causing plants to form a coherent wave phenomenon, called monami, and enhancing the vertical transport of momentum (Ghisalberti & Nepf 2009; Okamoto & Nezu 2009). A number of researchers have also used flume and field experiments with plants of different shapes to study the impact of plant shape on the flow field. For instance, the effect of leaf area on the flow was analysed to reveal that the existence of leaves increases the flow velocity gradient between the inside and the outside of the vegetation canopy (Sand-Jensen & Mebus 1996; Hendriks et al. 2008; Wang et al. 2015). Nepf (2012) compiled a comprehensive review of studies on the influence of vegetation distribution density on the velocity profile, turbulent flow characteristics and material transport, based on flume experimental results. It is pertinent to mention that the turbulence intensity of canopy-scale vortices in the upper and outer layers of vegetation is high, which is the main reason for the mass and momentum exchange between the upper and outer layers of vegetation (Ghisalberti & Nepf 2005), while the turbulence intensity of stem-scale vortices in the lower layer is low (Brown & Roshko 1974; Nepf et al. 2007).

Because of the limitations of the flume experimental method, such as experimental equipment, space, materials and the precision of the measurement method, the phenomena and physical mechanism of the flow–vegetation interaction cannot be fully revealed. Therefore, it is essential to use numerical models to study the flow–vegetation interaction. Numerical methods established in earlier studies equate the effects of plants to an increase in the drag force on the flow or the bed (Jordanova & James 2003; Dijkstra & Uittenbogaard 2010; Rominger, Lightbody & Nepf 2010). The most common approach is to model the action of plants into a stress model, in which the role of vegetation is simplified as the change of shear stress in the boundary layer or the mixing layer near the bed to realise its effect on the flow (Liu & Shen 2008; Liu et al. 2013; Luhar & Nepf 2016). Naot, Nezu & Nakagawa (1996) combined an empirical model with a 3-D turbulent algebraic stress model to simulate the resistance effect of plants on flow, and studied the physical significance and effective expressions of the new variables generated due to the addition of plants. Zhang et al. (2005) established a two-dimensional (2-D) k–ε turbulence numerical model under the action of vegetation, which can be widely applied to natural river courses and wetlands with vegetation. Wang & Wang (2011) evaluated the changes of aquatic plants at different velocities and obtained the drag force formula of flow containing emergent and submerged plants.

Numerical models that account for the effects of vegetation to increase drag force or resistance to flow are widely used in engineering calculations, because of their high computational efficiency and simple physical mechanism. However, since these models cannot simulate the real morphology of individual plants, they cannot simulate the turbulence structure, nor the movement and deformation of plants due to the flow–vegetation interaction. In recent years, in order to solve the problem, many studies have analysed the influence of plants on the flow by directly simulating the real morphology of plants. The direct simulation method is more convenient to study the physical mechanism of the flow–vegetation interaction (Stoesser et al. 2009; Maza, Lara & Losada 2015; Boothroyd et al. 2016; Wolski & Tymiński 2020), and mass and momentum transport (Mayaud, Wiggs & Bailey 2016; Kim, Kimura & Park 2018; Liu et al. 2018). Presently, most of these models are 3-D rigid vegetation models, in which rigid cylinders were used to simulate the effects of plants on the flow (Stoesser, Kim & Diplas 2010; Huai, Xue & Qian 2015; Etminan, Lowe & Ghisalberti 2017). Neary et al. (2012) established a numerical model of rigid emergent plants based on the large-eddy simulation (LES) method to study the influence of the stems of emergent plants on the turbulence characteristics and sediment transport. They argued that the effect of vegetation is to enhance the stability of the riverbed or form a non-uniform bed morphology. Monti, Omidyeganeh & Pinelli (2019) established a numerical model of open-channel flow including submerged rigid plants based on the LES method and immersed boundary method (IBM). They found that, in the canopy area, the velocity is related to the local bed shear stress.

As compared to the two types of vegetation models mentioned above, there are still a few other numerical models that could directly simulate the flexible vegetation, and most of them are 2-D models (Zeller et al. 2014; Leclercq & de Langre 2016). The reason for this is accredited to the fact that the physical characteristics of flexible vegetation are difficult to describe by numerical models, and the flow–vegetation interaction greatly increases the computational complexity and difficulty of programming. Favier, Revell & Pinelli (2014) established a 2-D flexible flow–vegetation interaction model based on the lattice Boltzmann method and IBM. They studied the physical mechanism of the two-way interaction between incompressible oscillating free-surface flow and flexible flaps, and clarified the movement rules of the plants in the flow. Based on this, O'Connor & Revell (2019) found that the coherent fluctuations observed in vegetation actually couple the responses of the flow to the array of plants, rather than being a purely flow-driven instability. The vegetation movement is not only related to the flow, but also affected by their own natural frequency.

Besides, despite there having been a few attempts to develop a 3-D flexible plant model, they remain in their infancy. Tschisgale et al. (2021) established a sheet flexible plant model to simulate seagrass movement in a flow. The model could satisfactorily simulate the coherent fluctuations of plants. However, the shortcoming of the model is that the movement of vegetation was limited to the xz plane, without any movement in the spanwise direction. Marjoribanks et al. (2014) established an N-pendula plant model, which can be used to simulate single-stem plants with high flexibility. However, the model could not simulate flexible plants with complex morphology, and the accuracy of the model needs to be verified when the plant deformation is large.

Based on the existing numerical model of flexible vegetation flow, it can be found that there are still some shortcomings in current studies of flexible flow–vegetation interaction. Most numerical models are 2-D, and there are errors in simulating the flow and vegetation movement. In a few studies of 3-D flexible vegetation models, the plants were usually simplified to column-like cantilever structures (Marjoribanks et al. 2014; Wang et al. 2020), i.e. elastic rods with one end fixed. However, this kind of structure cannot accurately simulate the characteristic motion of highly flexible plants that can undergo large deformation. There have also been studies to model the plants as flaps (Tschisgale et al. 2021), but such structures generally restrict vegetation movement to the streamwise directions and fail to describe the spanwise motion of vegetation. In addition, the shapes of these two kinds of plant models are still too simple to describe the interactions between leafed plants and flows. To fill the gaps, a 3-D flexible flow–vegetation interaction numerical model is develop based on the LES method and IBM. In this model, the plant adopts the structure of pellet–rope series, which can effectively simulate the vegetation movement characteristics with soft stems having large deformations. In addition, the modified model can effectively simulate the mechanical properties of clustered leaves, which can make up for the shortcomings of the simple shape of vegetation models in previous studies. Furthermore, the degrees of freedom (DOF) of the plant movement in the model are not restricted, which can effectively simulate the 3-D vegetation characteristic motion in the streamwise, spanwise and vertical directions, making the simulation results of the flow–vegetation interaction more accurate.

A flume experiment was designed to validate the results of the 3-D numerical model. It is found that the model can predict the velocity distribution covering the inside and outside of the vegetation canopy, the flow turbulence in the streamwise and spanwise directions, and the vegetation movement in three directions. The model is more reliable than the previous 2-D numerical models.

In addition, we also compare and analyse the simulation results of the model with those of the 3-D rigid vegetation model. We analyse the flow–vegetation interaction from two aspects: the influence of the flow on the state of the motion of flexible vegetation, and the influence of the vegetation tilt and deformation on the flow patterns. In the first aspect, the effects of the flow velocity on the vegetation offset angle and swaying amplitude are analysed. In the second aspect, the effects of the vegetation tilt and deformation on the flow velocity, vortex structure and turbulent kinetic energy (TKE) are examined. Finally, the influence of the deformability of vegetation on its resistance to the flow is analysed by comprehensively exploring the flow–vegetation interaction characteristics. In essence, the results of this study can effectively improve the theoretical results in terms of the effects of flexible vegetation on the velocity distribution and vortex structure.

2. Numerical method

In this study, LES is used as the flow solver. The direct-forcing IBM method based on the Cartesian coordinate system is used to simulate the geometric motion characteristics of the plants and to establish the mechanistic relationship between the flow and the vegetation (Yang & Stern 2015; Fadlun et al. 2000).

2.1. Flow solver: large-eddy simulation (LES)

The dimensionless LES equations are obtained by filtering of the 3-D incompressible Navier–Stokes equations. The continuity and the Navier–Stokes equations are as follows:

(2.1)\begin{gather}\frac{{\partial {{\bar{u}}_i}}}{{\partial {x_i}}} = 0,\end{gather}

(2.2)\begin{gather}\frac{{\partial {{\bar{u}}_i}}}{{\partial t}} + \frac{\partial }{{\partial {x_j}}}({\bar{u}_i}{\bar{u}_j}) ={-} \frac{{\partial \bar{p}}}{{\partial {x_i}}} + \upsilon \frac{{{\partial ^2}{{\bar{u}}_i}}}{{\partial {x_j}\partial {x_j}}} - \frac{{\partial {{\bar{\tau }}_{ij}}}}{{\partial {x_j}}},\end{gather}

where u_i and u_j are the ith and jth components of the instantaneous dimensionless velocity vector (i, j = 1, 2, 3), respectively, x_i is the spatial location vector in the ith direction, p is the dimensionless pressure, υ is the coefficient of kinematic viscosity of the fluid, which is the inverse of the Reynolds number in the program, and the overbar represents time averaging. The subgrid-scale (SGS) stress ${\bar{\tau }_{ij}}$ results from filtering of the nonlinear convective fluxes. This term reflects the influence of the SGS turbulence on the large-scale turbulence structures. The SGS stress ${\bar{\tau }_{ij}}$ is calculated from the eddy viscosity relationship as

(2.3)\begin{equation}{\tau _{ij}} = {\bar{u}_i}{\bar{u}_j} - \overline {{u_i}{u_j}} ={-} {\upsilon _{SGS}}\left( {\frac{{\partial {{\bar{u}}_i}}}{{\partial {x_j}}} + \frac{{\partial {{\bar{u}}_j}}}{{\partial {x_i}}}} \right) + \frac{1}{3}{\delta _{ij}}{\bar{\tau }_{kk}},\end{equation}

where δ_ij = 1 when i = j, and υ_SGS is the SGS viscosity, being computed from the dynamic SGS model proposed by Germano et al. (1991).

In this study, the second version of a code called the ‘large eddy simulation on curvilinear coordinates’ (LESOCC2), which was first developed at the Institute for Hydromechanics, Karlsruhe Institute of Technology, Germany (Breuer & Rodi 1994; Fröhlich & Rodi 2002), is used for the simulations. In this code, the governing equations were discretised by the finite-volume method on non-staggered curvilinear grids. The details of the discretisation of the LESOCC2 are available elsewhere (Fang et al. 2014, 2018).

2.2. Vegetation movement solver: immersed boundary method (IBM)

The ‘direct-forcing immersed boundary method’ is used to describe the fluid–solid interaction between the flow and the plants (Peskin 1972; Mohd-Yusof 1997; Fadlun et al. 2000). Each plant is modelled as a series of pellets composed of the IBM boundary grids (figure 1).

Figure 1. A section of the pellet that creates a single plant. Boundary grids are the grids filled with hatch marks.

The momentum equation of the boundary grids can be written as follows:

(2.4)\begin{equation}\frac{{\partial \boldsymbol{u}}}{{\partial t}} + (\boldsymbol{u} \boldsymbol{\cdot}{\nabla})\boldsymbol{u} ={-} {\nabla}p + (\upsilon + {\upsilon _{\textrm{SGS}}}){\nabla ^2}\boldsymbol{u} + \boldsymbol{f},\end{equation}

where u is the velocity vector of boundary grids, ${\nabla}^{2}$ is the Laplace operator and f is the forcing source term. The terms in bold represent vectors. Equation (2.4) can be discretised in time using the direct-forcing method as follows:

(2.5)\begin{equation}\frac{{{\boldsymbol{u}^{n + 1}} - {\boldsymbol{u}^n}}}{{\Delta t}} = \textrm{RH}{\textrm{S}^{n + {1 / 2}}} + {\boldsymbol{f}^{n + {1 / 2}}},\end{equation}

where RHS is the sum of the convective and viscous terms and the pressure gradient in (2.4), and uⁿ is the calculated velocity of the boundary grids at the nth time step. According to the force imposed on the boundary grids at the nth time step, their state of motion can be analysed and the velocity of the grids at the (n + 1)th time step v^{n +1} can be deduced (see § 2.3). According to (2.5), letting u^{n +1} = v^{n +1}, the body force f can be calculated as

(2.6)\begin{equation}{\boldsymbol{f}^{n + {1 / 2}}} ={-} \textrm{RH}{\textrm{S}^{n + {1 / 2}}}\textrm{ + }\frac{{{\boldsymbol{v}^{n + 1}} - {\boldsymbol{u}^n}}}{{\Delta t}}.\end{equation}

2.3. Movement and force analysis of plants

In the numerical model, a single plant is composed of solid sphere pellets connected with disembodied ropes (figure 2a). Each pellet is created by multiple IBM boundary grids. This structure can be used to simulate a submerged plant with a thin, soft stem and clumped leaves, such as Ceratophyllum demersum (figure 2b).

Figure 2. (a) Shape of the plants in the numerical model. The fine lines in the figure do not exist in the numerical model, but are realised by constraining the relative motion between the pellets. (b) C. demersum. By varying the number and diameter of the pellets contained in a single plant in (a), different shape characteristics of the plants can be simulated.

The pellets obstruct the flow, changing its momentum and receiving the reaction force of the flow. In addition, they are also subject to buoyancy, gravity and pulling force used to restrain the relative movement of the adjacent pellets. The rotation of an individual pellet is severely limited by the existence of the rope tensions. To be explicit, when the straight line of the tension does not pass through the centre of the pellet due to rotation, the torque generated by the tension dominates all the torques applied to the pellet. It causes the pellet to rotate again to the state where the tensions pass through the centre of the pellet. Therefore, a pellet cannot rotate to a large angle or at a high angular velocity. This is consistent with the actual observation in the physical model study described in § 3. Based on this, the influence of the torque is ignored in the force analysis of the pellet. By analysing the momentum of a single pellet, the following equation can be obtained:

(2.7)\begin{equation}\frac{{\textrm{d}{\boldsymbol{V}_l}}}{{\textrm{d}t}} = {\boldsymbol{F}_l} + {\boldsymbol{T}_l} - {\boldsymbol{T}_{l - 1}},\end{equation}

where V_l is the velocity vector of the lth pellet (from the top to the bottom), F_l is the vector sum of the force of flow on the lth pellet (that is, the force calculated in § 2.2), the buoyancy and the gravity of the pellet, and T_l is the tension of the lth rope (from the top to the bottom), taking the direction of the pulling force of the rope on the lth pellet as the positive direction. The force analysis of the lth pellet is shown in figure 3.

Figure 3. Dynamics of the lth pellet and its adjacent pellets. The blue arrows represent the force vectors, the red arrows represent the unit vectors in the direction of the ropes, and the green arrows represent the velocity vectors. The physical parameters represented by the notation in the figure are explained in § 2.3.

Since (2.7) is not closed, the dynamic state of a single pellet cannot be accurately solved, and the motion constraint equations of the adjacent pellets need to be added. According to the motion constraints of the adjacent pellets, the following equations can be obtained:

(2.8)\begin{gather}{\boldsymbol{V}_l}\boldsymbol{\cdot }{\boldsymbol{e}_{l - 1}} = {\boldsymbol{V}_{l - 1}}\boldsymbol{\cdot }{\boldsymbol{e}_{l - 1}},\end{gather}

(2.9)\begin{gather}\frac{{\textrm{d}{\boldsymbol{V}_l}}}{{\textrm{d}t}}\boldsymbol{\cdot }{\boldsymbol{e}_{l - 1}} = \frac{{\textrm{d}{\boldsymbol{V}_{l - 1}}}}{{\textrm{d}t}}\boldsymbol{\cdot }{\boldsymbol{e}_{l - 1}},\end{gather}

(2.10)\begin{gather}|{{\boldsymbol{T}_{l - 1}}} |= \frac{{{{(|{{\bf (}{\boldsymbol{V}_{l - 1}} - {\boldsymbol{V}_l}) \times {\boldsymbol{e}_{l - 1}}} |)}^2}}}{{{R_{l - 1}}}},\end{gather}

where e_{l −1} is the unit vector in the direction of the (l − 1)th rope, and R_{l −1} is the length of the (l − 1)th rope. According to (2.7)–(2.10), the forces and instantaneous acceleration of the pellet, a_l = dV_l/dt, can be obtained. Assuming that the force on the pellet remains constant within a time step, the equations of motion of the pellet can be discretised in time as follows:

(2.11)\begin{gather}\boldsymbol{V}_l^{n + 1} = \boldsymbol{V}_l^n + \Delta t\,\boldsymbol{a}_l^n,\end{gather}

(2.12)\begin{gather}\boldsymbol{x}_l^{n + 1} - \boldsymbol{x}_l^n = \Delta t\,\boldsymbol{V}_l^n + \frac{{\Delta {t^2}\,\boldsymbol{a}_l^n}}{2},\end{gather}

where $\boldsymbol{x}_l^n$ is the displacement of the lth pellet at the nth time step. According to (2.12), the position of the pellets at each time point can be obtained.

Because the model developed in this study cannot simulate the collision between clumped leaves, it is only suitable for plants with a single stem, and not for plants with bifurcate stems. Therefore, the present model can be improved as a future scope of this study involving the collision of leaves. Nevertheless, the model still has a wide adaptability, because the number and diameter of clumped leaves for each plant can be set to different values.

3. Experimental set-up

3.1. Physical model study

In order to verify the accuracy of the numerical model, a physical model study was designed to simulate flow with vegetation. The experimental results were compared with the numerical model results. The experiment was carried out in a long rectangular tilting flume with a length of 16 m and a width of 0.5 m having a bed slope of 0.0025. A flow stabiliser and circulation device were provided upstream and downstream of the flume, respectively. Previous studies have shown that the influence of submerged vegetation on the flow velocity is affected by the relative submergence of plants, that is, the ratio of flow depth H to vegetation height h (Nepf & Vivoni 2000). In natural rivers, most of the submerged plants are found in the range of shallow submergence (1 < H/h < 5) (Chambers & Kalff 1985; Duarte 1991). In order to simulate the actual situation in the experiment, we set the flow depth and the vegetation height as 0.2 and 0.1 m, respectively, creating a relative submergence H/h = 2. Each plant was made of five smooth solid polypropylene pellet spheres (mass density ρ = 0.97 × 10³ kg m⁻³ and diameter D_s = 10 mm) in series. Each pair of pellets were connected by a polyester thread having a diameter of 0.3 mm and a length of 10 mm. The simulated vegetation was located in the middle of the flume. There were five rows of plants in the streamwise direction and nine rows in the spanwise direction. A total of 45 plants were set up in the experiment, consisting of 225 pellets. The spacing between the individual plants was set to half the plant height as 5 cm, as was done in previous studies (Favier et al. 2017). In this study, the simulated aquatic plants represented by C. demersum could not survive in the fast flow due to their soft stems, because they generally exist in rivers and lakes with low flow velocities (Hilt et al. 2018). In the physical model study, the incoming flow discharge was set to 0.01 m³ s⁻¹ (average velocity $\bar{u} = 0.1\;\textrm{m}\;{\textrm{s}^{ - 1}}$ and Reynolds number $Re = \bar{u}H/\upsilon = 20\,000$, where υ is the coefficient of kinematic viscosity of water).

In the experiment, the swaying and displacement of the plants were captured by a Canon EOS 50D camera with an imaging frequency of 30 Hz and were measured by an image processing program from the swaying angle of each plant in the captured images. The flow was measured by an acoustic Doppler velocimeter (ADV) at four vertical distances, 0.025, 0.075, 0.125 and 0.175 m. The sampling frequency of the ADV was set as 100 Hz for the data collection, and the acoustic frequency was 10 MHz. The velocities at these locations corresponded to the flow velocities near the bed, inside and outside of the vegetation canopy, and near the free surface. In the streamwise direction, the flow velocities at six locations were measured to analyse the velocity changes before, during and after the passage of flow through the vegetation. In the spanwise direction, in order to prevent the ADV probe from colliding with the plants, the measuring points were located at the middle of each two rows of the plants. The set-up of the physical model study is shown in figure 4.

Figure 4. Top view of the flume of the physical model study showing the zone of simulated vegetation within the dashed box. Circles in the dashed box represent the individual plants.

3.2. Numerical simulation

Numerical simulations were carried out to simulate the flow–vegetation interaction in the vegetation-covered zone of the flume (figure 4). In order to avoid disturbances near the inlet boundary and so that the simulation of the flow–vegetation interaction remained unaffected, a flow buffer zone 0.2 m long was kept upstream of the vegetation zone. In addition, to study the flow characteristics after passing through the vegetation zone, a space of length 0.2 m was allowed downstream of the vegetation. The overall size of the calculation domain was 0.6 m (x) × 0.5 m (y) × 0.2 m (z). The calculation domain and the simulated vegetation in the numerical simulation are shown in figures 5(a) and 5(b). Previous studies showed that the grid scale of the LES model should be between the Kolmogorov scale η (also called the dissipative scale, l_d) and the energy-containing scale l_e (Kolmogorov 1941; Zhang et al. 2005). These scales can be obtained from the dimensional analysis as

(3.1)\begin{gather}{l_d} = \eta \sim {\left( {\frac{{{\upsilon^3}}}{\varepsilon }} \right)^{1/4}},\end{gather}

(3.2)\begin{gather}{l_e} \sim \frac{{{{u^{\prime}}^3}}}{\varepsilon },\end{gather}

where u′ is the root mean square of the velocity fluctuations and ε is the TKE dissipation rate, which can also be obtained using the dimensional analysis as

(3.3)\begin{equation}\varepsilon \sim \frac{{{u^4}}}{\upsilon },\end{equation}

where u is the time-averaged flow velocity. Using (3.1)–(3.3), it can be estimated that the side length of the grids should be between 10⁻³ and 10⁻⁵ m. Complying with this requirement, considering the computational efficiency and in view of the fact that the pellets constituting the plant should not contain too few IBM boundary grids, the side length of the grid was set to be 1 mm in this study.

Figure 5. (a) Computational set-up of flow with vegetation and (b) the xz plane of two single plants with grids used in the simulation.

The boundary conditions of the momentum equation were as follows. In the x direction, the Dirichlet boundary was used upstream, implying the flow takes place from the boundary grids in the calculation domain with a constant bulk channel velocity U_b = Q/WH = 0.1 m s⁻¹, where Q is the flow discharge, W is the channel width and H is the flow depth. In addition, in order to study the relationship between the flow velocity and the plant movement, five additional numerical simulations were carried out with the bulk flow velocity U_b = 0.05, 0.06, 0.07, 0.08 and 0.09 m s⁻¹. The corresponding Reynolds numbers ranged over 10 000–18 000. The wall function was adapted at the bottom and the spanwise direction (y direction) in order to match with the physical model study, and the Manning roughness coefficient of the channel was the same as that in the physical model study (set as n = 0.045 s m^−1/3). The free surface was assumed to be a rigid lid with a slip condition, which is usually used to simulate the free surface of flow with minute fluctuation.

According to the experimental set-up, the methods of the physical model study and the numerical simulation, and considering the tolerance of flexible vegetation to flow velocity, nine groups of experiments are conducted in this study, as shown in table 1. Among them, the RVMRE20(V) case is a rigid vegetation model of the numerical simulation in which the vegetation remains stationary in the initial vertical state under the flow condition of Re = 20 000, while the RVMRE20(I) case is another rigid vegetation model of the numerical simulation under the flow condition of Re = 20 000 in which the vegetation remains stationary in an inclined state. The RVMRE20(I) case is set up to compare with the FVMRE20 case to distinguish the effects of the vegetation tilt and swaying on the flow conditions. To ensure that the effect of vegetation tilt on flow is consistent in the two cases, the inclination angles of the plants in the RVMRE20(I) case are set according to the time-averaged values of their swaying angles, ${\bar{\theta }_x}$, simulated by the FVMRE20 case (calculated in § 4.2).

Table 1. List of the experiments conducted in this study.

In the following sections, the accuracy of the numerical model is verified and the flow–vegetation interaction is analysed in detail according to the results of the nine groups of experiments. Specific verification and analysis are carried out as follows. The results of the PMRE20 and FVMRE20 cases are compared to verify the simulation effects of the numerical model on the flow–vegetation interaction (§ 4). In terms of the effects of the flow on the vegetation movement, six groups of numerical simulations of the flexible vegetation model (FVMRE20–FVMRE10 cases) with different flow velocities are analysed (§ 5.1). In terms of the impact of the vegetation motion on the flow, the simulation results of the rigid vegetation model of the numerical simulations (RVMRE20(V) and RVMRE20(I) cases) and the flexible vegetation model of the numerical simulation (FVMRE20 case) are compared. The influence of vegetation deformability on the flow velocity field, turbulence structure and energy transmission is examined (§ 5.2). According to these, the difference of flow resistance between the flexible and the rigid vegetation cases can be further obtained. Based on the analysis of the above results, the law governing the variation of vegetation canopy height caused by variabilities of the flow conditions and the influence of the vegetation deformability on the drag force acting on the vegetation due to the flow are summarised (§ 5.3).

In order to facilitate the model validation and the analysis of the results in the following sections, 45 plants are sorted in ascending order according to the spanwise distance (that is, y coordinate) as the main condition, and the streamwise distance (that is, x coordinate) as the secondary condition. The simulated vegetation zone in the physical model study, the calculation domain of the numerical simulations (RVMRE20(V), RVMRE20(I) and FVMRE20–FVMRE10 cases) and their dimensionless coordinates are shown in figure 6(a–c).

Figure 6. (a) Elevation view of the simulated vegetation zone (dashed box in figure 4) in the physical model study. The hollow circles represent the pellets that simulate the plants and the squares with crosses at the centre represent the ADV measuring points. The x and z are the streamwise and vertical distances, respectively, taking the dimensionless coordinate as z/D_s with the initial position of the vegetation canopy as 0 (that is, the vertical distance of 10D_s from the bed). (b) Top view of the simulated vegetation zone in the physical model study. All markings are the same as in (a). The y is the spanwise direction, taking the dimensionless coordinate as y/D_s with the initial position of the middle row of plants along the width of the flume (plants numbers 21–25) as 0. (c) Elevation view of the numerical simulation calculation domain. In the numerical simulation and the physical model study, the coordinate axes and directions are the same, and their values correspond to each other.

4. Model validation

The most important factors affecting the interaction between the flow and the flexible vegetation are the Kelvin–Helmholtz (KH) instability and the resulting vegetation movement. The main factor affecting the intensity of the KH instability is the flow velocity gradient at the flow–vegetation interface. Therefore, whether the numerical model can accurately predict the velocity difference between the inside and the outside of the vegetation canopy and accurately simulate the vegetation movement process proves to be of the utmost importance. Owing to the ADV-probe-induced flow disturbance, it is not appropriate to compare the difference between the physical model and the numerical model for the simulation of physical quantities associated with velocity fluctuation in this study. In this section, the simulation effects of the PMRE20 and FVMRE20 cases on the velocity difference inside and outside the vegetation canopy and the vegetation movement are mainly compared.

4.1. Flow velocity

Figure 7 shows the time-averaged streamwise velocity ${\bar{u}_x}$ of flow at different locations measured in the PMRE20 case and the computed values of ${\bar{u}_x}/{U_b}$ at the corresponding positions in the FVMRE20 case. The results show that, compared with the physical model study, the prediction error of this numerical model is less than 5%. According to the results of the FVMRE20 case, the numerical model has a good simulation effect on the flow velocity in the low-velocity zone below the vegetation canopy height and in the high-velocity zone above the vegetation canopy height.

Figure 7. Comparison of the dimensionless time-averaged streamwise velocity ${\bar{u}_x}/{U_b}$ of flow at all 192 ADV measuring points in the PMRE20 case and the computed value of ${\bar{u}_x}/{U_b}$ at the corresponding positions in the FVMRE20 case. The solid line is the y = x auxiliary line, representing the zero error. The two dotted lines are y = 1.05x and y = 0.95x auxiliary lines, representing 5% error.

In order to verify whether the numerical model can effectively calculate the difference of the flow velocity between the inside and outside of the vegetation canopy, figure 8 compares the time-averaged streamwise velocity ${\bar{u}_x}$ of flow obtained from the PMRE20 and FVMRE20 cases at 24 measuring points on the y/D_s = 2.5 plane for different x/D_s values. As compared with the PMRE20 case, the prediction errors of the FVMRE20 case of ${\bar{u}_x}/{U_b}$ above the vegetation canopy height are 0.76%, 0.34%, 0.29%, 5.38%, 1.12% and 4.63%, respectively. The prediction errors of the flow velocity below the canopy height are 3.73%, 9.51%, 9.49%, 3.66%, 10.65% and −6.55%, respectively. The calculation error of the FVMRE20 case for the difference between ${\bar{u}_x}$ inside and outside of the vegetation canopy is less than 0.1U_b, that is, the maximum error is less than 10% of the bulk channel velocity.

Figure 8. (a–f) Vertical distributions of the dimensionless time-averaged streamwise velocity ${\bar{u}_x}/{U_b}$ on the y/D_s = 2.5 plane for different x/D_s. The triangle symbols represent the measured results of the physical model study. The solid lines represent the simulation results of the FVMRE20 case.

The velocity distribution characteristics of vegetated flow have been studied extensively. Among them, the most representative one is Ghisalberti & Nepf (2002), who gave the velocity distribution in the mixing layer near the vegetation canopy height under different flow conditions. Their study showed that the dimensionless flow velocity $\hat{U} = (U - \bar{U})/\Delta U$ presents a uniform hyperbolic tangent distribution in the mixing layer (figure 9), where U is the time-averaged flow velocity, $\bar{U} = ({U_1} + {U_2})/2$, with U ₁ and U ₂ being the low and high stream velocities in the mixing layer, respectively. Since the mixing layer thickness is different under different flow conditions, the vertical coordinate should also be made dimensionless to facilitate comparison, namely, $\hat{z} = (z - \bar{z})/\theta$, where z is the initial vertical coordinate, and $\bar{z}$ is the vertical coordinate, where $U = \bar{U}$, and θ is defined as

(4.1)\begin{equation}\theta = \int_{ - \infty }^\infty {\left[ {\frac{1}{4} - {{\left( {\frac{{U - \bar{U}}}{{\Delta U}}} \right)}^2}} \right]\textrm{d}z} .\end{equation}

Figure 9. Time-averaged flow velocity distributions for the FVMRE20–FVMRE10 cases. The comparison between the simulated velocity profiles and the hyperbolic tangent flow velocity curve (Ghisalberti & Nepf 2002) is satisfactory.

Comparing the dimensionless flow velocity distribution of the FVMRE20–FVMRE10 cases in the mixing layer with the hyperbolic tangent flow velocity curve obtained by Ghisalberti & Nepf (2002), the simulation effect of this study on the flow velocity distribution with vegetation under different flow conditions can be verified (figure 9).

In figure 9 under various flow conditions (FVMRE20–FVMRE10 cases), the velocity distribution in the mixing layer is consistent with the hyperbolic tangent velocity curve in Ghisalberti & Nepf (2002). It can be seen that the velocity distribution simulated by the model in this study corresponds with the prediction results of the velocity distribution of vegetation flow in previous studies.

4.2. Vegetation movement

Since the upstream flow boundary has a constant discharge in the numerical simulation, the movement of vegetation is not due to the change of discharge, but due to the velocity fluctuation. The velocity fluctuation is highly random. The fluctuations that result from any two physical experiments and simulation cases cannot be exactly the same. Therefore, the physical model study and the numerical simulation are not able to get exactly the same simulation results of vegetation movement. Hence, we only compare the movement features of vegetation simulated in the FVMRE20 case with those obtained from the physical model study in a period of time from the perspective of statistics (figure 10a,b).

Figure 10. (a) Plant movement at a certain moment simulated in the FVMRE20 case. In order to avoid visual confusion caused by the overlapping projections of multiple rows of plants, only the plants on the y = 0 plane are shown. (b) Plant movement at a certain moment obtained from the PMRE20 case. The coordinates of the plants are shown in figure 4(b).

Figure 11(a,c,e,g,i) shows the dimensionless streamwise offset Δx/D_s of the top of plants numbers 21–25 simulated in the FVMRE20 case and obtained from the physical model study. It can be observed that the average values of Δx/D_s simulated in the FVMRE20 case are consistent with those obtained from the physical model study within an error of range 10%. Since the average values of Δx/D_s can approximately represent the positions of the equilibrium force of the pellets in the average flow velocity field, the numerical model can be accurate and effective in simulating the force on the pellets. Comparing the variation range of Δx/D_s simulated in the FVMRE20 case with that obtained from the physical model study, it can be found that the variation range of Δx/D_s of the five plants simulated in the FVMRE20 case is generally larger than that obtained from the physical model study. This is accredited to the measurement frequency of the plant movement trajectory in the FVMRE20 case (f_{FVMRE 20} = 1000 Hz), which is much higher than that in the physical model study (f_PM = 1 Hz). This makes it more likely to capture the extreme values of Δx/D_s in the FVMRE20 case.

Figure 11. (a,c,e,g,i) Fluctuations of the dimensionless streamwise offset Δx/D_s of the top of plants numbers 21–25 over a period of dimensionless time 0 < t/(D_s/U_b) < 400 simulated in the FVMRE20 case and obtained from the physical model study. The blue dotted lines and red triangles represent the simulation results of Δx/D_s in the FVMRE20 case and the data plots of the physical model study, respectively. The blue and red solid lines are the mean values of Δx/D_s simulated in the FVMRE20 case and obtained from the physical model study, respectively. (b,d,f,h,j) Probability density distributions PD of the fluctuations of Δx/D_s of the top of plants numbers 21–25 simulated in the FVMRE20 case and obtained from the physical model study. The solid and dashed lines represent the PD curves of the fluctuations of Δx/D_s simulated in the FVMRE20 case and obtained from the physical model study, respectively.

Besides, figure 11(b,d,f,h,j) presents the probability density (PD) distributions of the fluctuations of Δx/D_s of the top of plants numbers 21–25 simulated in the FVMRE20 case and obtained from the PMRE20 case. As can be seen, for the first four plants (numbers 21–24), Δx/D_s obtained from the PMRE20 case is slightly smaller than that obtained from the FVMRE20 case, while for plant number 25 the PD curves obtained in the two cases are quite consistent. According to the PD curves, it can be calculated that, compared with the PMRE20 case, the errors of the mathematical expectation of Δx/D_s of plants numbers 21–25 in the FVMRE20 case are 4.86%, 8.33%, 5.29%, 5.34% and −2.37%, respectively. The errors of the median values are 4.60%, 8.55%, 6.77%, 6.96% and −1.02%, respectively.

Figure 12 compares of the distributions of the streamwise offset angle θ_x of plants numbers 21–25 relative to their initial positions simulated in the FVMRE20 case and obtained from the physical model study over a period of dimensionless time 0 < t/(D_s/U_b) < 400. It is evident that, in the θ_x distribution, the median and the mean values are almost the same in both the FVMRE20 simulation case and the physical model study. In addition, the distances between the median value and the upper and lower quartiles are almost the same, as are the distances between the maximum value and upper quartile, and between the minimum value and lower quartile. These indicate that, under the condition of constant unidirectional flow, the plants sway uniformly on both sides of the force balance point with a similar amplitude. The probabilities of swaying to both the sides (including upstream and downstream) are mostly the same. Comparing the simulation results of the FVMRE20 case with the results obtained from the physical model study, trends similar to those in figure 11 can be obtained. It indicates that, except for plant number 21, the mean values of θ_x of the other plants are quite similar in the two experiments (numerical and physical), with a maximum error of 9.24%, while the range of extreme values is slightly smaller in the physical model study, with a maximum error of 14.76%.

Figure 12. Comparisons of the distributions of the streamwise offset angles θ_x of plants numbers 21–25 simulated in the FVMRE20 case (the blue rectangular boxes) and obtained from the physical model (PM) (the orange rectangular boxes). The crosses in the boxes represent the median values θ_xm, the horizontal lines in the middle of the boxes represent the average value ${\bar{\theta }_x}$ and the small solid circles represent the statistically significant outliers. Since all the data were simulated in the FVMRE20 case and obtained from the physical model, the outliers are reliable values and have some significance in indicating the maximum swaying amplitude of the plants. The average value ${\bar{\theta }_x}$ in the box plot can represent the offset angle corresponding to the plant force balance position. The difference between the upper and lower edges of the straight line represents the maximum swaying amplitude Δθ_x of the plant.

According to the comparison of multiple indices, the results of the FVMRE20 simulation case and the physical model study are in good agreement. This numerical model can effectively simulate the movement of flexible plants under the flow–vegetation interaction.

5. Numerical results

In order to comprehensively analyse the flow pattern of various places, according to the relative positions of flow and vegetation, we select several characteristic locations in the spanwise and vertical directions to conduct flow pattern analysis. In the spanwise direction, we mainly analyse flow patterns on the plane where the middle row of plants is located (y/D_s = 0) and the plane between the plants adjacent to it (y/D_s = 2.5). In the vertical direction, we mainly compare the flow patterns inside the vegetation (z/D_s = −5 in the rigid vegetation model, and lower in the flexible vegetation model), at the flow–vegetation interface (z/D_s = 0 in the rigid vegetation model, and a value selected according to the actual location of the interface in the flexible vegetation model) and on the plane outside the vegetation (z/D_s = 5). These characteristic locations are represented by abbreviations given in table 2.

Table 2. Coordinates of characteristic locations (and planes) and their abbreviations.

5.1. Impact of flow conditions on vegetation movement

The flow velocity variation directly affects the forces acting on the plants – see § 2.3 and Luhar & Nepf (2011) – and also the state of the vegetation movement. With an increase in flow velocity, the turbulence intensity increases and vortex structure in the flow changes, which in turn affect the amplitude of the swaying of the plants. In order to study the difference of the states of the plant movement under different flow velocities, we compare the statistical distributional features of the streamwise and spanwise offset angles, θ_x and θ_y, of the plants simulated by the six types of numerical simulations (FVMRE20–FVMRE10 cases) in figures 13(a) and 13(b).

Figure 13. Movement of five plants (numbers 21–25) under different flow velocities. Distributions of (a) the streamwise swaying amplitude θ_x and (b) the spanwise swaying amplitude θ_y of the top of the plants. The six colours represent six cases with Reynolds numbers varying from 10 000 to 20 000, corresponding to the FVMRE10–FVMRE20 cases. The symbols used in the figure are same as in figure 12.

It is apparent from figures 13(a) and 13(b) that, as the flow velocity increases, the streamwise force balance angles θ_x of the five plants increase, that is, the plants need to tilt to a larger angle to balance the increased horizontal force (Wilson et al. 2003), while the spanwise force balance angles θ_y remain unchanged. Although the streamwise average force on the plants increases with the flow velocity (Luhar & Nepf 2011), it has little effect on the spanwise average force. The maximum swaying amplitudes in the streamwise Δθ_x and the spanwise Δθ_y directions of the five plants increase clearly with an increase in flow velocity. This is attributed to the fact that, as the flow velocity increases, the turbulence in flow is intensified, and thus the maximum swaying amplitude of the plants enhances.

Under the same flow velocity condition, the streamwise force balance angles θ_x of the five plants decrease with an increase in their initial x coordinates, consistent with the results of the existing 2-D flexible vegetation models (Favier et al. 2017; O'Connor & Revell 2019). However, the maximum swaying amplitudes in the streamwise Δθ_x and the spanwise Δθ_y directions of the downstream plants are not significantly different from those of the upstream plants. This indicates that, in the vegetation zone, the average force on the downstream plants is significantly weaker than that on upstream plants at various flow velocities. On the other hand, there is no significant difference in turbulence intensities between the upstream and the downstream vegetation zones in small vegetation patches.

Analysing the relation of the flow velocity $\bar{u}$, the position of the average force balance angle $\bar{\theta }$ of plants and the maximum swaying amplitude Δθ, the relation of the flow velocity, stress on the plants and turbulence intensity in the vegetation zone can be obtained. Luhar & Nepf (2011) obtained the analytical solution of the streamwise offset angle ${\bar{\theta }_x}$ by analysing the relationship between the stress and the deformation of flexible submerged vegetation.

According to the force balance analysis of the buoyancy, drag force and plant resistance to flow in Luhar & Nepf (2011), a simplified formula for calculating the vegetation streamwise offset angle ${\bar{\theta }_x}$ can be obtained:

(5.1)\begin{equation}\textrm{sin}\,{\bar{\theta }_x} = CU_b^2\ \textrm{co}{\textrm{s}^2}\,{\bar{\theta }_x},\end{equation}

where C is a value related to the plant's own properties (such as stem diameter, plant length, etc.) and the drag force coefficient C_D on flow (Blevins 1984). In this study, the value of C is related only to the location of plants. For the same plant, the value of C does not vary with Re. Therefore, the simulation results of the FVMRE10 case are used to determine the value of C for plants numbers 21–25, and the theoretical values of vegetation offset angles under other Re conditions are calculated from (5.1). The theoretical values are compared with the simulation results of FVMRE12–FVMRE20.

Figure 14(a) shows the variations of the streamwise offset angle ${\bar{\theta }_x}$ and the maximum swaying amplitude Δθ_x of plants numbers 21–25 with Reynolds number Re, and compares the results with the predicted results of ${\bar{\theta }_x}$ obtained by Luhar & Nepf (2011). It can be seen that, for plant number 21, the experimental results of this study match well with those predicted by Luhar & Nepf (2011) at various Reynolds numbers. However, for plants numbers 22–25, the matching is found to be good only when Re is small (Re < 18 000). When Re = 20 000, the streamwise offset angle ${\bar{\theta }_x}$ of this study is greater than that predicted by Luhar & Nepf (2011). This is attributed to the fact that, in Luhar & Nepf (2011), the drag force on the plant is directly related to the bulk flow velocity U_b in the channel. For plant number 21, the flow velocity at the plant location is fairly consistent with U_b. Therefore, the method of Luhar & Nepf (2011) can accurately predict the drag force on the plant and calculate the accurate value of ${\bar{\theta }_x}$ based on it. However, for downstream plants, due to the flow obstruction by the upstream plants, the flow velocity at the location of the plants reduces nonlinearly with Reynolds number. Therefore, it is erroneous to calculate the drag force on plants based on the bulk velocity U_b. In essence, the method of Luhar & Nepf (2011) is not suitable for the estimation of ${\bar{\theta }_x}$ of each plant in a plant cluster. Based on the results of the effects of plant clusters on the flow velocity reduction in this study (§ 5.2.1), the method of Luhar & Nepf (2011) can be improved to obtain an analytical solution of ${\bar{\theta }_x}$ for different plants in a plant cluster.

Figure 14. (a) Mean value of the streamwise offset angles ${\bar{\theta }_x}$ (force balance angle) of the five plants (numbers 21–25) as a function of Reynolds number Re and (b) maximum swaying amplitude Δθ_x of the five plants as a function of Reynolds number Re. Different symbols correspond to different plant numbers. The dotted lines in (a) are the predicted values of ${\bar{\theta }_x}$ obtained from the model of Luhar & Nepf (2011), and the dotted lines in (b) are the linear fitted lines.

Figure 14(b) shows the variations of the maximum swaying amplitude Δθ_x of plants numbers 21–25 with Reynolds number Re. It is evident that the maximum swaying amplitude Δθ_x of the plants increases linearly with Re. The gradients of the straight lines obtained by fitting linear curves are almost the same, indicating that the modes of the swaying amplitude Δθ_x of the upstream and downstream plants changing with Re are essentially the same. For a given Re, the values of the maximum Δθ_x of the five plants have no fixed relationship. The maximum Δθ_x difference between two adjacent plants is within 20°. This reveals that, at various Re, the turbulence intensity in the flow within the vegetation zone is roughly uniform and does not vary significantly in the downstream region in small vegetation patches (see § 5.2.2 for details).

5.2. Impact of vegetation movement on the flow structure

5.2.1. Averaged flow velocity field

The influence of the flexible vegetation movement on the flow velocity field can be divided into two parts. The first part is related to the change in flow velocity caused by the vegetation tilt under the action of flow impulse. It can be analysed by examining the simulation results obtained from the RVMRE20(V) and RVMRE20(I) cases. The other part is related to the disturbance of the flow caused by the swaying of vegetation around the average force balance position. It can also be analysed by inspecting the simulation results obtained from the RVMRE20(I) and FVMRE20 cases.

Figure 15(a–c) depicts the dimensionless streamwise time-averaged velocity structures ${\bar{u}_x}/{U_b}$ on the plane y/D_s = 0 (Y1) simulated numerically in the RVMRE20(V), RVMRE20(I) and FVMRE20 cases, respectively. Examining figure 15(a,b), it can be seen that the vegetation tilt eliminates the ribbon distribution of ${\bar{u}_x}$ within the vegetation zone and makes the velocity distribution within the vegetation more uniform. The canopy height of inclined vegetation decreases, but does not significantly change the ${\bar{u}_x}$ inside and outside the vegetation canopy. However, given that the effective volume occupied by the vegetation zone decreases, it can be predicted that the resistance of the vegetation to the flow decreases (specifically analysed in § 5.3). In essence, vegetation tilt does not significantly change the difference of ${\bar{u}_x}$ between the inside and the outside of the vegetation canopy, and therefore it does not significantly change the momentum and mass transport inside and outside of the vegetation canopy.

Figure 15. Dimensionless time-averaged streamwise velocity structures ${\bar{u}_x}/{U_b}$ on the Y1 plane simulated numerically in the (a) RVMRE20(V) (vertical rigid vegetation with Re = 20 000), (b) RVMRE20(I) (inclined rigid vegetation with Re = 20 000) and (c) FVMRE20 (flexible vegetation with Re = 20 000) cases. Owing to vegetation movement, the characteristic location Z2 (that is, the flow–vegetation interface) is not a fixed plane in (b and c).

Examining figure 15(b,c), it is evident that the swaying of vegetation further homogenises the distribution of ${\bar{u}_x}$ within the vegetation zone. At the same time, the average velocity ${\bar{u}_x}$ inside the vegetation decreases and ${\bar{u}_x}$ above the vegetation canopy increases, causing the difference between ${\bar{u}_x}$ inside and outside the vegetation canopy to enhance. Therefore, the swaying of the flexible vegetation around the force balance point is the main factor to intensify the energy and mass transport between the inside and outside of the vegetation canopy.

Figure 16 shows the vertical distributions of the double-averaged (spatially averaged in the spanwise direction) streamwise velocity ${\bar{u}_x}$ at various places in the vegetation zone, and upstream and downstream of the vegetation zone for different x/D_s. It is obvious that the flow velocity field is affected before the flow enters the vegetation zone, regardless of the vegetation being rigid or flexible. If the local banded distribution of velocity caused by plant shape is not considered, the velocity distributions obtained from the three simulation cases are similar at x/D_s = 17.5. The blocking effect of the vegetation on the flow can be reversed upstream, and it may be equivalent to the increase in the shear layer thickness near the bed surface.

Figure 16. (a–d) Vertical distributions of the dimensionless width-averaged streamwise velocity ${\bar{u}_x}/{U_b}$ within and downstream of the vegetation zone for different x/D_s. The solid lines represent the simulation results of the RVMRE20(I) case, the blue dotted lines represent the simulation results of the RVMRE20(V) case and the red broken lines represent the simulation results of the FVMRE20 case.

The flow velocity field inside the vegetation zone is more complex (figure 16b,c). It is noticeable that the velocity distributions computed from the RVMRE20(I) and FVMRE20 cases are consistent with those of Nepf (2012) and Ghisalberti & Nepf (2002). The effects of the vegetation on the flow can be equivalent to forming a mixing layer, in which the flow velocity ${\bar{u}_x}$ does not show a monotonic trend with the vertical distance. Unlike vertical rigid vegetation, the ${\bar{u}_x}$ with inclined rigid vegetation declines with an increase in the mixing layer thickness, but the ${\bar{u}_x}$ inside the mixing layer does not decrease significantly. Therefore, there is no significant change in velocity difference between the inside and the outside of the vegetation canopy. Besides, the vegetation tilt even tends to decrease the flow velocity above the canopy. It follows that the effect of vegetation tilt is to reduce the momentum loss of the flow compared to the vertical state (as predicted by Nepf (1999)).

However, the swaying of flexible vegetation on the basis of tilt significantly increases the mixing layer thickness, the flow velocity above the canopy, and the velocity gradient inside and outside the vegetation canopy, and reduces the flow velocity within the vegetation. This means that the effect of swaying of vegetation is to intensify the KH instability near the canopy, and the mass and energy transfer inside and outside of the vegetation canopy. The oscillation of vegetation significantly increases the momentum consumption of the flow. It can be predicted that the swaying of vegetation increases the drag force on the flow, and thus the obstruction to flow is more obvious.

After the fluid has flowed out of the vegetation zone, the influence of the vegetation gradually weakens, and the effect of the vegetation on the flow can be equivalent to changing the boundary layer characteristics. In the three simulation cases, the results of the FVMRE20 case show that the swaying of vegetation has a more lasting effect on the flow. This indicates that the effect of swaying of vegetation is to significantly intensify its disturbance on the flow (see §§ 5.2.2 and 5.2.3 for details).

5.2.2. Turbulence structure

Vorticity is the curl of the flow velocity vector, which is an important hydrodynamic parameter to characterise the intensity and direction of the swirl in a flow (Wallace & Foss 1995). The vorticity can be calculated as follows:

(5.2) \begin{equation} \boldsymbol{\omega } = {\nabla} \times \boldsymbol{u } = \left( {\frac{{\partial {u_z}}}{{\partial y}} - \frac{{\partial {u_y}}}{{\partial z}}} \right)\boldsymbol{i} + \left( {\frac{{\partial {u_x}}}{{\partial z}} - \frac{{\partial {u_z}}}{{\partial x}}} \right)\boldsymbol{j} + \left( {\frac{{\partial {u_y}}}{{\partial x}} - \frac{{\partial {u_x}}}{{\partial y}}} \right)\boldsymbol{k},\end{equation}

where ω is the vorticity vector, u, u_x, u_y and u_z are the velocity vector and its components in the x, y and z directions, respectively, and i, j and k are the unit vectors in the x, y and z directions, respectively.

Figure 17(a–c) depicts the dimensionless vorticity structures ω_z/ω ₀ on the Z1 plane simulated in the RVMRE20(V) (z/D_s = −5), RVMRE20(I) (z/D_s = −7.5) and FVMRE20 (z/D_s = −7.5) cases, respectively, where ω_z is the vorticity component about the z axis and ω ₀ = U_b/D_s. From the simulation results of the RVMRE20(V) case (figure 17a), it is obvious that there are two rows of alternately generated vortices downstream of a single plant, forming a commonly known wake structure (Nikora 2010), in which the turbulence only occurs directly downstream of plant bodies. The two rows of vortices rotate in opposite directions and form a vortex street. Within the vegetation zone, the two rows of vortices move downstream with a minimal spreading in the spanwise direction. Between any two adjacent rows of the plants, ω_z vanishes. Downstream of the vegetation zone, ω_z decreases rapidly in a short distance. Besides, the two rows of vortices in opposite directions mix with each other and dissipate rapidly. Then, the flow pattern returns to the initial stable state within a short distance. The vorticity distribution of flow in vegetation obtained from the RVMRE20(I) case is consistent with the results obtained from the RVMRE20(V) case, that is, the vorticity is distributed in a ribbon and does not diffuse in the spanwise direction. Unlike the simulation results of the RVM cases, the vortices inside the vegetation simulated in the FVMRE20 case are larger, but more broken (figure 17c). The vortices are irregularly directed, and the magnitudes of ω_z of the two groups of adjacent vortices with opposite directions are different. The vegetation movement forms a jet structure, that is, not only is the turbulence distributed downstream of the plant bodies, but also the momentum is transmitted to the surrounding water body. In the spanwise direction, the magnitudes of ω_z are uniformly distributed; however, there is no ribbon distribution of the vorticity, as observed in figures 16(b) and 17(a), similar to that simulated in the RVM cases. This indicates that the vegetation movement causes vortex diffusion and mixing in the spanwise direction, even in the gap between the two adjacent rows of the plants.

Figure 17. Dimensionless vorticity structures ω_z/ω ₀ simulated in (a) the RVMRE20(V) case, (b) the RVMRE20(I) case and (c) the FVMRE20 case on the Z1 plane.

Since the vorticity distributions in the vegetation zone obtained from the RVMRE20(V) and RVMRE20(I) cases are similar, this implies that the plant tilt does not affect the vortex structure of flow in the vegetation zone. Therefore, the effects of the flexible vegetation on the vortex structure of flow are mainly realised by the swaying. Figure 18 compares the vortex structure of flow obtained from the RVMRE20(I) and FVMRE20 cases. Figure 18(a) shows that, when there is no swaying of vegetation, the vortex formed is mainly a small hairpin vortex (HP vortex). The vortex size is similar to the plant diameter. The KH instability at the fluid–vegetation interface is not strong enough to form large-scale KH vortices. According to figure 18(b), when the vegetation sways, the vortex scale significantly increases, and KH vortices with the same scale as the plant spacing are formed at the fluid–vegetation interface. These KH vortices are connected with the HP vortices and form the KH–HP vortex structures, in conformity with Tschisgale et al. (2021). Unlike Tschisgale et al. (2021), the spanwise swaying of the vegetation is considered in this study, making the computation of the spanwise transfer of the TKE more accurate (see § 5.2.4). The spanwise transfer of the TKE significantly increases the difference of velocity and pressure pulsations in different sections in the spanwise direction. Therefore, the spanwise scale of the KH vortices computed in this study is no larger than the plant spacing, being significantly smaller than those computed by Tschisgale et al. (2021).

Figure 18. Vortex structure of flow in the vegetation zone obtained from (a) the RVMRE20(I) and (b) the FVMRE20 cases. The vortex structure in the figure is represented by the iso surface of pressure fluctuation. The colour of the background flow field indicates the instantaneous velocity. The colour of the pressure pulsation iso surface indicates the velocity fluctuation.

The distribution of the HP vortex structure and the KH–HP vortex structure and their influence on the flow can be directly displayed by quantitatively comparing the eddy size at different positions under the condition of swaying and stationary vegetation. The eddy size in the inertial subrange is defined as the Taylor microscale λ_T, which is given by

(5.3)\begin{equation}{\lambda _T} = {\left( {\frac{{15\upsilon \overline {u^{\prime}u^{\prime}} }}{\varepsilon }} \right)^{0.5}},\end{equation}

where ε is the TKE dissipation rate (the same as in (3.1)), which can be estimated by using Kolmogorov's second hypothesis (Dey et al. 2012):

(5.4)\begin{equation}k_w^{5/3}E({k_w}) = C{\varepsilon ^{{2 / 3}}},\end{equation}

where k_w is the wave number, E(k_w) is the energy spectrum function that represents the TKE k′ to contain eddies of size l having wave number k_w (= 2π/l) and C is a constant, which is usually set as 0.5 (Monin & Yaglom 2007). The energy spectrum function E(k_w) can also be regarded as the TKE k′ energy density of eddies with wave number k_w, which satisfies

(5.5)\begin{equation}k^{\prime} = \int_0^\infty {E({k_w})\ \textrm{d}{k_w}} .\end{equation}

Figure 19(a–d) compares the energy spectra E(k_w) of the FVMRE20 and the RVMRE20(I) cases at the height of the vegetation canopy directly above a row of plants (Y1) and in the middle of two rows of plants (Y2). It can be observed that Kolmogorov's −5/3 scaling law is satisfied in the inertial subrange of flow with vegetation for both the rigid and the flexible cases. It can therefore be concluded that the value of $k_w^{5/3}E({K_w})$ is relatively constant in the inertial subrange, being independent of k_w (as shown in figure 19e–h). Accordingly, (5.2) and (5.3) can be used to compare the size of the eddies at different positions for stationary and moving vegetation flows.

Figure 19. (a–d) Energy spectra E(k_w) and (e–h) estimations of the TKE dissipation rate ε in flow with stationary (RVMRE20(I) case) and moving (FVMRE20 case) vegetation. Here ${\tau ^{\prime}_{xx}}$ represents the streamwise Reynolds normal stress, i.e. ${\tau ^{\prime}_{xx}} = \overline {u^{\prime}u^{\prime}} /U_b^2$.

Examining figures 19(e) and 19(g), it is apparent that, near individual plants, the TKE dissipation rate ε in the FEVRE20 case is nearly the same as that in the RVMRE20(I) case, while the magnitude of the streamwise flow velocity fluctuation (i.e. streamwise Reynolds normal stress $\overline {u^{\prime}u^{\prime}}$) is about 3–4 times that of the RVMRE20(I) case. It can be calculated from (5.2) that the averaged eddy size λ_T of the vortices near the flexible plants is about 1.73–2 times that of the rigid vegetation. On the other hand, from figures 19(e) and 19(f), it is evident that the TKE dissipation rate ε at the gaps between two rows of plants in the FEVRE20 case is about 0.46 times that of the ε near the plants in the RVMRE20 case. Notably, the $\overline {u^{\prime}u^{\prime}} $ of the former is 2–3 times that of the latter. It can be calculated that the averaged eddy size λ_T of the vortices in the gaps between the flexible plants is about 2–2.6 times that of rigid vegetation. In figure 19(h), the magnitude of $\overline {u^{\prime}u^{\prime}} $ is very small in the gaps between the rigid plants, and there are no vortices with eddy size λ_T similar to those in figure 19(e–g).

From the analysis of the above phenomenon, it is apparent that the eddy size of the vortices in the flexible vegetation canopy is significantly larger than that in the rigid vegetation canopy. In the flow with rigid vegetation, vortices caused by the vegetation disturbance only exist near the plants and the wake region downstream of the plants, called the HP vortex structure; while in the flow with flexible vegetation, in addition to the HP vortex structure, there also exist larger-scale vortices widely distributed in the whole canopy, called the KH vortex structure. This is consistent with the intuitive result in figure 18. The existence of the KH vortex structure significantly increases the dissipation rate of TKE and intensifies the turbulence in the flow, accelerating the momentum transfer and the energy loss in the flow.

5.2.3. Turbulent kinetic energy (TKE) budget

The TKE represents the turbulence intensity of flow and is directly related to the momentum, energy and mass transport in the boundary layer and mixing layer. The time-averaged TKE k′ can be estimated as follows:

(5.6)\begin{equation}k^{\prime} = {\textstyle{1 \over 2}}(\overline {u^{\prime}u^{\prime}} + \overline {v^{\prime}v^{\prime}} + \overline {w^{\prime}w^{\prime}} ),\end{equation}

where u′, v′ and w′ are the velocity fluctuations in the x, y and z directions, respectively. The quantities with overbars are time-averaged quantities.

Figure 20(a–f) compares the dimensionless TKE structures $k\prime /U_b^2$ simulated in the RVMRE20(V), RVMRE20(I) and FVMRE20 cases on the Y1 and Y2 planes. It can be observed that, for vertical rigid vegetation, the TKE is mainly distributed in the plant spacing downstream of the individual plants. The range of influence of the turbulence caused by the plant obstruction to the flow is less than twice the plant diameter. The turbulence neither spreads in the spanwise direction (Neary et al. 2012; Monti et al. 2019) nor affects the zones above the canopy height. Unlike the case of vertical rigid vegetation, the TKE distribution in the flow with inclined rigid vegetation is more uniform (figure 20c,d). Vegetation as a group influences the turbulence in flow rather than through the obstruction of individual plants in it. Akin to vertical rigid vegetation, the TKE in the flow with inclined rigid vegetation is mainly distributed directly downstream of the individual plants. The TKE does not diffuse in the spanwise direction, but tends to diffuse above the vegetation canopy. The TKE in the flow with inclined rigid vegetation is obviously smaller than that with vertical rigid vegetation. This indicates that the vegetation tilt significantly reduces the disturbance and obstruction to the flow by reducing the canopy height and effective waterward area (Nepf 1999). The TKE distribution within the flexible vegetation is more uniform, and the influence of the individual plants is further weakened as compared to that within inclined rigid vegetation (figure 20e,f). A strong turbulent band is formed in the mixing layer. The trend of the TKE diffusion in the spanwise direction and above the canopy height is also stronger than that of the two kinds of rigid vegetation.

Figure 20. Dimensionless TKE structures $k^{\prime}/U_b^2$ simulated in (a) the RVMRE20(V) case on the Y1 plane with y/D_s = 0, (b) the RVMRE20(V) case on the Y2 plane with y/D_s = 2.5, (c) the RVMRE20(I) case on the Y1 plane, (d) the RVMRE20(I) case on the Y2 plane, (e) the FVMRE20 case on the Y1 plane and (f) the FVMRE20 case on the Y2 plane.

Figure 21 depicts the vertical distributions of the width-averaged TKE k′ and its variation along the streamwise direction within and downstream of the vegetation zone for different x/D_s. It can be seen that the distribution of the TKE in the downstream direction for vertical rigid vegetation is mostly unchanged, and all of the TKE is distributed below the vegetation canopy height. Within the vertical rigid vegetation, the peaks of the width-averaged TKE are in the range ${k^{\prime}_{max}} = (0.01-0.02)U_b^2$. The peak values of the width-averaged TKE in the flow with inclined rigid vegetation are similar to those with vertical rigid vegetation. However, since the equivalent spacing between the plants shrinks after the vegetation tilt, the disturbance of the upstream vegetation to the flow can be transmitted to the downstream plants. As a result, the TKE is more evenly distributed within vegetation. In addition, the TKE also diffuses to the outside of the vegetation canopy, but the diffusion distance is limited. To sum up, the disturbance of inclined rigid vegetation to the flow is weaker than that of vertical rigid vegetation. The peak values of the TKE in the flow with flexible vegetation are much higher than in those with rigid vegetation, reaching ${k^{\prime}_{max}} = (0.08-0.09)U_b^2$. In addition, the vegetation disturbance to the flow also spreads to a larger area outside the vegetation canopy. In conclusion, the vegetation tilt does not change the magnitude of the TKE numerically, but makes the distribution range of the TKE smaller, and weakens the disturbance of vegetation to the flow. However, the effect of swaying of the vegetation is to increase the TKE diffusion in the spanwise direction and towards/beyond the canopy. It also causes an enhancement of the peak value of TKE, and thus enhances the disturbance of the vegetation to the flow.

Figure 21. (a–f) Vertical distributions of the dimensionless width-averaged TKE $k^{\prime}/U_b^2$ within and downstream of the vegetation zone for different x/D_s. The meanings of the different types of lines in the figure are the same as in figure 16.

The analysis of the TKE budget terms, such as TKE production, diffusion, dissipation and pressure energy diffusion rates, in the flow with vegetation allows one to quantify the influence of the presence of the vegetation and its movement on the surrounding flow. The expression of the TKE budget was derived by Raupach, Antonio & Rajagopalan (1991) for the flow over canopies. The time-averaged TKE budget in an open-channel flow is given by

(5.7)\begin{equation}P + {T_t} + {T_p} - \varepsilon = 0,\end{equation}

with

(5.8a–c)\begin{equation}P ={-} \overline {{{u^{\prime}}_i}{{u^{\prime}}_z}} \frac{{\partial {{\bar{u}}_i}}}{{\partial z}},\quad {T_t} ={-} \frac{{\partial \overline {{u_z}k^{\prime}} }}{{\partial z}},\quad {T_p} ={-} \frac{1}{{{\rho _w}}}\;\frac{{\partial \overline {p^{\prime}{{u^{\prime}}_z}} }}{{\partial z}},\end{equation}

where P is the TKE production rate, T_t is the TKE diffusion rate, T_p is the pressure energy diffusion rate, p′ is the pressure fluctuation and ε is the TKE dissipation rate (the same meaning as ε in (3.3)). According to Kolmogorov's similarity hypothesis (Kolmogorov 1941), the ε value is closely related to the kinematic viscosity and the strain rate of the velocity fluctuation in a small-scale universal equilibrium range. The relationship can be expressed as follows:

(5.9)\begin{equation}\varepsilon = 2\upsilon \overline {{{s^{\prime}}_{ij}}{{s^{\prime}}_{ij}}} ,\end{equation}

with

(5.10)\begin{equation}{s^{\prime}_{ij}} = \frac{1}{2}\left( {\frac{{\partial {{u^{\prime}}_i}}}{{\partial {x_j}}} + \frac{{\partial {{u^{\prime}}_j}}}{{\partial {x_i}}}} \right).\end{equation}

Figure 22 shows the vertical distributions of the width-averaged TKE production rate and dissipation rate within the vegetation zone. Evidently the TKE does not have large-scale vertical diffusion within the vertical rigid vegetation. The TKE is produced only in the wake flow of the individual plants and dissipates downstream. The TKE production rate and dissipation rate in the flow with inclined rigid vegetation are basically balanced, and there is only a small range of vertical diffusion near the mixing layer (figure 22b). On the whole, the TKE production rate and dissipation rate of vegetation are not significantly affected by the vegetation tilt, but their distribution range is reduced. The distributions of the TKE production rate and dissipation rate within flexible vegetation are quite different from those within rigid vegetation (figure 22c). The TKE is mainly generated in the mixing layer near the canopy height, and diffuses vertically to the interior of the vegetation, and dissipates mainly below the mixing layer. In addition, due to the swaying of vegetation, the TKE production and dissipation rates increase significantly. In addition, the distribution range also increases. The above phenomenon indicates that the vegetation tilt does not significantly affect the energy transfer and dissipation within the vegetation zone, but does decrease the energy dissipation of the whole system due to the reduction of the equivalent volume occupied by the vegetation. However, the swaying of the vegetation does significantly increase the energy transfer and dissipation.

Figure 22. Vertical distributions of the dimensionless width-averaged TKE production rate $P/(U_b^3{D_s})$ (solid lines) and the TKE dissipation rate (dotted lines) within and downstream of the vegetation zone in (a) the RVMRE20(V) case, (b) the RVMRE20(I) case and (c) the FVMRE20 case. The black lines represent the simulation results on the x/D_s = 37.5 section, and the red lines represent the simulation results on the x/D_s = 42.5 section.

The difference of TKE dissipation rate ε shown in figure 22 in flows with moving and stationary vegetation is consistent with the trend in the energy spectrum function E(k_w), shown in figure 19. There is little difference in the TKE dissipation rate ε near the individual plants at canopy height, which is consistent with the results shown in figure 22(b,c). However, in the gap between plants, the TKE dissipation rate ε increases significantly due to the vegetation movement (figures 19f and 18h). Overall, the vegetation movement leads to a significant increase in the TKE dissipation rate ε. This is in agreement with the conclusion obtained by TKE budget analysis. The above phenomenon indicates that the vegetation movement leads to an additional flow energy loss.

5.3. Impact of deformability of vegetation on its resistance to flow

According to the comparison of the influences of flexible vegetation and rigid vegetation on the flow velocity, the turbulence structure and the TKE in § 5.2, it is apparent that vegetation tilt has little influence on the flow velocity difference between the inside and the outside of the vegetation canopy and the turbulence structure in the flow, but it has a significant influence in reducing the disturbance to the flow. The swaying of the vegetation increases the velocity difference between the inside and the outside of the vegetation canopy. It can aggravate the KH instability of the mixing layer, forming the KH–HP vortex structure and intensifying the energy and mass exchange between the inside and the outside of the vegetation zone. In addition, it significantly increases the disturbance of vegetation to the flow. Thus, it can be inferred that the vegetation tilt can reduce the resistance of vegetation to the flow (Nepf 1999), while the effect of the vegetation swaying is opposite in nature. In order to verify this inference, this section compares the energy loss of flow through the vegetation zone computed from the RVMRE20(V), RVMRE20(I) and FVMRE20 cases, and analyses the influence of the vegetation deformability on its resistance to the flow.

The total flow energy E_tot is defined as the sum of the kinetic energy E_k, gravitational potential energy E_g and pressure potential energy E_p of the flow. Therefore, the time-average total energy of the flow ${\bar{E}_{tot}}$ can be expressed as

(5.11)\begin{equation}{\bar{E}_{tot}} = {\bar{E}_k} + {\bar{E}_g} + {\bar{E}_p},\end{equation}

with

(5.12a–c)\begin{equation}{\bar{E}_k} = {\textstyle{1 \over 2}}\rho \overline {{u_i}{u_i}} = {\textstyle{1 \over 2}}\rho ({\bar{u}_i}{\bar{u}_i} + \overline {{{u^{\prime}}_i}{{u^{\prime}}_i}} ) = k + k^{\prime},\quad {\bar{E}_p} = \bar{p},\quad {\bar{E}_g} = \rho gz,\end{equation}

where g is the gravitational acceleration and z is the vertical distance.

Examining the average energy in the flow in each section from upstream to downstream before the flow enters the vegetation zone (in the process of flowing through the vegetation and flowing out of the vegetation), the energy loss in the flow through the vegetation can be obtained. Then, the difference of the obstruction effects of flexible vegetation and rigid vegetation on the flow can be compared.

Since the average gravitational potential energy of the fluid particles at each section is zero, it is not listed and analysed separately. The average total energy and its components simulated in the RVM and FVMRE20 cases at the seven sections x/D_s = 10–42.5 as well as the variation of the average total energy along the channel are listed in table 3.

Table 3. Average values of the components of the dimensionless total flow energy at seven sections from upstream to downstream. The $\Delta {\bar{E}_{tot}}$ values represent the difference between the average total energy of the current section and that at the section x/D_s = 10. The zero point of pressure potential energy is taken as the average pressure potential energy at the section x/D_s = 10. Here, ρ_w is the mass density of water.

From tables 3(a) and 3(b), it is revealed that the effect of vegetation tilt is to slightly increase the kinetic energy of the average flow velocity field k due to an increase in the velocity difference between the inside and the outside of the canopy, but it can have a slight influence in reducing the TKE k′ all over. On the whole, it has little effect on the total kinetic energy of the flow through vegetation. In the process of flow through vegetation, if the vegetation is tilted, the loss of pressure potential energy of the flow weakens. In general, in the RVMRE20(I) case, the energy loss of the flow through the vegetation zone is reduced by about one-third as compared with the RVMRE20(V) case. On the other hand, from tables 3(b) and 3(c), it can be seen that, due to the swaying of vegetation, the difference in velocity between the inside and the outside of the vegetation canopy increases, and the disturbance to the flow is aggravated. Therefore, when the flow passes through the flexible vegetation zone, its kinetic energy of the average velocity field k and the TKE k′ are both significantly larger than those in the rigid vegetation. In addition, the pressure potential energy loss $\Delta {\bar{E}_{tot}}$ of the flow passing through the flexible vegetation is more than 4.5 times that through the inclined rigid vegetation. This is also caused by an increase in the KH instability around the vegetation canopy. In general, due to the swaying of the vegetation, the energy loss of the flow through the flexible vegetation zone is about 4.46 times greater than that through the rigid vegetation zone. The above results show that the effect of vegetation tilt is to decrease its obstruction to the flow, while the effect of its swaying is to increase its drag force on the flow, which is consistent with the foregoing inference.

6. Conclusions

In this study, a 3-D numerical model of the flow–vegetation interaction based on LES and IBM has been established. By modelling plants into a structure of pellet–rope series, highly flexible plants subject to large deformation are simulated for the first time. In addition, the pellet–rope series structure can effectively simulate the stress characteristics of clustered leaves in flow, which is closer to the shape and mechanical characteristics of some specific plants than the elastic rod and flap structures commonly used in the previous studies. Furthermore, as a 3-D model does not limit the DOF of plant movement as compared to the 2-D model, this study can take comprehensive consideration of the transverse transport of momentum and energy loss caused by the movement of vegetation in the spanwise direction. The flow field distribution and the vortex structure analyses at different spanwise positions of the vegetation canopy can be more comprehensive. A physical model study was carried out in an experimental flume to validate the simulation results of the numerical model. The results show that the numerical model is effective in simulating the velocity field and the vegetation movement caused by the flow–vegetation interaction.

According to the results of the numerical simulations and the physical model study, it is found that the vegetation has a complex effect on the flow velocity field. During the flow through the vegetation zone, the streamwise velocity distribution varies across the flow depth. This is mainly reflected in the gradual increase in velocity difference between the inside and the outside of the vegetation canopy and the gradual increase in peak velocity above the vegetation canopy. It is also revealed that simplifying the effects of the vegetation on the flow and increasing the bed resistance to the flow lead to an erroneous flow velocity distribution. It is necessary to introduce a plant model with mechanical properties consistent with the actual conditions to study the flow–vegetation interaction. The main conclusions of this study are follows.

The flow Reynolds number directly affects the vegetation movement state. As the Reynolds number increases, both the force balance angle of the plants and the swaying amplitude in the streamwise direction increase. The force balance angle of the plants in the downstream zone is smaller than that in the upstream zone, but the swaying amplitude remains unchanged. In the spanwise direction, an increase in Reynolds number leads to an increase in swaying amplitude of the plants. From the perspective of the variational trend, the conclusion of this study is consistent with that of the previous studies on the analytical solution of the tilt angle of flexible vegetation (Luhar & Nepf 2011). However, in this study, we can obtain the velocity distribution at different places inside the vegetation more accurately, thus improving the calculation effect of the models in previous studies on forces on the individual plants.

The influence of the flexible vegetation on the flow can be divided into two parts: vegetation tilt and swaying around the force equilibrium point. The vegetation tilt does not significantly change the velocity difference between the inside and the outside of the vegetation canopy. However, due to the decrease in the effective canopy height, the thickness of the mixing layer decreases, and therefore the flow velocity outside the canopy also decreases slightly. The swaying of the vegetation increases the mixing layer thickness and the velocity difference between the inside and the outside of the vegetation canopy.

Vegetation tilt has little effect on the vortex structure and the vortex intensity in the flow. Owing to the regular distribution of vegetation, within rigid vegetation, the vorticity mainly presents a ribbon distribution, and the vortex structure formed is mainly HP vortices at a stem diameter scale. The swaying of the vegetation has a more severe impact on the vortex structure of the flow. Owing to the significant increase of velocity difference between the inside and the outside of the vegetation canopy, the large-scale KH–HP vortex structure is formed near the canopy height, which intensifies the energy and mass transport between vegetation and inside. Unlike the previous flexible vegetation models, this study can accurately simulate the spanwise swaying of vegetation; thus, the simulation of the KH–HP vortex structure scale is more reliable.

The vegetation tilt does not significantly change the TKE in the vegetation zone, but does decrease its distribution range, indicating that the effect of vegetation tilt is to reduce its disturbance to the flow. There is a small range of TKE diffusion within the inclined rigid vegetation, indicating that the effect of vegetation tilt is to slightly intensify the energy transport between the inside and the outside of the vegetation canopy. The swaying of vegetation increases the magnitude and the distribution range of the TKE and accelerates the TKE diffusion from the mixing layer to the interior of the vegetation zone. This indicates that the disturbance of vegetation to the flow can be significantly intensified by the swaying of the vegetation.

Based on the analysis of the flow velocity distribution, vortex structure, TKE distribution, and its production and dissipation rates for both flexible vegetation and rigid vegetation, it is recognised that the effect of vegetation tilt is to weaken its obstruction effect on the flow, while the effect of its swaying is opposite in nature. The energy loss of the flow through the vegetation zone also supports this inference. It can be predicted that the effects of both vegetation tilt and swaying on the flow patterns are affected by factors such as Reynolds number, distribution density and relative density of vegetation. Owing to the limited number of experiments carried out in this study, it is still not possible to obtain a general law of whether the obstruction effect of flexible vegetation on the flow increases or decreases under the combined action of vegetation tilt and swaying, which needs further research.