Numerical study of submerged bending vegetation under unidirectional ﬂow

Submerged vegetation commonly grows and plays a vital role in aquatic ecosystems, but it is also regarded as a barrier to the passing ﬂow. Numerical simulations of ﬂow through and over submerged vegetation were carried out to investigate the effect of vegetation density on ﬂow ﬁeld. Numerical simulations were computationally set up to replicate ﬂume experiments, in which vegetation was mimicked with ﬂexible plastic strips. The ﬂuid e structure interaction between ﬂow and ﬂexible vegetation was solved by coupling the two modules of the COMSOL packages. Two cases with different vegetation densities were simulated, and the results were successfully validated against the experimental data. The contours of the simulated time-averaged streamwise velocity and Reynolds stress were extracted to highlight the differences in mean and turbulent ﬂow statistics. The turbulence intensity was found to be more sensitive to vegetation density than the time-averaged velocity. The developing length increased with the spacing between plants. The snapshots of the bending vegetation under instantaneous velocity and vorticity revealed that ﬂexible vegetation responded to the effects of eddies in the shear layer by swaying periodically. The ﬁrst two rows of vegetation suffered stronger approaching ﬂow and were prone to more streamlined postures. In addition, the origin of tip vortices was investigated via the distribution of vorticity. The results reveal the variation of ﬂow properties with bending submerged vegetation and provide useful reference for optimization of restoration projects. © 2023 Hohai University. Production and hosting by Elsevier B.V. This is an open access article under the CC BY license (http:// creativecommons.org


Introduction
Submerged vegetation is a diverse group of rooted plants with their vertical extent less than the water depth.It commonly grows in fluvial and lacustrine areas, providing habitats for aquatic species, but it is also regarded as a barrier to the passing flow and removed from wetlands (L opez and García, 2001).
The resistance of submerged vegetation patches (VPs) apparently reduces the conveyance of waterways.However, the retarded flow in the vegetated zone remarkably inhibits the near-bed suspicion of sediment and significantly improves water clarity (Kosten et al., 2009;Wang et al., 2015a).The submerged vegetation also plays a vital role in aquatic ecosystems by uptake of nutrients and production of oxygen.Thus, it has been widely utilized for ecological restoration of wetlands (Vance, 2001;Taddeo and Dronova, 2018).As all the processes above are associated with ambient flow properties, a better understanding of vegetated hydrodynamics is required for optimization of restoration projects.
Experimental studies on flow through and above submerged vegetation have generally focused on KelvineHelmholtz (KeH) vortices at the top of the canopy (Ghisalberti andNepf, 2004, 2005;Siniscalchi et al., 2012;Huai et al., 2019).These vortices are generated due to the vertical discontinuity of vegetated drag, forming a coherent structure analogous to a mixing layer.The mixing layer grows downstream until shear production is offset by canopy dissipation and dominates vertical exchange processes between the vegetated zone and overlying zone.For a fully-developed flow field with submerged vegetation, hydrodynamic properties are complex in time and space due to the continuous alternation of sweep and ejection at the interface (Nepf and Ghisalberti, 2008).Comprehensive experimental investigations upon the shear layer and instantaneous flow field are difficult and costly.
Computational fluid dynamics (CFD) tools provide an efficient alternative by solving NaviereStokes equations with a numerical framework.There are two main approaches to process the vegetated turbulence.One is the Reynoldsaveraged-Navier Stokes (RANS), which includes different turbulence models to calculate the time-averaged Reynolds stress.By adding extra vegetation-induced terms to the momentum and turbulence transport equations (Fischer-Antze et al., 2001;Neary, 2003;Nicholas and McLelland, 2004), RANS models could successfully simulate the time-averaged flow field on a relatively coarse grid.However, their precision is significantly dependent on the drag coefficient of plants, an empirical parameter.As instantaneous information is lost in time averaging, the turbulent parameters of RANS simulations are sometimes unreasonable in comparison with experimental measurements (Defina and Bixio, 2005).To overcome the difficulties in modeling time-dependent turbulence, the large-eddy simulation (LES) model has been introduced to reproduce the production, propagation, and dissipation of KeH vortices at the top of the submerged vegetation (Stoesser et al., 2010;Yan et al., 2017;Gong et al., 2019;Zeng et al., 2022).With the LES model, instantaneous three-dimensional (3D) turbulent flow fields are almost directly resolved while small-scale motions are modeled.The vegetation area in the computational domain could be shaped with body-fitted grids.Thus, the interaction between flow and vegetation is explicitly computed instead of the use of a drag coefficient.Thus, LES results could elucidate large-scale coherent structures within the shear layer and agree with measured first-and second-order statistics.
Aquatic plants in most experiments and simulations have been idealized as rigid cylinders for simplicity (Tanino and Nepf, 2008;Stoesser et al., 2010;Wang et al., 2015b;Anjum and Tanaka, 2019).However, natural submerged plants exhibit a wide range of rigidity, and their reconfiguration is determined by fluid mechanics and biomechanics (Luhar and Nepf, 2011).With an increase in flow rate and flexibility, natural vegetation normally experiences the following three patterns: (1) gently swaying, (2) strong and coherent swaying, and (3) prone (Nepf and Vivoni, 2000).The Cauthy number (Ca; the ratio of drag force to rigidity force) of the first two states is small (Ca z O(1)), corresponding to bending plants categorized by Nikora (2010), while Ca of the last state is large (Ca is significantly greater than 1), corresponding to tensile plants.Recent studies have revealed that this classification significantly influences flowevegetation interactions (Luhar and Nepf, 2011;Okamoto et al., 2016;Pu et al., 2019).For example, the plant-flapping-induced turbulence is more prevalent with tensile vegetation.In order to simulate flow with flexible vegetation, reconfiguration should be considered in numerical models.Kutija and Hong (1996) introduced the theory of elastic cantilever beam to calculate the deflected height of bending vegetation.This approach was further developed and refined by following studies, but most of them only dealt with a static and steady status at each time step.Marjoribanks et al. (2014) developed a dynamic immersed boundary method with the LES model for bending vegetation using the EulereBernoulli beam equation and for tensile vegetation using the N-pendula model.Chen and Zou (2019) used an immersed boundary method to insert vegetation with large deflections into a RANSeVolume of Fluid (VOF) NaviereStokes solver for waveevegetation interactions.The immersed boundary method is robust for fluidestructure interaction (FSI), but the control volume of Lagrangian grids is difficult to quantify for a swaying plant, which might reduce the accuracy of simulations.In these coupled models, vegetation motion is generally computed with the finite element method (FEM), whereas fluid motion is resolved with the finite volume method or the finite difference method.The momentum term should be transformed between the rigidity force in the vegetation and the drag force in the fluid.Thus, one time-step lag was introduced between the two solvers.
This study aimed to investigate the variation of flow properties with bending submerged vegetation with increasing Ca.The LES model was used to simulate the unsteady and turbulent flow behavior at the top of the canopy.A global FEM solver for both the fluid motion and vegetation motion was used to achieve synchronous coupling.Simulations were conducted with an analogous laboratory experiment, in which plastic strips were used to mimic natural submerged vegetation, and observed data were used to validate the FSI model.

Numerical framework
The COMSOL Multiphysics package was used to resolve interactions between flow and bending vegetation, and it has been well verified and successfully applied in various cases (Baghel et al., 2020;Houda et al., 2017;Villalobos-Lara et al., 2020).The 3D flow field was simulated with the CFD module that uses LES to resolve turbulence.The velocity and pressure terms in NaviereStokes equations are divided into resolved and unresolved scales: where u and p are the resolved velocity and pressure, respectively; u 0 and p 0 are the unresolved velocity and pressure, respectively; r is the fluid density; m is the coefficient of viscosity of expansion; and t is the time.The momentum equation (Eq.( 2)) is then projected onto the finite-element sub-spaces of resolved velocity and pressure scales.The test functions for the resolved velocity and pressure scales are denoted as v and q, respectively.The projection U could be expressed as where f is the applied traction force on the boundary vU of the spatial domain U.The unresolved scales were predicted with the residual based variational multiscale method (Hashiguchi, 2012).The deformation of flexible vegetation was calculated in the module of structural mechanics using Hooke 0 s law as the basic theory for linear elastic bodies: where s is the strain tensor; ε is the second-order tensor; and C is the constitutive tensor, which is the fourth-order tensor.In the experiment, the vegetation was mimicked using thin plastic strips with negligible spanwise displacement, which could be simplified as plants in the numerical model.The stress at a certain level of plants could be computed based on the membrane stress (s m ) and bending stress (s b ): where k is a parameter ranging from À1 to 1; D is the plane stress constitutive matrix; d is the grain diameter; j is the membrane strains; j i is the initial membrane strains; x is the secant coefficient of thermal expansion and can be temperature dependent, T is the temperature; T m is the midsurface temperature, through which the influence of thermal strains is included; T ref is the reference temperature, DT is the temperature difference between the top and bottom surfaces of the shell; b h is the coefficient of hygroscopic swelling; c m is the midsurface moisture concentration; c m0 is the moisture concentration; c m0ref is the strain-free reference concentration; N i is the initial membrane forces; c is the bending part of the strain tensor (actually the curvature); c i is the initial value of the curvature; and M i is the initial bending and twisting moments.Other parameters were determined with the Structural Mechanics Module Manual.FSI coupling was conducted on the boundaries of fluid mechanics formulated with an Eulerian description and solid mechanics formulated with a Lagrangian description (Malaeb et al., 2018).The NaviereStokes equations above provide a solution of the flow field u fluid .Then the force exerted on the solid boundary is the negative of the reaction force on the fluid: where n is the outward normal to the boundary, and I is the identity matrix.The transformation of the force from the fluid frame to the material frame was carried out according to Eq. ( 9): where F is the force in the material frame; and d v and d V are the mesh element scale factors for the fluid frame and material frame, respectively.The two frames were unified in the FEM solver, and the geometric change of the fluid domain was automatically computed with the arbitrary LagrangianeEulerian (ALE) method.The moving mesh (ALE) in COMSOL is a dynamic mesh adaptation method that handles the dynamics of the deforming geometry and moving boundaries with moving grids (Sun et al., 2011).

Laboratory experiments and numerical setup
Numerical simulations were set up to replicate the flume experiments of Zeng and Li (2014).In their experiments, detailed measurements of the flow field were performed with The spanwise spacings between neighboring columns in both cases were identical (3 cm), and the streamwise spacing in Case 2 was doubled to be 6 cm.ADV measurements were carried out along the centerline of the flume, and gauging points were placed between vegetation columns.For a better comparison, the two cases were conducted in the same inflow conditions.Discharge at the inlet and water depth at the tail gate were fixed to be 13.89L/s and 0.246 m, respectively.Fig. 1 shows the computational domain with vegetation in Case 2. The bottom corner at the leading edge of the vegetated section was set as the origin of the coordinate, in which x, y, and z denote the streamwise, spanwise, and vertical directions, respectively.The computational domain has a span of L ¼ 7.9 m in the streamwise direction, consisting of a vegetated section and two gravel-bed sections.The gravel bed was simplified as a flatbed, and no-slip boundary conditions were used on the bed and side walls.A low Reynolds number model was used to compute near-wall velocities.Frictionless rigid lid conditions were used on the water surface that was relatively stable based on the measured data.In both cases, the cross-sectional averaged velocity (U ) was 0.31 m/s.Tetrahedral and hexahedral meshes were used to discrete the computation domain, and the number of elements was approximately 7.7 Â 10 6 .Simulations were initially run for a flow-through time (t L ¼ L/U ) to develop the flow and then continued for a period of 10t L to collect time-averaged flow statistics.Simulations were performed with 64 cores (2 Â AMD EPYC 7452) on Paracloud 0 s supercomputer, and each simulation consumed approximately 4 200 CPU hours.

Results and discussion
The numerical results of the two cases were validated with the experimental data for the time-averaged flow properties.Afterwards, turbulent statistics were used to characterize the mixing layer at the top of the canopy in each case.Finally, the simulated instantaneous flow field was used to identify the motion mode of the flexible vegetation.

Time-averaged flow characteristics
In numerical simulations, flow could pass through the vegetation array along the spanwise gaps between rows, rather than just divert to the overlying layer, because the vegetation was explicitly resolved with fine grids.Two time-averaged streamwise velocity (u) profiles along the vegetated section were extracted as shown in Figs. 2 and 3.In both cases, the approaching flow decelerated due to the submerged plants, and the flexible plants were pushed in bending postures.The deflections of leading rows were greater than those of downstream rows due to the stronger bending force exerted by the incoming flow.Above the top of the canopy, flow velocity increased rapidly before reaching the water surface.In contrast, under the canopy, flow velocity was nearly uniform regardless of the water depth.In Case 1, the existence of the submerged vegetation significantly retarded the ambient flow (Fig. 2  damping effect of vegetation was strong enough to overcome half of the spanwise distance between vegetation arrays, which was five times the width of plants.Similar velocity distributions could also be found in the experiments and simulations of flow with submerged vegetation (Liu et al., 2008;Stoesser et al., 2009).The developing process of incoming flow also occurred in the forepart of the vegetated zone, i.e., a highspeed wedge penetrated the low-speed zone from the leading edge.Fig. 3 clearly demonstrates this process.The passage of the upstream flow was narrowed by the vegetation arrays.Fig. 3 shows the time-averaged streamwise velocity profiles in Case 2. The most notable distinction between the two cases was the developing length in the transitional section above the vegetation.As the vegetation density decreased, the velocity gradient at the upper layer significantly declined, and the flow rate in the vegetated zone slightly increased.As a result, the mean deflection of the vegetation elements in Case 2 was greater than that in Case 1.However, other flow features did not remarkably changed, compared to the two cases.
Fig. 4 shows the velocity distribution in the horizontal cross-section at the half height of the vegetation (z ¼ 0.075 m).The upstream flow was pushed into narrow gaps between the vegetation arrays.The flow rate rapidly increased in the entrance zones and gradually decreased along the sub-channels.The loss of the streamwise velocity was not only caused by the resistance of vegetation but also induced by the vertical transport of momentum and the dissipation of turbulent energy.The streamwise variation of velocity was minimized after a certain developing length that was longer in Case 2 due to the lower vegetation density.The difference between the velocity distributions of the gaps could be neglected.Thus, the vertical profile of the streamwise velocity in the middle gap of the fully-developed zone was extracted to validate the numerical model.Fig. 5 compares the simulated vertical profiles of the streamwise velocity at different locations in the fully-developed zone with the measured data in the two cases.Gauging points were set along the centerline of the flume, and the velocity distributions at different locations were not significantly different, because the time-averaged flow field was fully developed.The velocity in the vegetated zone showed a slow growth in the vertical direction, whereas that in the upper layer showed a rapid growth.Thus, an inflection point appeared at the top of the bending vegetation.As shown in Fig. 5, the measured data were evenly distributed around the simulation results, showing that the numerical model could accurately predict the time-averaged streamwise velocity.

Turbulent flow characteristics
The turbulent flow through the submerged vegetation was characterized with the structure of coherent vortices at the stem top, which could be classified as shear layers.Experimental studies have revealed that the growth of shear layers ceases once the production of the shear-layer-scale turbulent kinetic energy is balanced by turbulent dissipation within the canopy.The thickness of the shear layer is roughly constant beyond the balance point.The shear-layer-scale turbulence could be assessed by the Reynolds stress ( À u 0 w 0 , with u 0 and w 0 denoting the fluctuating streamwise and vertical velocities, respectively) or turbulent kinetic energy.Figs. 6 and 7 show the Reynolds stress normalized by the shear stress at the stem top (u*) of the two selected lines in cases 1 and 2, respectively.Peak Àu 0 w 0 values in both cases appeared at the interface between the submerged vegetation and overlying water.The turbulence intensity gradually increased from the leading edge of the vegetated section, demonstrating the development of the shear layer.As the shear layer is dominated by the streamwise and vertical instabilities, the development length can be quantified via the analysis of the streamwise variation, and the vortex size can be evaluated with the vertical distribution of À u 0 w 0 .
Unlike the mean flow field, the turbulence field was significantly impacted by the vegetation density.In Case 1 with dense vegetation, the range and intensity of the shear layer at the stem top were significantly greater than those in Case 2. The high values of the normalized Reynolds stress ( À u 0 w 0 =u *2 ) (greater than 0.8) were shown as a continuous band in Case 1 but shown as discrete spots in Case 2. The discrepancy was mainly caused by the difference of the spacings between the neighboring rows in the two cases.In Case 1, the spacing was sufficiently narrow for vortices to overcome dissipation before reaching the next tip.The canopy dissipation (W) defined by Ghisalberti and Nepf (2004) is where a is the frontal area of the vegetation per unit volume, C D is the drag coefficient of the canopy, and v 0 is the fluctuating spanwise velocity.The frontal area index (a) is directly proportional to vegetation density, while the drag coefficient (C D ) is typically positively correlated with vegetation density at the same Reynolds number.The superposition of the two parameters led to stronger canopy dissipation in Case 1.In addition, the effective frontal area decreased with the increase of deflection, further reducing canopy dissipation in Case 2. As a result, the growth of the Reynolds stress was terminated before the trailing edge of the vegetated section, indicating a fully-developed state of the shear layer because of sufficient dissipation.However, the relatively weak canopy dissipation could not balance the turbulent energy in the range of the vegetated section, thereby causing the value of À u 0 w 0 in Case 2 to continuously increased, even outside the tail edge.Fig. 8 compares the vertical distributions of the absolute Reynolds stress in the two cases with the experimental data.
The Reynolds stress was extracted at the end of the vegetated section as the shear layer was fully developed at the trailing edge in Case 2. The consistency of the measured Reynolds stress at different points was apparently lower than that of the measured data of the time-averaged velocity.Nevertheless, the simulation results of the two cases agreed well with the experimental measurements.The fluctuation of the two curves was caused by the vibration of vegetation and/or the insufficient averaging time span.The simulated peak Reynolds stress decreased from 0.002 47 m 2 /s 2 to 0.001 13 m 2 /s 2 because the vegetation density was reduced by 50% in Case 2. It should be noted that the peak Reynolds stress in the two cases was not exactly located at the mean height of bending plants (the dashed line in Fig. 8).For flexible vegetation, the vertical intrusion of shear-layer vortices into the canopy is promoted by the vegetation vibration, but the Reynolds stress is reduced in this process, which is the so-called monami (Ackerman and Okubo, 1993).In Case 1 with dense vegetation, vegetation vibration and subsequent intrusion were inherited due to the narrow spacing.As a result, vortex intrusion was extremely   low in Case 1 but significantly noticeable in Case 2. Furthermore, the reduction of the peak Reynolds stress was partly induced by the deeper vortex excursion.Although under the same inflow conditions, vegetation showed different properties due to the decrease in vegetation density.It could be inferred that the effects of monami would be further enhanced with low vegetation density.

Fluidestructure interaction
The interaction between flow and rigid vegetation mainly reflects the influence of vegetation upon flow.It can be simplified as a one-way coupling problem, i.e., the impact of flow on rigid vegetation is not considered.However, flow through flexible vegetation is a two-way FSI coupling problem, in which vegetation deforms within the passing flow and flow properties are altered due to the deformation.In the simulation scheme of this study, the CFD solver was strongly coupled with the structural solver, which allowed for the computation of the interaction between flow and swaying vegetation.
Previous results revealed that the reconfiguration of vegetation elements in cases with sparser vegetation was relatively greater.Thus, Case 2 was selected to study the bending vegetation under instantaneous streamwise velocities at three instances (Fig. 9).Unlike the time-averaged results, the postures of the plants varied, which was dependent on their locations in the flow field.Most of them bended into streamlined shapes, while others bounced back to the upstream direction.The plants at the leading edge deflected significantly and vibrated violently due to the strong incoming flow, and the swing amplitude of the side plants slightly decreased because of the retarded near-wall flow.
The instantaneous flow field was much more chaotic than the time-averaged flow field, especially in the zone above the canopy.The flow rate of the overlying water normally accelerated when the underlying vegetation was pushed into streamlined postures as shown in the velocity distributions in Figs.9(a) and (b).Fig. 9(c) shows that the second row was recurved towards the first row, and some plants moved backward against the flow, thereby slightly decelerating the overlying flow.Nevertheless, the velocity profiles at the three instances were highly turbulent.Fig. 10 demonstrates that turbulent structures were mainly generated at the tips of the plants.Boundary layers formed along the upstream surface of the plants and subsequently separated at the tips, shedding continuous structures above the canopy.It should be noted that the negative spanwise vorticity was absent at the top of plants.This indicated that the vibration was not significant, compared to the clockwise and counterclockwise vortices generated by the advance and return motions of the tips.

Conclusions
In this study, flow through submerged flexible vegetation was investigated using a fluidestructure interaction technique in COMSOL packages.The numerical model was verified through a typical flume experiment with submerged vegetation in two densities under the same inflow conditions.The comparison of the time-averaged velocity and Reynolds stress in fully-developed regions showed that the simulation results agreed with the measured data in both cases.The numerical results revealed that submerged vegetation significantly altered the flow field and vegetation motion.The distribution of the streamwise velocity with an inflection point at the top of the deflected vegetation was found whereas the location of the peak Reynolds stress was slightly lower than the vegetation top due to the vortex-driven oscillation.Although the influence of vegetation density on the velocity distribution was not significant, the turbulence level was significantly reduced in the case with sparser vegetation due to the weaker equivalent stiffness and enhanced vortex-driven oscillation vibration.Meanwhile, the shear layer continuously increased along the top of the canopy until a balance point was reached.In the case with denser vegetation, the developing length was shorter, and the turbulence intensity was stronger.The tip vortices induced by vegetation motion was insignificant.This indicated that the flow properties in this study were similar to those of rigid submerged vegetation.
an acoustic Doppler velocimeter (ADV).The length and width of the flume were 12.50 m and 0.31 m, respectively.The vegetated section, located in the middle part of the flume, had a width of 0.30 m and a length of 2.40 m.Two gravel-bed sections were placed at the immediate upstream and downstream sides of the vegetated area.Flexible vegetation was represented by plastic strips with rectangular cross-sections, and they were fixed on a perforated baseboard in a regular pattern.The plastic strips had a height of h v ¼ 0.15 m, a width of b v ¼ 0.005 m, a thickness of t v ¼ 0.001 5 m, and a flexural rigidity of E I ¼ 4.75 Â 10 À3 N/m 2 .Two different vegetation densities were considered: 1 111 elements per square meter for Case 1 and 556 elements per square meter for Case 2.
Fig. 1.Vegetated section of computational domain in Case 2.

Fig. 2 .
Fig. 2. Profiles of simulated time-averaged streamwise velocity along selected lines through plants and between plants in Case 1 with dashed box denoting vegetated zone.

Fig. 4 .
Fig. 4. Simulated time-averaged streamwise velocity in horizontal plane at half height of erected vegetation in Case 1.

Fig. 3 .
Fig. 3. Profiles of simulated time-averaged streamwise velocity along selected lines through plants and between plants in Case 2 with dashed box denoting vegetated zone.

Fig. 5 .
Fig. 5. Vertical profiles of time-averaged streamwise velocity in fully-developed zone in Case 1 and Case 2.

Fig. 6 .
Fig. 6.Simulated profiles of normalized Reynolds stress along selected lines through plants and between plants in Case 1.

Fig. 7 .
Fig. 7. Simulated profiles of normalized Reynolds stress along selected lines through plants and between plants in Case 2.

Fig. 8 .
Fig. 8. Vertical profiles of Reynolds stress at the end of vegetated zone at different locations in Case 1 and Case 2 with dashed line representing average height of bending plants.