Simulations of moving effect of coastal vegetation on tsunami damping

A coupled wave–vegetation simulation is presented for the moving effect of the coastal vegetation on tsunami wave height damping. The problem is idealized by solitary wave propagation on a group of emergent cylinders. The numerical model is based on general Reynolds-averaged Navier–Stokes equations with renormalization group turbulent closure model by using volume of fluid technique. The general moving object (GMO) model developed in computational fluid dynamics (CFD) code Flow-3D is applied to simulate the coupled motion of vegetation with wave dynamically. The damping of wave height and the turbulent kinetic energy along moving and stationary cylinders are discussed. The simulated results show that the damping of wave height and the turbulent kinetic energy by the moving cylinders are clearly less than by the stationary cylinders. The result implies that the wave decay by the coastal vegetation may be overestimated if the vegetation was represented as stationary state.

Many numerical and experimental approaches have been developed in recent years to help understanding the tsunami wave interactions with coastal vegetation. The coastal vegetation was idealized by a group of rigid cylinders in most investigations. Huang et al. (2011) performed both experiments and a numerical model by considering solitary wave propagating on emergent rigid vegetation and found that dense vegetation may reduce the wave transmission because of the increased wave reflection and energy dissipation into turbulence in vegetation. By using both direct simulation and 5 macroscopic approach, Maza et al. (2015) simulated numerically the interaction of solitary waves with emergent rigid cylinders based on the arrangement of laboratory experiments of Huang et al. (2011). Previous approaches (e.g. Anderson et al., 2011;Huang et al., 2011;Maza et al., 2015;Wu et al., 2016) assumed that the idealized mangrove vegetation is stationary and neglected the plant motion with the wave. In the reality, however, vegetation stems may move like cantilever or whip driven by waves (Paul et al. 2012). 10 For most of natural vegetation is deformable, which reduces flow resistance, thus it is robust way to be not neglecting vegetation motion. Accordingly, this study tries to provide a better physical model in the numerical simulation by considering vegetation motion coupled with waves to investigate wave damping performance by mangroves. A direct numerical model based on computational fluid dynamics (CFD) is adopted in this paper for simulating the wave damping characteristics including both the stationary and moving vegetation. 15

Numerical model description
Among a number of open source CFD codes available, IHFORM (Higuera et al., 2013;2014) is specially designed for coastal engineering applications. IHFORM model was used in Maza et al. (2015) for direct simulation of the solitary wave interacting with the stationary vegetation. Alternatively, the model Flow-3D (Flow Science, Inc., 2012) is applied in this paper to conduct the numerical simulations to involve the motion of the vegetation accompanied by wave. Flow-3D provides 20 exclusively the FAVOR (fractional area/volumes obstacle representation) technique (Hirt, 1993) and the general moving object (GMO) model that is capable of simulating the rigid body motion dynamically coupled with fluid flow. The FAVOR technique retains rectangular elements with a simple Cartesian grid system and shown to be one of the most efficient methods to treat the immersed solid bodies (Xiao, 1999). The free water surface tracking in the model is accomplished by using volume of fluid (VOF) method (Hirt and Nichols, 1981). 25 Referring to previous literature, the problem is idealized by solitary wave passing on a group of emergent rigid cylinders.
Considering the fluid to be incompressible, the continuity and momentum equations for a moving object formulated with area and volume fraction functions are given as Nat. Hazards Earth Syst. Sci. Discuss., doi:10.5194/nhess-2016-353, 2016 Manuscript under review for journal Nat. Hazards Earth Syst. Sci. Published: 23 December 2016 c Author(s) 2016. CC-BY 3.0 License.
where VF is the fractional volume open to the flow and A is the fraction area for the subscript direction, ui is the mean velocity component in the ith direction, the subscripts = 1, 2, 3 represent x-, y-and zdirections, respectively, p is the pressure intensity, ρ is the fluid density, g is the gravitational acceleration, µ is the absolute viscosity, t ν is the eddy viscosity, k is the turbulent kinetic energy, and ij δ is the Kronecker delta function such that ij δ = 1 when i = j; ij δ = 0, when i ≠ j. The eddy viscosity t ν is related to the effect of the space and time distribution of the turbulent 5 motion, which is solved here by using the renormalization group method (RNG k-ε model) proposed in Yokhot and Orszag (1986). It is noted that the above governing equations are rendered to standard RANS equations as both VF and A are set to unity.
Comparing with the continuity equation for stationary obstacle problems, / (1) is equivalent to an additional volume source term and exists only in mesh cells around the moving object boundary. It can be calculated using 10 where cell V is volume of a mesh cell, obj S , n  and obj V  are respectively surface area, unit normal vector and velocity of the moving object in the mesh cell. The relative transport equation for the VOF function F is given using According to kinematics, general motion of a rigid body can be divided into a translational motion and a rotational motion. 15 Equations of motion governing the two separate motions are where T, J, and ω  are total torque, moment of inertia and angular acceleration about the fixed axis. And the velocity of any point P on the rotating cylinder is calculated by denotes distance from the fixed end C of the cylinder to point P.
In computing the coupling of fluid and rigid body interaction, the velocity and pressure of fluid flow are first solved. The hydrodynamics forces on the rigid body are then obtained and used to calculate the velocity of the rigid body. Then the 5 volume and area fractions are updated according to the new position of the rigid body, and the source term can be calculated using equation (3). The flow field is computed repeatedly until the convergence is achieved.
As for the boundary conditions, the normal stress is in equilibrium with the atmospheric pressure while shearing stress is zero on the free surface. All of the solid surfaces were treated using the no-slip boundary condition. The variation of the turbulent energy and the turbulent energy dissipation on the free surface boundary was set as zero in the normal direction. 10 The solution of solitary wave derived from Boussinesq equations was employed as the incident wave. Two different uniform computational meshes around the cylinder field, 0.002 m and 0.001 m respectively, were used to test the numerical accuracy and the sensitivity to grid size. Fig. 2 shows that FAVOR technique resolved successfully the geometry of cylinders using these two computational grids constructed. Fig. 3 shows the comparison of free surface evolution between the present numerical results and experimental measurements for an incident wave height Hi = 0.05 m.

Validation
The results obtained by the direct simulation using IHFORM in Maza et al. (2015) were also shown in the figure. Fig. 4  25 shows the maximum wave height at each wave gauge probe between the numerical results and experimental measurements.
The results depict that the present numerical simulations are in a very good agreement with the laboratory experiments. The results also show almost no differences between both computational meshes used.

Results and discussion
The above comparisons demonstrated the present numerical model is capable of simulating accurately the wave evolution  respectively as the solitary wave crest passing through gauges G3 to G6 for an incident wave height Hi = 0.05 m. For the stationary cylinders, it can be observed that the water velocity reduces rapidly as wave crest impinging on the array of cylinders. The moving cylinder has angular rotation and immerses in the water while the wave crest passing over, thus 15 resulting less water velocity reduction than the stationary cylinders. Fig. 6 shows the comparison of the horizontal velocity profile as wave crest passing through each gauge for the moving and stationary cylinders, which confirms that the motion effect of the cylinder leads to less water velocity reduction.

Surface elevation evolution
The wave surface evolutions for the moving and stationary cylinders respectively are shown in Fig. 7. Before the solitary 20 wave reaches the cylinders, the wave height almost keeps with the same, and the surface evolution has no much difference between the stationary and moving cylinders. However, surface evolution appears difference while and even after the wave propagation over the stationary and moving cylinders. The weakly wave reflection can be found at the front row of the stationary cylinders while it is not obviously for the moving cylinders. The result of free surface evolution depicts that the stationary cylinders are working better than the moving cylinders for the wave height damping as the waves pass through the 25 cylinders. It can be seen that the wave height ratio (H/Hi) decays rapidly for the stationary cylinders but mildly for the moving cylinders. 30 As shown in Fig. 9, the maximum wave height damping ratio for the moving cylinders is approximately 26% yet it is nearly 61% for the stationary cylinders after the wave passed through the cylinders. Here the wave height damping ratio is

Turbulent kinetic energy dissipation evolution
The turbulent kinetic energy will be generated and dissipated during the wave interacting with the group of cylinders. The 5 turbulent kinetic energy dissipation is the important mechanism of the wave height damping. The turbulent kinetic energy (k) and dissipation of the turbulent kinetic energy ( ε ) are obtained from the RNG k-ε turbulent closure model while solving general RANS equations. Fig. 10 shows snapshots of the distribution of the turbulent kinetic energy dissipation (DTKE) along the cylinder array for the moving and stationary situations respectively when the wave crest passing over gauges G3 to G6 for an incident wave 10 height Hi = 0.05 m. It can be seen that the turbulent kinetic energy dissipation happens as the wave crest impinging on the cylinders. For stationary cylinders, the largest turbulent energy dissipation occurs in the region between gauges G3 and G4 as the wave crest passing through the gauge G4. But for the moving cylinders, the largest turbulent energy dissipation occurs as the wave crest passing through the gauge G5 and the dissipation is less than the stationary cylinders case.
The total DTKE shown in Fig. 11 is calculated from the total computed meshes of the numerical tank at the time when 15 wave crest passing each gauge. It shows that the total DTKE increases rapidly to a maximum value after wave crest passed the front cylinders. After having its maximum value at gauge G4 for the stationary cylinders and gauge G5 for the moving cylinders, the total DTKE decreased. It is reasonable to see that the larger incident wave height has the larger DTKE. It is also observed that the stationary cylinders have larger DTKE before the wave passing gauge 5 but the moving cylinders have larger DTKE after the wave passed gauge G5. This is because that a little turbulent kinetic energy could be generated and 20 dissipated by moving cylinders during the restoring process to the stand position after the wave passed over. Nevertheless, the integral DTKE for the moving cylinders is less than that of the stationary cylinders.

Conclusions
A numerical simulation based on the general RANS equations and RNG k-ε turbulent model was implemented to investigate the moving effect of coastal vegetation on the damping of tsunami wave. The vegetation was idealized by a group 25 of emergent, rigid cylinders. The FAVOR technique and moving object (GMO) model provided in Flow-3D code were employed in this paper for simulating the coupling of fluid and rigid body interaction. The evolutions of wave height, flow field and turbulent kinetic energy dissipation rate for both stationary and moving cylinder are investigated. Due to the moving effect of the cylinders accompanied by the wave, the numerical results showed that their wave height damping and the turbulent kinetic energy dissipation rate were markedly less than those of stationary cylinders. This result leads to an 30 Nat. Hazards Earth Syst. Sci. Discuss., doi:10.5194/nhess-2016-353, 2016 Manuscript under review for journal Nat. Hazards Earth Syst. Sci.  Vol. 4, No. 7, 1510-1520., 1992. 15 20 25 30 Nat. Hazards Earth Syst. Sci. Discuss., doi:10.5194/nhess-2016-353, 2016 Manuscript under review for journal Nat. Hazards Earth Syst. Sci.          Nat. Hazards Earth Syst. Sci. Discuss., doi:10.5194/nhess-2016-353, 2016 Manuscript under review for journal Nat. Hazards Earth Syst. Sci. Published: 23 December 2016 c Author(s) 2016. CC-BY 3.0 License.