Numerical investigation of nonlinear wave passing through finite circular array of slender cylinders

ABSTRACT Cylinders extend from the plat bottom of an open channel piercing the water free surface to mimic a patch of emergent vegetation and are often used as models in numerical simulations or experiments related to flow in the presence of vegetation. This work focuses on investigating the flow field structures induced by nonlinear waves passing through a cylinder patch with a diameter of DG and different solid volume fractions . A 3D non-hydrostatic numerical model is used to simulate the nonlinear wave passing through the emergent circular array of slender cylinders. The instantaneous vortical structures within the cylinder patch are clearly predicted by the numerical model. The residual flow reveals abundant flow characteristics of a long term-averaging wave motion. The 3D numerical model effectively presents the spatial and temporal wave process evolution, and contributes to highlighting vegetation flows under nonlinear waves.


Introduction
In estuarine and coastal reaches, vegetation plays a significant role in the attenuation of storm surges and waves. It has been established that coastal forests can serve as barriers against tides, storm surges and tsunami waves (Tanaka, Sasaki, Mowjood, Jinadasa, & Homchuen, 2007). In vegetation flow investigations, the rigid vegetation is often modeled by solid vertical cylinders, by the effective approximation of the stems.
Waves interacting with solid bodies are often studied in ocean engineering; for example, the floating bodies supported by columns. An analytical model provides a powerful means of revealing certain critical parameters and identifying characteristic wave motions, including wave diffraction, attenuation and near trapping (Evans & Porter, 1997;Linton & Evans, 1990;Malenica, Taylor, & Huang, 1999;Maniar & Newman, 1997;Meylan & Taylor, 2009). Moreover, the analytical method has been extended to handle a large number of cylinders (Kashiwagi, 2017). The analytical method uses linear wave theory and neglects the viscous effect, which is applicable to cases with a small Keulegan-Carpenter (KC) number. For wave propagation through vegetation, the analytical method is of little practical use, because the viscous effect is more important and the number of vegetation stems is large.
CONTACT Jingxin Zhang zhangjingxin@sjtu.edu.cn In recent decades, real-life case studies using computational fluid dynamics have been promoted (Akbarian et al., 2018;Fotovatikhah et al., 2018;Mou, He, Zhao, & Chau, 2017;Wu & Chau, 2006). For the vegetation flow studies, numerical models have also been promoted to simulate the flow through vegetation patches. The numerical models consider eddy viscosity by solving the Navier-Stokes equations. Using the VOF method, nonlinear water waves can be simulated accurately. The computational cost of modeling flows through vegetation is often enormous, because of the large number of rigid stems and high grid resolution required to mimic solid obstructions. In order to save on computational costs, a porous model is a popular approach in the numerical simulation of waves through vegetation patches (Gao, Falconer, & Lin, 2011;Li & Yan, 2007;Li & Zeng, 2009;Li & Zhang, 2010;Ma, Kirby, Su, Figlus, & Shi, 2013;Maza, Lara, & Losada, 2013). In the porous model, solid obstructions are modeled by adding a bulk drag force to the momentum equations. Although the porous model is highly efficient in modeling vegetation flows, the simulated flow structures within the vegetation patch are not accurate. Compared to the porous numerical model, the direct numerical simulation of a flow passing through the solid cylinder patch is feasible for visualizing turbulent flow structures. Cui and Neary (2008) applied large eddy simulation (LES) to investigate the flow in partially of fully vegetated channels containing a uniform layer of vegetation, by means of adding source terms into the momentum equations to model the drag force. Okamoto and Nezu (2010) promoted an LES model for investigating the turbulent flow and dissolved mass advection and diffusion in fully vegetated channels, in which the non-slip solid boundary condition was used to model the vegetation with a high-resolution grid. Nicolle and Eames (2011) used 2D cylinders to model a finite porous obstruction, focusing on the local and global effects of an isolated group of cylinders. Baranya, Olsen, Stoesser, and Sturm (2012) used a 3D RANS model based on nested grids to simulate the flow around circular piers. Compared to the direct modeling of a unidirectional flow through cylinders, simulations of waves through cylinders are rarely reported. Large amplitude of free surface elevation and high resolution of the solid wall boundary flow simulation are the main challenges of numerical simulating waves through a patch with a large number of cylinders. Although the direct simulation is lower efficient than the porous model, it is highly contributed to the investigation of the flow field within the cylinder patch, for example, calculating the drag force on individual cylinders.
In this study, a full 3D numerical model was used for the direct simulation of a nonlinear water wave passing through a finite group of cylinders with different solid volume fraction (SVF) values ( ). The model was developed by introducing non-hydrostatic pressure into the hydrostatic model; that is, the non-hydrostatic model (Casulli, 1999;Casulli & Zanolli, 2002;Fringer, Gerritsen, & Street, 2006;Zhang, Sukhodolov, & Liu, 2014). The model is highly efficient in capturing the free surface, and lower on the computational cost by means of the pressure splitting scheme. This work focuses on the highlights of the 3D flow structures induced by waves through a finite circular array of slender cylinders. The eddy viscosity was considered and the solid obstructions were accurately modeled by a non-slip solid wall boundary condition on the cylinder surface.

Numerical model
Compared with the shallow water equation model with the hydrostatic assumption, a non-hydrostatic component p n is introduced into the momentum equations beyond the hydrostatic component p h = ρg(ζ − z). The governing equations are derived in the σ coordinate, which is highly efficient in simulating the moving water elevation and the spatial varying bed in the natural large-scale fluid flows (Chau & Jiang, 2002, 2004. The z coordinate is transformed to σ coordinate using σ = (z − ζ )/(h + ζ ). The coordinate transformation is sketched in Figure 1. The 3D governing equations are expressed in a conservative form: where g is the gravitational acceleration, ζ is the free surface elevation, h is the still water depth, D = h + ζ is the total water depth and Q is the source intensity for wave generation. The velocity u, v and w are replaced by variables. The velocity in the vertical σ coordinate is obtained by the following formula: The S-A model (Spalart & Allmaras, 1994) is used to calculate the eddy viscosity coefficient υ t . The working variableυ is calculated by solving the following transport equation: and the eddy viscosity υ t is then calculated by where f v1 = χ 3 /(χ 3 + c 3 v1 ), and χ ≡υ/υ is the aspect ratio ofυ to the molecular viscosity υ. The parameter f w is calculated as The magnitude of the strain rate tensor is calculated as |S| = 2 ij ij with ij = 1 2 (∂u i /∂x j − ∂u j /∂x i ). The constants in the above equations are c b1 = 0.1355,σ = 2/3, c b2 = 0.622, κ = 0.41, c w1 = c b1 /κ 2 + (1 + c b2 )/σ , c w2 = 0.3, c w3 = 2.0, and c v1 = 7.1.

Wave generation method
Using the source term to generate water waves, two wave trains propagate in opposite directions from the source site (illustrated in Figure 2). Taking the x-axis as the wave direction, the source intensity is calculated as where U(x s , z, t) is the horizontal velocity of flows, x s is the adding source term site and x is the grid size where the source term is fixed. In this work, the velocities under the third order Stokes water wave are substituted into Equation (8) to generate the waves, which are expressed as Equation (9): where k is the wave number, D is the total water depth and σ is the vertical transformed coordinate. The coefficients are calculated as F 1 = ka/2 sinh kh, F 2 = 3(ka) 2 /4sinh 4 kh and F 3 = 3(ka) 3 [11 − 2 cosh 2kh]/ 64sinh 7 kh. The wave height a is determined by using the relation between a and the wave height H 0 , which is expressed as follows: In Equation (9), the wave speed C is determined by the following dispersion equation: The phase function θ in Equation (9) is calculated as In the numerical model, a sponge zone is used to damp the waves approaching the flume ends. A linear damping factor μ(x) is applied: where x 0 is the beginning location of the sponge region and x l is the sponge region length. In the sponge regions, the fluid velocities are damped by multiplying with the factor provided in Equation (13) at every numerical simulation step.

Numerical scheme
The non-hydrostatic model is solved by means of a splitting strategy of the pressure terms (Zhang et al., 2014). The finite volume method (FVM) was used to discretize the governing equations based on unstructured grids. A second-order total variation diminishing scheme (OSHER scheme) is adopted to discretize the advection terms. The in-house codes are paralleled by the OpenMP library.

Geometry of cylinder array
The array with a diameter D G was selected to have the same circular geometry with a diameter D S as its constituent components. The SVF was set by means of changing the number of cylinders . The cylinder number was determined by N c = (D G /D S ) 2 . The cylinders were uniformly fixed in the array in a staggered pattern. The distance D 0 of two neighbor cylinders changed with the varying SVF. All of the parameters are tabulated in Table 1.
The incident wave was generated in front of the cylinder array. Approaching the inlet and outlet of the open channel, the sponge zone was used to dampen incident waves. The channel width was W = 14D G or 500D S , while the length was L = 22D G or 800D S . In order to reduce the waves reflected against the side walls, the channel width was specified as sufficiently large to ensure the diffracted waves around the cylinders propagated into the sponge zone before arriving at the side walls.

Simulation implementation
To directly simulate the flow through cylinders, a fine grid resolution is needed to capture the flow separation induced by the solid obstruction. Within the cylinder array, the flow is affected by the complex geometry. The characteristic spatial scales involve the cylinder diameter D S , gap scale D 0 and cylinder arrangement pattern. The grid scale was designed to be of 0.1D S on the cylinder surface. The boundary layer was meshed in a sub-region with one-D S thickness adjacent to one cylinder. The mesh scale in the normal cylinder surface direction was varied from 0.05 to 0.12D S , with a stretching ratio of 1.1. Beyond the cylinder array, the grid scale was enlarged from 0.1 to 0.2 h. At a high SVF, the wake pattern began to exhibit properties similar to the wake flow induced by a solid cylinder, where the diameter D G becomes one critical length scale apart from the local water depth. A schematic of the computational mesh is illustrated in Figure 3. The incident wave was generated at approximately 6D G upstream of the cylinder patch. The numerical wave flume setup is sketched in Figure 4. The parameters of the incident third-order Stokes wave are listed in Table 2.

Wave modeling
The wave generation method was first validated by the numerical simulation of nonlinear waves propagating over a flat bottom. In the numerical flume, the still water depth was h = 0.05m and the target traveling wave height was H 0 = 0.01m. The horizontal space was discretized with x = 0.01m and the vertical space was discretized with a relative coordinate σ = 0.004 ∼ 0.1, where the bottom boundary layer was meshed using the finest grids. The historical simulated water elevation at three times the wave length, far from the wave generation site, was recorded and compared to the analytical third order Stokes wave. Figure 5 illustrates the wave modeling validation. The calculated water elevation coincided strongly with the analytical results.

Unidirectional flow through cylinders
In order to demonstrate the applicability of the turbulence model, the experiments on the unidirectional flow through a finite patch of cylinders, as carried out by Zong and Nepf (2011), were replicated by the numerical model. The experiments were conducted in a laboratory flume under uniform flow conditions. The case used for validation had the following parameters: flow depth h = 0.133m and free-stream velocity U 0 = 0.098m/s. The vegetation was simulated by rigid slender cylinders with diameters of D s = 0.006m, at bulk SVFs of = 0.03 Figure 6. Calculated non-dimensional distances + n of the nearest grid points around the most upstream cylinder at the middepth. and 0.1. The cylinder patch geometry in this study was the same as in the experiments. Figure 6 shows the grid resolution near the wall, which is selected at the most upstream cylinder. Figures 7 and 8 illustrate a comparison of velocities between the measurements and simulations. The mean velocities (ū,v) simulated agreed well with the measurements. Nevertheless, the predicted fluctuation velocities (v rms ) deviated gradually from the measurements along the transverse direction. The simulated velocity fluctuation was obviously weaker than the experimental measurements outside the wake (about y/D > 0.5), which indicated the background velocity fluctuation being poorly modeled by the RANS model. It is predicted to improve the simulations by increasing the turbulence intensity at the inlet or adopting advanced numerical model, e.g. LES.
The cylinder patch geometry for simulating the wave through the cylinders was the same as that of the unidirectional flows. The orbital velocity amplitude listed in Table 2 approximated the uniform velocity in the experiments. The maximum Reynolds number of the oscillated flow around the cylinders was approximately equal to that of the unidirectional flow. It was reasonably assumed that the grid resolution for mimicking the solid cylinders in the unidirectional flow simulation was achievable for the oscillated flow simulation.

Wave force on cylinders
The flow field structure of a wave wrapping a single cylinder can be quantitatively identified by the KC number. The characteristic flow structure varies with the critical KC number (Venugopal, Varyani, & Barltrop, 2006).   At a small KC number ( < 4), the fluid viscosity can be neglected, and no flow separation or vortex shedding are created strongly on the solid wall. At a high KC number, the fluid viscosity effect must be considered in the numerical model in order to capture the flow separation and vortex shedding. In the present study cases, the KC number of wave wrapping for a single slender cylinder was approximately 6 for a shorter incident wave (W 1 ) and approximately 10 for a longer incident wave (W 2 ). The viscosity model is necessary for simulating the wave passing through the cylinder patch.
Compared to the wave wrapping a single cylinder, the wave process is significantly more complicated within the cylinder patch. Statistical analysis is a powerful tool for identifying the complicated flow constituents. The temporal drag and lift force on every slender cylinder were recorded and used to analyze the hydrodynamic force. The instantaneous total drag and lift force integrated through the water depth were normalized by a characteristic inertial force. The dynamic force coefficients were normalized by the maximum orbital velocity of the incident wave at the free surface, the diameter of a single cylinder and the still water depth; that is, C Di = F xi /0.5ρU 2 0 hD S and C Li = F yi /0.5ρU 2 0 hD S . In the statistical analysis, only the forces on the cylinder at the array center were calculated and analyzed. Figure 9 illustrates the temporal records of the drag coefficient and frequency spectrum in study cases C 1 W 1 ∼ C 2 W 2 . The disturbance of the C Di values was distinctly observed in study cases C 1 W 2 and C 2 W 2 compared to those in study cases C 1 W 1 and C 2 W 1 . The stronger flow separation at a higher KC number was determined to contribute to the drag force disturbance. The relevant frequency spectra of the C Di coefficients are displayed in Figure 9 for the different study cases. The first frequency peak is equal to that of the incident wave, while the second peak corresponds to the flow separation induced by the nonlinear waves.
As opposed to the drag force, the lift force is mainly induced by the alternate vortex shedding from the cylinder surface. The vortex shedding generation and evolution are confined to both the incident wave conditions and array geometry. Figure 10 illustrates the temporal records of the lift force coefficients of the same cylinder and the relevant frequency spectrum. In the study case C 1 W 1 , the lift coefficient C Li approximates that of a single cylinder owing to the low SVF (C 1 ). At a high SVF (C 2 ), the lift coefficient C Li reveals stronger disturbances owing to the intense flow interaction affected by dense solid obstructions. Apart from the SVF, the incident wave condition is critical to identifying the lift force characteristics. The non-dimensional parameter KC number is generally used to quantize the strength of the inertial force against the viscous force. The KC number of the single cylinder in case C 1 W 1 was approximately 6, in which the flow separation at the solid wall was not very strong for inducing the spanwise wave loads. With an increase in the KC number (KC ≈ 10 for condition W 2 ), the vortex shedding was distinctly formed, resulting in a higher C Li amplitude in case C 1 W 2 compared to case C 1 W 1 . However, this deduction cannot be used for a higher SVF because the vortex is severely suppressed by the bleed flow through the cylinder patch (Nicolle & Eames, 2011). The lift force constituents can be clearly identified by means of frequency analysis. Figure 10 illustrates the frequency spectra of all of the study cases. Apart from study case C 1 W 1 , where the vortex shedding is suppressed, the main frequency peak was shifted to twice that of the first peak. The main frequency to wave frequency ratio appeared within a reasonable range of [1.64, 3.07] in previous research (Moberg, 1988) for wave loads on a single cylinder with a moderate KC number. The frequency peak amplitude for the high KC number (W 2 ) cases was larger than that for the lower KC number (W 1 ) cases. The frequency analysis again clarified the critical KC number for identifying the vortex shedding generation with a wave wrapping the vertical cylinders. A modulation of the lift coefficient C Li was observed in cases C 1 W 2 and C 2 W 2 . The analysis revealed the vortex shedding being strengthen with the KC number increasing, i.e. the wave length increasing. The vortex shedding by the bulk of cylinders was deduced to be stronger in cases with a longer incident wave length (W 2 ). Owing to the frequency characteristics of a uniform flow around a single cylinder, i.e. a characteristic Strouhal number at 0.2, the lowest frequency in the spectrum justly approximated to the case of a similar incident flow around a single cylinder with the present patch diameter. It was reasonable to deduce that the modulation of the lift coefficient C Li was induced by the vortex shedding from the bulk of the cylinder patch. At a high SVF, the porous obstruction approaches to a solid one, in which case the vortex shedding is expected to be stronger. Figures 11 and 12 reveal the spatial distribution of the averaged C Di and C Li coefficients in order to visualize the heterogeneity of the wave loads on cylinders within the patch. The root-mean-square (RMS) of the C Di and C Li where N is the sample number of the temporal records. Figure 11 reveals that, at a low SVF, the mean drag force increased through the cylinder patch and reached peak values near the patch rear. With the increase in SVF, the drag force reached maximum values over a shorter penetration length into the patch, and then gradually decreased downstream. The comparison of the drag force within the patch clarified the SVF effect on the wave dampening. The results reveal that the C Di values under longer incident waves were frequently smaller than those under shorter incident waves for a constant SVF. The smaller C Di values indicated a lower wave dampening efficiency or stronger wave transmission through the cylinder patch. Figure 12 illustrates the C Li distribution within the patch for different SVF and incident waves. A larger C Li indicates a stronger flow separation, which results in intense vortex shedding. Under longer waves (W 2 ); that is, a larger KC number, the C Li values are larger than those under shorter waves (W 1 ).

Flow field through the cylinder array
The flow field visualization provides direct insights into the wave process through the cylinder array. Figure 13(a) shows an instantaneous snapshot of the water elevation around the cylinder patch. In the cylinder array lee, Figure 11. RMS of drag force coefficients C Di . partially resonant waves were formed where the crest line was nearly perpendicular to the incident wave direction. The resonant wave length increased downstream, resulting in a triangle-like wake field. Despite different SVFs and incident waves, the resonant wave trains were similar in the wake field. Within the cylinder array, the flow structures were distinctly affected by the SVFs and incident waves. Figure 13(b,c) reveal the instantaneous wave forms within the patch for cases with W 1 and W 2 , respectively. The instant wave fields were used to analyze the instantaneous hydrodynamic characteristics. The instantaneous second invariants of the velocity gradient tensor   (Q) were used to visualize the turbulent flow structures. Figure 14 illustrates the instantaneous isosurface of Q ( = 0.5), where the isosurface was flooded by the vorticity magnitude. With an increasing wave length (higher KC number), the separation flows were significantly more intense, resulting in larger vortical structures. With the increase of SVF, the separated vortex was distinctly suppressed within the cylinder array for both wave conditions. Under these wave conditions, the vortex motion was distinctly limited in the upper water body. Near the bottom, no distinct horseshoe vortex was formed around the cylinders.
Furthermore, the instantaneous streamlines were used to represent the spatial and temporal evolution of the wave process through the patch. Figure 15 illustrates the instantaneous streamlines at the middle water depth. The flow violently fluctuated where the velocity direction turned around (Figures 15(a-d)) within the array. In the cases of shorter waves (W 1 ), two velocity turning zones were common within the array scope, which outlined a limited space for the separated flow evolution. In the cases of longer waves (W 2 ), only one velocity turning zone appeared within the cylinder patch. The separated flow and shedding vortex developed significantly more adequately during the longer half-wave cycle. In the study cases, the instantaneous flow structures highlighted the critical identification of the vortical structures according to KC number. With the increasing SVF, the denser solid obstruction destructed the large turbulent flow structures into smaller ones.
The flow signatures can be visually represented by the instantaneous vorticity field. Figure 16 shows the instantaneous vertical vorticity at the middle depth. For uniform flows through cylinders, vortex shedding is characteristic within the cylinder patch at a low SVF, even with a sufficiently low Re number (approximately 100) (Nicolle & Eames, 2011). The characteristic Reynolds number calculated for the maximum orbital velocity in the study cases was approximately 600, which was sufficient for inducing vortex shedding in uniform flows. However, a high Reynolds number flow cannot be sustained for the oscillated flows, which induces a non-fully developed separated flow to generate an intense shedding vortex. For a larger KC number, the separated flow was developed more fully, and vortex shedding could be observed significantly (Figures 16(b,d)). The wake signatures in the near and far field are characteristic coherent lumps of positive and negative vorticity under uniform incident flow conditions (Nicolle & Eames, 2011). However, the vorticity signature does not survive when the wave passes through the cylinder array and is rapidly attenuated.
The mean flow structures represent the residual flow of the oscillated flows passing through the cylinder array. The residual flow is important for analyzing the longterm mass transportation forced by oscillated flows. Figures 171819-20 reveal the 3D characteristics of the mean wave process through the cylinder array, in which the upper plane is at σ = −0.05, the middle plane is at σ = −0.5 and the lower plane is at σ = −0.95. Figure 17 illustrates the mean streamlines adjacent to the array center for C 1 W 1 at different water depths. Near the free surface, four steady circulations were formed around the single cylinder located on the central line, but only two distinct steady circulations were observed around the cylinders deviating from the central line. The steady circulation pattern was clear at the middle water depth; however, the flow structures were not symmetric about the cylinder. The steady circulation structure was significantly smaller in the plane near the bottom, and appeared only in the down wave cylinder face because of the skew of the oscillated flow approaching the bottom. Figure 18 displays a similar but enlarged steady circulation structure through the water depth compared to those illustrated in Figure 17. For case C 1 W 2 , the incident wave length was approximately twice that for C 1 W 1 . The large circulation structures indicate that the separated flows developed adequately at a higher KC number.    The mean flow structures at a high SVF are illustrated in Figure 19 for shorter waves (W 1 ) and Figure 20 for longer waves (W 2 ). Compared to C 1 W 1 , the steady circulations around the cylinders were far more distinct for C 2 W 1 (Figure 19). With the increasing incident wave length, the steady circulations between the two cylinders merged into one larger circulation (Figures 20(a,b)). However, the similar mean flow structures were not observed near the bottom, indicating that the bed friction was predominant in the bottom boundary layer. The mean flows were not uniform through the water depth, particularly in the mean flow direction. For shorter waves, the mean flow was directed upstream near the free surface, but downstream at the middle water depth and near the bottom for both SVFs (Figures 17  and 19). For longer waves, the mean flow was directed upstream until the middle water depth, and changed to downstream when approaching the bottom for both SVFs (Figures 18 and 20). The mean flow direction indicated a net mass transportation path within the array under oscillated flows.

Summary
A 3D numerical model was used to investigate nonlinear waves passing through a finite circular cylinder patch. The present model is verified to be highly efficient in simulating waves through a large number of cylinders using a high grid resolution of the solid wall boundary layer. By quantifying the streamlines, vorticity field and instantaneous flow velocities, the simulation provides a clear illustration of the wave motion within and around the patch. The hydrodynamic force analysis reveals a clear frequency spectrum. The wave loads on a single cylinder within the patch are dominated by two main vortex systems, one is formed by the separation flow from the cylinder wall, and the other one is composited by the vortex shedding induced by the bulk of the cylinders. The bulk vortex shedding means the flow separation induced by an obstruction with the same patch diameter but with a porosity. The wave load characteristics directly relate to the KC number and the cylinder patch porosity. The mean wave loads on cylinders present a spatial distribution within the patch depending on the SVF and the incident wave conditions. Analysis of the mean flow, i.e. the residual wave motion, reveals abundant spatial flow structures induced by the nonlinear waves through the cylinder patch. The mean flow reveals a long-term pathway for the mass transportation forced by waves. The net flow flux varies through the water depth, which is distinctly different from the uniform flow through cylinders. The flow structures are characterized by both the SVF and KC numbers.
Although the present work highlighted full spatial hydrodynamic process when non-linear waves propagating through a patch of cylinders, the model scale was limited in the scope of experimental model rather than the prototype in natural environment. In actual field reality, the vegetation patch geometry is much more complex, and the hydrodynamic conditions are also much more complicated. Meanwhile, the computational cost on direct simulating vegetation stems is very huge. The numerical model is expected to be applied to the studies of the vegetation flow and related environmental issues in natural rivers, lakes and estuaries.