Steady Flow Over a Finite Patch of Submerged Flexible Vegetation

An immersed boundary‐finite element with soft‐body dynamics has been implemented to study steady flow over a finite patch of submerged flexible aquatic vegetation. The flow structure interaction model can resolve the flow interactions with flexible vegetation, and hence the reconfiguration of vegetation blades to ambient flow. Flow dynamics strongly depend on two dimensionless parameters, namely vegetation density and Cauchy number (defined as the ratio of the fluid drag force to the elastic force). Five different flow patterns have been identified based on vegetation density and Cauchy number, including the limited reach, swaying, “monami” A, “monami” B with slow moving interfacial wave, and prone. The “monami” B pattern occurred at high vegetation density and is different from “monami” A, in which the passage of Kelvin‐Helmholtz billows strongly affects the vegetation interface. With soft‐body dynamics, blade‐to‐blade interactions can also be resolved. At high vegetation density, the hydrodynamic interactions play an important role in blade‐to‐blade interactions, where adjacent vegetation blades interact via the interstitial fluid pressure. At low vegetation density, direct contacts among vegetation blades play important roles in preventing unphysical penetration of vegetation blades.


Introduction
Submerged aquatic vegetation (SAV) plays important roles in riverine and estuarine processes.By attenuating flow and modulating turbulence, SAV promotes sediment stability and retention of particles (Nepf & Ghisalberti, 2008;Triska et al., 1989;Vargas-Luna et al., 2015).SAV also provides food and shelter for ecologically and economically important species (Heck & Thoman, 1984;Schein et al., 2012) and improves water quality by removing nutrients, providing oxygen, and reducing turbidity (Dennison et al., 1993;Duarte, 1995;Orth et al., 2010).SAV, therefore, forms a critical link between the physical habitat and the biological community (Bornette et al., 1998;Moss, 1990;Rios & Bailey, 2006), and understanding the physical-biological couplings between SAV structure/function and flow regime is of critical importance for natural resources management (Nepf, 2012;Nikora, 2010).Recently, Reaver et al. (2019) identified that the continuous, periodic SAV canopy motion ("monami") inhibits algae growth, motivating interest in the role of blade-to-blade interactions on ecosystem structure and function.
For flexible vegetation, as flow velocity increases, vegetation blades bend, and the drag by SAV reduces.Deflected vegetation height strongly affects hydrodynamics (Jarvela, 2002;Luhar & Nepf, 2011, 2013), since the projected area scales with effective vegetation height.Kouwen and Unny (1973) proposed an empirical function of deflected plant height that depends on the Cauchy number (defined as the ratio of the elastic restoration to the hydrodynamic force).Stephan and Gutknecht (2002) investigated three species of SAV and found the deflected plant height to be an appropriate parameter to describe the drag force.However, in dense vegetation canopies, flow sheltering likely plays a role in the deflection height, since stems directly behind other stems experience lower velocity and drag force, and thus exhibit different bending properties compared to a single blade.Furthermore, when vegetation blades come into contact with each other, blade-to-blade interactions may also alter the deflected height.As such, attributing observed deflection height solely to bending stiffness may not be appropriate.The influence of interactions among multiple blades, including their hydrodynamic interaction and direct contacts, have yet to be investigated; however, they are likely important for describing the drag by the whole canopy.transport of momentum, nutrients, and sediment (Ghisalberti & Nepf, 2006).In addition, theoretical studies (Cáceres-Euse et al., 2022;Singh et al., 2016) and 3D numerical models (Fischer-Antze et al., 2001;Li & Xie, 2011) have been implemented to understand the properties of Kelvin-Helmholtz billows in submerged vegetation canopy.LES models can resolve turbulent scales of interest and have been applied to study flow over submerged rigid vegetation (Chakrabarti et al., 2016;Okamoto & Nezu, 2010;Stoesser et al., 2009), and Kelvin-Helmholtz instabilities at the top of the vegetation canopy (Cui & Neary, 2008).With efficient fluid-structure-interaction models, LES models have also been applied to study interactions between flexible vegetation and time-varying flows (Ikeda et al., 2001;Tschisgale et al., 2021).The penetration depth of turbulent motions depends on the flexibility of the vegetation blades.In addition, blade-to-blade interactions may affect the motions of small-scale turbulent eddies by altering the interstitial space among vegetation blades.Despite these advances in laboratory measurement and numerical simulation techniques, the effects of interactions between vegetation blades have not yet been thoroughly investigated or modeled.
Immersed boundary method (Kadoch et al., 2012;Peskin, 2002) (IBM) has been widely used to study flow structure interactions.O'Connor and Revell (2019) implemented a lattice Boltzmann immersed boundary finite element model to study steady flow over an array of wall-mounted flexible flaps at low Reynolds number.The model was first validated with laboratory experiment data of oscillatory flow over a row of flexible flaps (Favier et al., 2017;Revell et al., 2017) and then applied to investigate the coherent waving interactions between the vegetation and flow.Tschisgale and Fröhlich (2020) developed an efficient and robust semi-implicit immersed boundary method for the fluid-structure interaction of slender flexible structures in viscous fluid.They validated the model with extensive 2D and 3D configurations of slender rods and then applied the model to study flow over an abstracted aquatic flexible vegetation canopy (Tschisgale et al., 2021).The high-fidelity large-eddy simulation model revealed the well-known coherent turbulent structures, including velocity streaks, Kelvin-Helmholtz vortices in the mixing layers as well as hairpin vortices in the outer flow region, and hence the turbulent generation mechanisms by vegetation.He et al. (2022) conducted large-eddy simulations to study turbulent aquatic canopy flows with flexible stems and carried out analyses of turbulent energy budget.The authors showed that the wake production can be as large as one-half of the shear production near the canopy top.Both LES studies (He et al., 2022;Tschisgale et al., 2021) focused on an infinite vegetation canopy in a turbulent channel, which neglected the spatial heterogeneity at the canopy edge as in a finite vegetation canopy.Wang et al. (2022) and Ni et al. (2023) conducted high-fidelity simulations of flow over an array of flexible vegetation with finite patch sizes.Wang et al. (2022) used plants consisting of five spheres constrained by a fixed distance between adjacent spheres, finding that both the flow velocity and vegetation spacing affect the wavelength and amplitude of the coherent waving motion.Ni et al. (2023) used IBM to study "monani" modes and scales of an array of flexible vegetation in a laminar boundary layer and identified four different modes based on non-dimensional gap distance (defined as the ratio of the gap distance to the length of the vegetation blade) and Cauchy number.
Previous studies of flexible vegetation canopy did not include direct interactions among vegetation blades.To include the direct blade-to-blade contacts, a collision detection and response algorithm needs to be implemented.Chen et al. (2022) studied the flapping dynamics of vertically clamped three-dimensional flexible flags in a Poiseuille flow, in which a repulsive force was introduced to prevent the collision of flags.Zhang et al. (2020) adopted the same approach to study fluid-structure interactions of single and dual wall-mounted 2D flexible filaments in a laminar boundary layer and then adopted the framework to study flow over a large array of wall-mounted flexible plates and investigated different types of flow instabilities of different dynamic modes (Zhang et al., 2022).In these studies, the artificial short-range repulsive force modeled by a delta function acted as the fluid lubrication force.The solid objects interacted repulsively via the intervening fluid and did not directly collide when they were in close proximity.In natural ecosystems, the direct contacts among vegetation blades are of critical importance in the removal of algae attached on the vegetation blades (Reaver et al., 2019), and the modeling of the contact mechanics needs to be included.To model direct physical contacts among vegetation blades, a soft body dynamic approach with contact force (Rong et al., 2019;Shabana, 1997;Skrinjar et al., 2018) can be applied.Intriguingly, the computer graphics community has used the soft body dynamics approach with collision detection and response in physics engines (such as Bullet physics engine) to model dynamic systems such as grass (Chen & Johan, 2010;Lo et al., 2016;Selino & Jones, 2013), fur, and hair (Iben et al., 2013;Lengyel et al., 2001).This paper is organized as follows.First, we presented the coupled model framework of immersed boundary method and finite element method with soft-body dynamics in Section 2. We extensively validated the ANCF finite element method and the flow-structure interaction (FSI) method with Turek FSI benchmark case in Section 3. The validated model was then applied to study steady flow over a finite patch of submerged flexible vegetation with a detailed model setup in Section 4. By varying the vegetation and flow parameters, the dependence of flow dynamics on the important dimensionless parameters, namely the vegetation density and Cauchy number, were elucidated in Section 5.The importance of the hydrodynamic interactions was discussed in Section 6.Finally, conclusions were summarized in Section 7.

Methodology
In this study, we developed a novel multiphysics model framework that couples a soft-body dynamics ANCF finite element model with a fluid solver.Compared to previous studies of flow over submerged flexible vegetation (He et al., 2022;Wang et al., 2022), one of the advantages of the newly developed model framework is its ability to capture both flow-vegetation and blade-to-blade interactions, which allows the model to be applied to a wider range of vegetation density and Cauchy number compared to previous studies.

ANCF Formulation of Thin Vegetation Blades
In the absolute nodal coordinate formulation (ANCF), the global position vector (r) of an arbitrary point on a three-dimensional element can be written as where S is the element shape function matrix.x, y, and z are the local coordinates of the element, and e is the vector of the nodal coordinates, which is time dependent.
For a three-dimensional two-noded Bernoulli-Euler beam element (I, J), 12 nodal coordinates are used for each node (Figure 1).The coordinates of node I is = . (2) The local coordinates of the element can therefore be written as (3) The element shape function S is given as in which I is a 3 by three identity matrix.The shape functions are given as.
The length of the beam element in the initial undeformed configuration is represented by l, and the non-dimensional quantity ξ is defined as ξ = x/l.
The global velocity of an arbitrary point of the beam can then be written as (  ; ) = ̇(  ; ) = (  ) ̇ (9) in which the dot denotes the derivative with respect to time.The mass matrix of the beam element can be derived from the kinetic energy The mass matrix is given as For the three-dimensional beam element, with a given shape function matrix S, the mass matrix can be derived following the above equation, and remains the same when the beam element deforms.

Volume Penalization Immersed Boundary Method
In the volume penalization immersed boundary method, the solid objects, such as vegetation stems, are modeled as porous media with low permeability.The governing equations are given as follows: where ρ is the fluid density, ν is the kinematic viscosity of the fluid, and p is the dynamic pressure.x i , u i , and u s,i are the i-th component of the position vector, fluid velocity vector, and solid object velocity vector.i = 1, 2, 3 represents the streamwise, spanwise, and vertical directions, respectively.The penalty parameter ϵ can be interpreted as permeability.Since the flow velocity passing through the porous media increases with ϵ, a very small value is often implemented in the immersed boundary method.To find the optimal value for the penalty parameter ϵ, sensitivity analysis has been conducted following Kadoch et al. (2012) and Yu and Yu (2022).
The mask function α is defined as where Ω s is the solid phase, and Ω f is the fluid phase.The no-slip boundary conditions at the fluid-solid interface are no longer included.The solid object is represented by the mask function α.To avoid numerical instabilities, a smooth transition of the mask function α at the solid-liquid interface with values between 0 and 1 is required.
The thickness of the transition layer is defined by δ t .The thickness parameter δ t measures the spreading of the influences from the solid structure (such as forces or stresses) to the fluid and is determined based on the grid spacing.
In the present work, we use an analytical expression for the mask function, which can be easily resampled on the fluid grid.The 1D smoothed Heaviside function is implemented as The smoothed Heaviside function provides a good way to alleviate numerical instabilities at the sharp fluid-solid interface.The fluid-solid interface is located at x 0 , and the solid is at x ≥ x 0 .For 2D problems, the tensor product of the 1D Heaviside function is implemented at the fluid-solid interface as H(x, z) = H(x)H(z).The optimal value of the thickness parameter (δ t ) is determined by the grid spacing Δx.In this study, the value of δ t = 2Δx was implemented, which is consistent with the regularized delta function for the interpolation of the Eulerian fluid pressure to the Lagrangian points on the vegetation blades.Figure 2a shows the implementation of the mask function of two interacting blades.The beam centerlines are shown as blue dashed curves, and the mask function is shown in red based on the distance to the centerline.

Flow-Structure Interaction
The regularized delta function was implemented to smoothly interpolate the pressure from the Eulerian fluid field (x j ) to the Lagrangian points (x i ) on vegetation blades.For 1D problems, the pressure at point x i is given as The regularized delta function is given as In 1D problems, r = (x − x p )/Δx is the normalized distance of point x to the Lagrangian point x p , and Δx is the grid size.The tensor product of the 1D delta function is used for 2D and 3D problems.
The velocity of the solid object v i is projected from the Lagrangian points on the centerlines of the vegetation blades to the Eulerian fluid.We first identified the nearest points X p,j and X p,j+1 on the vegetation blade by the distance for the fluid cell at x i .The linear interpolation was then used to obtain the velocity.

Blade-To-Blade Interactions
Collision detection and response have been widely used in multibody dynamics with computer graphics and robotics applications to prevent unphysical interpenetrations.Collision detection algorithms are implemented to detect the intersection of two or more objects.In order to optimize the collision detection process, it is generally split into two phases, namely the broad phase and the narrow phase.The broad phase is responsible for finding pairs of objects that are potentially colliding and excluding pairs that are certainly not colliding.The narrow phase operates on the pairs of objects found by the broad phase and determines if the collisions actually occur.
Space partitioning coupled with bounding volumes that establish an upper bound for collision is commonly used in the broad phase, such as the axis-aligned bounding box (AABB) method and oriented bounding box (OBB) method.An AABB is simply a rectangular parallelepiped box whose faces are each perpendicular to one of the basis vectors.The AABB method is computationally cheap for simple geometries, such as boxes, cones, and cylinders, by using rectangular (cuboid in 3D) bounding boxes.Since the beam elements can be modeled as simple box geometries, AABB has been implemented in this study (Figure 2).
The narrow phase identifies the actual collisions on the pairs of bodies found by the broad phase that might be colliding.It is a refinement step, in which we compute the contact points that will be used to compute the collisional forces for each collision found in the narrow phase.In this study, we approximate the vegetation blade as chained-spheres, and hence the collision detection in the narrow phase was greatly simplified as collision detection of spherical particles.The distance between the centers of the two spheres was used for collision detection.
After a collision has been detected, the discrete element method (DEM) is implemented to compute the collisional forces.The main purpose of the collision force is to push colliding objects away from each other to prevent unphysical interpenetration.The cross-section area of the vegetation blade was uniform (Figure 2b), and hence, the thickness of the vegetation blade was used as the diameter of the chained-spheres.The particles in DEM were modeled as soft and identical spheres.Hookean contact model was implemented in this study to compute the collisional force, which is a linear function of the overlap distance δ p as in which x is the position vector at the center of the particle and D is the diameter of the particle.The subscript "i" and "j" represent the label of the particle.
For simplicity, only the normal force was implemented.The normal force has two terms, a contact force and a damping force.K n is the elastic coefficient for normal contact, and η n is the viscoelastic damping constant for normal contact.n ij is the unit vector along the line connecting the centers of the two spheres.v n is the normal component of the relative velocity.The timestep for DEM simulations is restricted by the elastic constant (or Young's Modulus of the material), which often leads to a small timestep.To accelerate the simulation, a value that is several order of magnitude smaller than the value based on material properties is often used.In this study, the elastic contact force is mainly used to prevent unphysical penetration of vegetation blades, and a reduced elastic coefficient will not affect the reconfiguration of vegetation blades.

Model Validation
The immersed boundary method has been validated previously (Yu & Yu, 2022), and we focused here on the validation of the structure solver and the fluid-structure interaction (FSI).

ANCF Structure Solver-Cantilever Beam
A cantilever beam is a rigid structural element that extends horizontally and is supported at only one end.In ANCF formulation, the dynamics of a one-dimensional beam element are governed by Newton's second law as follows: in which M is the mass matrix,  ë is the second order time derivative of vector e, Q e represents the restorative elastic force, Q p represents external point loads, and Q q represents external distributed load.
The solution to the above equation is a harmonic oscillator.To obtain the static solution, an artificial structural damping term is needed to dissipate the mechanical energy of a vibrating structure.The viscous damping term in the form of    =  ̇ is introduced, which is linearly proportional to the velocity of the beam  ( ̇) , with η as the damping coefficient and m as the mass of the beam element.
A concentrated point load is applied at the free end of the cantilever with a uniform cross-section area (Figure 3).The tip deflection (δ z ) is used to validate the model.For a cantilever beam with a small deflection, the bending curve of the beam is given as which gives the deflection at the free end E is the Young's modulus, and I is the moment of inertia.For large deflection, if a concentrated load F parallel to z axis is applied at the free end of a uniform cantilever beam, the whole bending curve of the beam can be solved using the following equation (Chen, 2010) The above equation is solved numerically to obtain the deflection at the free end.
The ANCF model has been validated against both small deflection and large deflection.Four elements were used to model the cantilever beam, and the ANCF model works reasonably well with as little as four elements for the one-dimensional cantilever beam with both small and large deflection.The error increased with the magnitude of the point load due to numerical errors since only four elements were used in the validation.

Flow-Structure Interaction
The widely used 2D benchmark Turek FSI (Turek & Hron, 2006) The computational grid was set to Δx = Δz = 2.5 mm, and 8 grids were used to resolve the thickness of the beam.The thickness parameter δ t was set to 0.5 mm.The penalty parameter ϵ = 2 × 10 −4 s.A total number of 100 elements were used to model the flexible beam, leading to the element size of 1 mm.The time step was set to 1 × 10 −4 s and the maximum CFL number was below 0.2.
Figure 4 shows the model validation.The vorticity field ω y at t = 10s shows the modulation of the vortex shedding by the flexible beam, in which the side walls limited the growth of the vortex.At steady state, due to the periodic vortex shedding from the circular cylinder, the flexible beam oscillated in the flow periodically as well.
The vertical tip position δ z (Figure 4b) was used to calibrate the model.The model predicted the amplitude of 32.02, compared to 34.38 mm in the benchmark (Turek & Hron, 2006).Our model performed reasonably well, with an error of 7%.

Model Setup
A finite patch of submerged flexible vegetation was placed in the middle of the computational domain (Figure 5).The total number of vegetation blades was set to 32, and the vegetation blades were uniformly distributed in The open-source package OpenFOAM (Weller et al., 1998) was used as the fluid solver.Inflow boundary condition was implemented at the inlet with a uniform velocity profile.Outflow boundary condition was implemented at the outlet, which is a mixed type of boundary condition.When the flow exits the computational domain, Neumann boundary condition of flow velocity is used.Flow velocity is set to zero when reverse flow back into the computational domain occurs.No-slip boundary condition was implemented at the bottom boundary, and free-slip boundary condition was used at the top boundary.A uniform grid was used with Δx = Δz = 0.5 mm, and the total number of grid points was between 2 and 2.5 million.Only 4 grid points were used to resolve the thickness of a single vegetation blade since it was modeled as a slender body, and we only focused on the flow over its length.The penalty parameter was set to ϵ = 5 × 10 −3 s, and the thickness parameter was set to δ t = 0.33 mm.A constant time step of Δt = 2 × 10 −4 s was used, and the maximum CFL number was around 0.4.The timestep is also sufficiently small for DEM calculations to avoid the unphysical penetration of vegetation blades, as the distance that one particle can travel in one timestep is much smaller than the diameter of the particle (thickness of the vegetation blade).
In total, 15 simulations were conducted by varying the vegetation density, Young's modulus, and inflow velocity.Table 1 summarizes the parameters for all simulations.The flow dynamics strongly depend on two key dimensionless parameters, namely vegetation density and Cauchy number.For 2D simulations, the vegetation density (ϕ, represented solid volume fraction) is defined as the ratio of the thickness of the vegetation blade (t) to the distance between two adjacent vegetation blades (s).The Cauchy number is defined as in which b is the width of the blade, l is the length of the blade, E is the Young's modulus, I = bt 3 /12 is the moment of inertia, t is the thickness of the blade, and ρ is the fluid density.Cauchy number is defined as the ratio of the Note.The last column represents the flow pattern in the five types, namely limited reach, sway, "monami" A, "monani" B, and prone.

Table 1 Summary of Parameters for All Simulations
hydrodynamic drag force on the blade to the elastic restoration force, and measures the rigidity of the plant.
Cauchy number can determine the magnitude of the bending motion for flexible vegetation.As Cauchy number increases, the reconfiguration of flexible vegetation to ambient flow becomes more pronounced.
By default, both OpenFOAM and the structure solver are three-dimensional solvers.To conduct 2D simulations, one grid cell in the spanwise y direction was implemented with empty boundary condition at the lateral boundaries in OpenFOAM.In the structure solver, we restricted the structure deformation to be 2D, in which twisting motion and displacement in spanwise y direction were disabled.
The Reynolds number (Re = Ul/ν with ν as the kinematic viscosity of water) for these simulations is sufficiently high, and a turbulent closure model is required.The standard two-equation RANS models such as k − ϵ or k − ω models do not include the slip production due to the relative motion between the vegetation blades and the ambient fluid, which could be important for turbulence generation.For simplicity, the Smagorinsky model was implemented as in LES.The eddy viscosity is given as where C s is the Smagorinsky constant and S ij is the strain rate tensor.The overline represents the resolved scale, and Δ is the filter-width length scale given as  Δ =  1∕3  with V cell representing the volume of the computational grid.
Natural submerged vegetation has different plant morphology, which also affects the hydrodynamics (Boothroyd et al., 2016;Xu & Nepf, 2020).For simplicity, in this study, the vegetation was simulated as a slender beam with a rectangular cross-section, which is representative of submerged vegetation with long thin leaves such as eelgrass (Zostera spp.) and tape grass (Vallisneria spp.).The thickness of the vegetation is on the order of 0.1-1 mm.The hydrodynamics of flow over submerged vegetation canopy are primarily controlled by three dimensionless parameters, namely Reynolds number, Cauchy number and Buoyancy number (Luhar & Nepf, 2011;Whittaker et al., 2015).In this study, the material density of the vegetation was set to be the same as the ambient fluid, and hence buoyancy was neglected.In the natural environment, the Cauchy number of vegetation canopy ranges from O(0.001) to O(1,000).In this study, the range of Ca was chosen to be from 4.22 to 480, focusing on more flexible vegetation.The modulus of elasticity of natural vegetation ranges from the order of MPa to GPa (Kouwen & Li, 1980).Reduced Young's moduli (compared to that of natural vegetation plants) were used to characterize a wider range of Ca.

Flow Patterns
We focused on flexible vegetation and did not include the rigid (or erect) pattern in this study, in which the vegetation blades do not change their tip position in time.Five distinct flow patterns have been identified, which depend on the vegetation density (ϕ) and Cauchy number (Ca).
At high vegetation density (ϕ = 0.133) and low Cauchy number (Ca = 4.22) in Case 7, the vegetation blades bent in flow direction initially at the leading edge of the canopy.The vegetation blades at the interior gently swayed with the flow due to interactions with the vortices formed at the top of the canopy (Figure 6a).The vortices propagated downstream and grew in size (Figure 6c) at the shear layer.At a later stage, a gap formed with trapped vortices (Figure 6b), and the vegetation canopy was divided into two regions.The vegetation blades at the leading edge formed a cluster and vibrated slightly as a whole.Vegetation blades behind the gap swayed gently and independently, without organized motions (Figure 6d).
At high vegetation density (ϕ = 0.133) and moderate Cauchy number (Ca = 16.88) in Case 6, gentle swaying of vegetation blades was observed in the vegetation canopy.Shear instabilities at the top of the vegetation canopy were evident as the presence of positive-negative pairs in the pressure field (Figure 7c).At later stage (t = 24s), the paring of vortex billow can be observed (Figure 7b), leading to the growth of billow size (Figure 7d).For the gentle swaying case, the wavy pattern at the interface was not significant compared to the cases with "monami".For cases 10 and 12, the Cauchy number with lower vegetation density, the wavy interface was more pronounced compared to Case 6 and, hence, they are classified as "Sway/Monami A" in Table 1.
At relatively low vegetation density (ϕ = 0.08) and moderate Cauchy number (Ca = 30) in Case 13, the deflection of the vegetation blade by the ambient current was more pronounced.A wavy interface was also evident at the top of the canopy due to shear instabilities.Kelvin-Helmholtz (KH) instabilities developed at the interface, and the billow size increased downstream due to the paring mechanisms (Figure 8).The vegetation blades adjusted  postures with the passage of the billow, and the wavelength of the interface was controlled by the size of KH billows, which increased downstream.The well-organized coherent billows led to a coherent wavy interface, known as "monami" phenomena.This pattern is categorized as "monami" A.
At high vegetation density (ϕ = 0.133) and high Cauchy number (Ca = 187.5),a wavy interface was evident in Case 5 (Figure 9).The vortex generated at the leading edge was primarily controlled by vortex shedding from the leading vegetation blade (Figures 9a and 9c) instead of shear instabilities at the interface.This pattern is categorized as "monami" B, as the wavy interface is not primarily controlled by Kelvin-Helmholtz billows.The interfacial wave propagated downstream slowly (Figures 9c and 9d).At large vegetation density, the interstitial fluid pressure can effectively prevent the unphysical penetration of vegetation blades, and no collisions among vegetation blades were observed.
At low vegetation density (ϕ = 0.08) and high Cauchy number (Ca = 120) in Case 11, vegetation blades were in the permanently prone position (Figure 10) and vegetation blades were laid on top of each other.Strong direct physical contact occurred as the interstitial fluid pressure was not strong enough to push two touching blades away to avoid unphysical penetration.Due to the formation of the shear layer (velocity profiles with two layer structure in Figures 10a and 10b) from the blockage of flow by the vegetation canopy, Kelvin-Helmholtz instabilities also developed at the top of the vegetation canopy, and the interface also adjusted with the propagation of KH billows downstream.However, the vortices formed at the top of the canopy were less coherent compared to "monami" A cases.In addition, the amplitude of the interfacial wave was smaller in prone cases because the overlying blades prevent the vertical momentum transfer of the fluid.

Blade Postures
Flexible vegetation blades interact with ambient flow and adjust their postures.When large turbulent eddies form at the top of the flexible vegetation canopy, spatial variations in vegetation posture can reflect the coherent structure of turbulence.To better understand the reconfiguration of the vegetation blade in a finite patch of submerged flexible vegetation, we analyzed the time history of the vertical tip position of several selected vegetation blades.Three blades were selected (Figure 11a) in Case 7 with high vegetation density (ϕ = 0.133 and Ca = 4.22).Blade 8 was located near the leading edge of the vegetation patch.The vegetation blade first got depressed by the   incoming flow, then it restored the tip position at t = 0.8s.Blade 8 slowly bent downward between t = 0.8 and 4s again by the ambient current.After t = 4s, vibrations of the blade were evident with a small amplitude.The vegetation blades near the leading edge in the convergence region (see Figure 6b) stayed together and acted like a single vegetation blade.Vegetation blade 16 was located near the gap with trapped vortices.Similar to blade 8, it bent down and restored its posture almost back to the initial position (2s < t < 4s), then it slowly reconfigured its posture toward the flow direction (4s < t < 10s) with the free-end tip approaching an asymptotic position of z tip ≈ 0.098 m.At the equilibrium stage, blade 16 swayed back and forth gently.Vegetation blade 24 was located near the trailing edge of the canopy.After it restored from the depression at t = 0.8s, the free-end tip gradually approached the asymptotic value around 0.099 m.At the equilibrium stage, the free-end tip motion of blade 24 was not well-organized, and the amplitude of the oscillatory motion was larger than that of blade 16 due to the growth of shear layer.
For Case 5 (Ca = 187.5 and ϕ = 0.133), interfacial waves propagated slowly downstream ("monami" B). Figure 11b shows the time history of the vertical tip position (z tip ) of three selected vegetation blades.Blade 8 was located near the leading edge of the vegetation canopy.Periodic oscillation can be observed for t > 9 s, with a period around 2 s.Several vegetation blades form a cluster and behave like a bluff body under the current.The vortex shedding from the vegetation blade (the first blade at the leading edge) induced the oscillatory motion of the vegetation blade near the leading edge.Similarly, the oscillatory motion was also evident in blade 27 near the trailing edge, which was induced by the vortex shedding from the cluster of vegetation blades at the trailing edge.Vegetation blade 16 was located at the center of the vegetation patch, and the large intermittent oscillation in z tip can be observed (18s < t < 21s) due to the propagation of the interfacial wave.
In Case 6 (ϕ = 0.133 and Ca = 16.88),gently swaying motion of the vegetation blades has been observed (Figure 7).Near the leading edge, the blade showed a periodic vibration at the equilibrium stage (Blade 8 in Figure 12a).The amplitude of the vibration was small compared to the swaying motion of the blades away from the leading edge.A peak around 6.4 Hz can be identified in Figure 12b, which corresponds to the natural frequency of the elastic blade.For blades 8 and 16, oscillations in the tip position were evident for t > 10s.The amplitude of the oscillation in beam 24 near the trailing edge was greater than that of beam 16 at the center of the patch.A peak around 1.6 Hz showed up for both beam 16 and beam 24, which corresponds to the frequency of the coherent turbulent eddies formed at the top of the canopy (Figure 7).
In Case 13 (ϕ = 0.08 and Ca = 30.0),"monami" occurred with well-organized turbulent eddies (Kelvin-Helmholtz billows) formed at the top of the canopy.After 4s, large oscillations can be observed in all three blades (Figure 13a).The power spectrum of z tip showed similar patterns.A peak around 2 Hz showed up for all three blades, which corresponds to the Kelvin-Helmholtz instability.This value is much higher than 0.2-0.3Hz in previous experiments by Okamoto et al. (2016) and 0.5 Hz in high-fidelity model simulation by Wang et al. (2022).Due to the high blockage ratio in 2D simulation, the boundary layer thickness (δ) of the shear layer is around 0.9 cm, and the velocity difference ΔU = 0.54 m/s.For the shear layer, the frequency of the Kelvin-Helmholtz instability (Ho & Huerre, 1984) is given by f KH = 0.032ΔU/δ = 1.92 Hz.Model results agree reasonably well with the empirical formula.Another strong peak appears around 4 Hz, and a minor peak shows up around 6 Hz, which are the high-order harmonics.The appearance of the high-order harmonics indicates strong non-linear interactions among the flow and flexible vegetation blades.In addition, a power-law relation between the power spectrum and frequency can be observed at f ≥ 5 Hz.
Compared to Case 13, at a high Cauchy number (Ca = 120.0), the vegetation blades were in the permanent prone position in Case 11.The free-end tip of the vegetation blade did not move freely like in "monami" or swaying cases (Figure 13b) because the vegetation blades cluster at the top of the canopy.At the equilibrium stage, all blades had a similar deflected height (after 20 s in Figure 13b), with small variations caused by the passage of KH billows at the top of the canopy.The power spectrums of z tip for the three selected vegetation blades show a similar pattern; however, we did not observe the strong principle frequency and high-order harmonics as in Case 13.Also, at high frequency, a slow decaying tail in the power spectrum can be observed for all three blades.
Figure 14 shows the flow regimes as a function of vegetation density (ϕ) and Cauchy number (Ca).At low vegetation density, we observed the transition from sway to "monami" to prone at ϕ = 0.08 when the Cauchy number increases.The transition from sway to "monami" A is smooth, as the motion of vegetation blades is mainly controlled by Kelvin-Helmholtz billows formed at the top of the vegetation canopy in both regimes.The Cauchy number at which the transition occurs increases with the vegetation density.At large vegetation density (ϕ ≥ 0.133), "monami" B occurs, in which the wavy pattern in vegetation canopy was not induced by Kelvin-Helmholtz instabilities at the fluid-vegetation interface.In "monami" B, several vegetation blades form clusters and move collectively when the interfacial wave propagates downstream slowly.Prone cases occur at relatively low vegetation density with a large Cauchy number.At high vegetation density (ϕ = 0.2), "monami" B occurs even at large Cauchy numbers around 500.

Hydrodynamic Interactions Among Vegetation Blades
For a patch of submerged flexible vegetation, the hydrodynamic interactions among vegetation blades play an important role in the blade posture and the deflected height.Adjacent vegetation blades interact through the interstitial fluid pressure in the gap between the blades, in to the direct physical contact (or collision).With hydrodynamic interactions, the collective behaviors of vegetation blades in a patch can be significantly different from that of a single vegetation blade.He et al. (2022) also observed that there were no blade-to-blade interactions.They attributed it to the force spreading at the fluid-solid interface in the immersed boundary method, which is mediated by the fluid to prevent direct contact among vegetation blades.
To investigate the importance of the hydrodynamic interactions, a series of simplified simulations with two flexible vegetation blades mounted vertically in a laminar channel were conducted.Therefore, no turbulent closure model was implemented.The length of the vegetation blade was set to l = 0.1 m, and the thickness of the vegetation blade was set to b = 4 mm.The Young's modulus was E = 5 × 10 5 N/m 2 .Zhang et al. (2020) showed the interspace l/l (defined as the gap distance between two blades d normalized by the length of the blade l) in the dual-filament system plays a critical role in the dynamic behaviors of the system.They mainly focused on the parameter range of 0.5 ≤ d/l ≤ 4. Also, the material density of the beam was much greater than that of the fluid.In our study, we focused on smaller values of l/l = 0.15, 0.25 and 0.5, respectively.The deflected heights of the two vegetation blades were compared with that of a single blade.A parabolic velocity profile was prescribed at the inlet given by u(z) = U 0 z(H − z)/H 2 , in which H is the half depth of the channel (Figure 16 window).The velocity U 0 was set to 0.4 m/s for all cases.The Cauchy number   =  2 0  3 ∕ 3 = 50 .The material density of the vegetation was set to be the same as the fluid.The half-depth of the channel H was set to be the same as l.No slip boundary condition was implemented at both the top and bottom walls, and outflow boundary condition was applied at the outlet.The outflow boundary condition is a mixed type boundary condition to prevent the reverse flow back into the computational domain.The computational domain was 20l, sufficiently large to ensure that the outflow boundary condition did not affect the interior flow around the vegetation blades.
Figure 15 shows the vorticity and pressure field at t = 30 s for all four cases.Compared to the case with a single beam, the deflected height is always larger in cases with two beams.The two beams did not interact through direct collision; however, the trailing beam bent in the flow direction by the fluid pressure in the gap between the two beams (Figures 15c-15h).The deflected height depends on the ratio of d/l.When the ratio is small to intermediate (d/l = 0.15 and 0.25), large vortices with scales comparable to the length of the beam can be observed with a strong recirculation zone.When the ratio is large (d/l = 0.5), large vortical structures are still evident; however, the secondary vortex is weaker compared to the case with smaller d/l.
The deflected height increases with the decreasing d/l ratio (Figure 16) for both beams, represented by the vertical tip position of the beam.The deflected height of the leading beam depends more strongly on the ratio d/l than  that of the trailing beam.At a large d/l ratio of 0.5, the deflected height of the leading beam is close to the one with a single beam, suggesting that at low vegetation density, the leading beam behaves similarly to an isolated single beam.The deflected height of the leading beam is greatest at low d/l = 0.15, which is around 15% greater than that of a single beam.For the trailing beam, the deflected height at quasi-equilibrium does not vary much with d/l.The deflected height increases by around 16%-19% compared to that of a single beam.For d/l of 0.25 and 0.5, oscillations in the tip position of the trailing beam can be observed due to the vortex shedding.The mean deflected height of the two beams is larger than that of the single beam.In a patch of submerged vegetation, the deflected height of the canopy can be enhanced by the hydrodynamic interactions among vegetation blades and hence depends on the vegetation density.

Three-Dimensional Flow Structures and Blade Motions
The fully coupled fluid-structure model can capture both the shear instabilities and the interactions of Kelvin-Helmholtz billows with the vegetation canopy.The most notable limitation of the present work is the restriction to two dimensions.Compared to three-dimensional flow, the blockage effect in two-dimensional simulation is significant because the out-of-plan flow can pass through the gaps between vegetation blades.When Kelvin-Helmholtz instabilities occur at the top of the vegetation canopy, the breakup of the 2D rollers and the redistribution of energy in the lateral direction due to turbulence are important turbulent generation mechanisms in three-dimensional simulations.Two-dimensional simulations tend to over-predict the vortical structures due to the lack of three-dimensional turbulence generation and hence turbulent dissipation.
On the other hand, 2D simulations are widely used to study hydrodynamic instabilities since the most unstable (or fastest growing) modes are often two-dimensional (Charru, 2011;Ho & Huerre, 1984).Shear instability at the flow-vegetation interface is one of the main mechanisms for the formation of interfacial patterns, making 2D simulation an ideal tool to identify different regimes in flow over flexible submerged vegetation canopy (Ni et al., 2023).Although the flow over submerged vegetation canopy tends to be highly three-dimensional due to turbulence (Nepf, 2012), the two-dimensional simplification can still provide valuable insight on the hydrodynamics and fundamental mechanisms in flow over submerged vegetation canopy (Singh et al., 2016;O' Connor & Revell, 2019).

Conclusion
An immersed boundary-finite element with soft-body dynamics has been implemented to study steady flow over a finite patch of submerged flexible aquatic vegetation.The finite element ANCF model and the fluid-structure interaction model have been extensively validated.The FSI model can resolve the flow interactions with flexible vegetation and hence the reconfiguration of vegetation blades to ambient flow.The soft-body dynamics model can capture the direct physical contacts among vegetation blades to prevent unphysical penetration.
Five different flow patterns have been identified, which depend on the vegetation density (ϕ) and Cauchy number (Ca).At high vegetation density and low Ca, a converging cluster of vegetation blades forms at the leading edge of the canopy, with trapped vortices separating the cluster from the rest of the canopy.Gently sway occurred when the Cauchy number was moderate.At high vegetation density and high Ca, a wavy interface formed ("monami" B) and propagated slowly downstream, which was not caused by Kelvin-Helmholtz (KH) billows due to shear instabilities.At low to moderate vegetation density, "monami" occurred at moderate Cauchy number ("monami" A), in which well-organized KH billows formed at the top of the canopy, causing the wavy interface.At low to moderate vegetation density, when Cauchy number was large, vegetation blades were in the permanent prone position.KH billows still formed at the top of the canopy; however, the vortices were less organized compared to the "monami" A cases.
At high vegetation density, vegetation blades interacted with adjacent blades without direct physical contacts.The hydrodynamic interactions via pressure can effectively push approaching vegetation blades away from each other and prevent unphysical penetration.However, at relatively low vegetation density, a soft-body dynamics model was necessary.With hydrodynamic interactions and direct physical contacts, vegetation blades in the canopy behaved collectively, and hence the deflected height was larger than the deflection height of a single blade with the same inflow velocity.In a flexible aquatic vegetation canopy, due to blade-to-blade interactions, the posture and hence deflected height of vegetation blades will depend on both vegetation density and Cauchy number.
Corrections to empirical/semi-empirical formulas based on a single vegetation stem may be needed when applied to predict the deflected height of a patch of vegetation.

Figure 2 .
Figure 2. (a) Mask function and AABB bounding box.(b) Collision detection and response using DEM.Black boxes are AABBs for the first blade, and blue boxes are AABBs for the second blade.
test was conducted.Because the material density of the vegetation is close to the density of water, we assume vegetation blades are neutrally buoyant and neglected the buoyancy force in this study.Only the FSI3 case was implemented to validate the model.The computational domain in the 2D benchmark case is 2.5 m by 0.41 m.The center of the circular cylinder is positioned at C = (0.2, 0.2) measured from the left bottom corner of the channel, and the radius of the circle is 0.05 m.The elastic beam has a length of 0.35 m and a thickness of 0.02 m.The tip is positioned at (0.6, 0.19), and the left end is fully attached to the circular cylinder.The density of the fluid and the beam were set to 10 3 kg/m 3 , and the kinematic viscosity of the fluid was 1 × 10 −3 m 2 /s.A large value of the kinematic viscosity was used in the benchmark test to keep the flow laminar.A parabolic velocity profile was implemented at the inlet with  () = 1.5 (ℎ−) (∕2) 2 and H = 0.41 m.The mean inflow velocity   was 2 m/s, and the maximum inflow velocity was 1.5   .The inflow was introduced into the computational domain slowly as a function of time (t) following

Figure 3 .
Figure 3. (a) Large deflection cantilever beam by a point load.(b)Model Validation with analytical solution.

Figure 11 .
Figure 11.Time history of the tip position z tip of three beams from (a) Case 7 and (b) Case 5.

Figure 12 .
Figure 12.(a) Time history of the tip position z tip of three beams from Case 6.(b) The corresponding energy spectrum of the tip position z tip .

Figure 13 .
Figure 13.Time history of the tip position z tip of three beams from (a) Case 13 with "monami" and (b) Case 11 in prone position.The corresponding energy spectrum of the tip position z tip from (c) Case 13 and (d) Case 11.

Figure 14 .
Figure 14.Flow pattern as function of vegetation density (ϕ) and Cauchy number (Ca).

Figure 16 .
Figure 16.Time history of the tip position z tip .