Numerical Modeling of the Interaction of Solitary Waves and Submerged Breakwaters with Sharp Vertical Edges Using One-Dimensional Beji & Nadaoka Extended Boussinesq Equations

Using one-dimensional Beji & Nadaoka extended Boussinesq equation, a numerical study of solitary waves over submerged breakwaters has been conducted. Two different obstacles of rectangular as well as circular geometries over the seabed inside a channel have been considered in view of solitary waves passing by. Since these bars possess sharp vertical edges, they cannot directly be modeled by Boussinesq equations. Thus, sharply sloped lines over a short span have replaced the vertical sides, and the interactions of waves including reflection, transmission, and dispersion over the seabed with circular and rectangular shapes during the propagation have been investigated. In this numerical simulation, finite element scheme has been used for spatial discretization. Linear elements along with linear interpolation functions have been utilized for velocity components and the water surface elevation. For time integration, a fourth-order Adams-Bashforth-Moulton predictor-corrector method has been applied. Results indicate that neglecting the vertical edges and ignoring the vortex shedding would have minimal effect on the propagating waves and reflected waves with weak nonlinearity.


Introduction
Boussinesq type equations are among the most practical mathematical models used in offshore engineering.These equations include nonlinear terms as well as dispersion terms.Thus, they are one of the most robust tools for hydrodynamic study of nearshore waves.During the years from 1871 to 1872, Boussinesq introduced these equations by adding dispersion effects to the shallow water equations originally known as Saint Venant.These equations have a hyperbolic structure with derivatives of high order in order to numerically model the dispersion-based physics.Peregrine [1] introduced what is known as the basic type of Boussinesq equations.Using Boussinesq equations for inviscid fluids, continuity equation with integral representation and applying respective boundary conditions, the basic Peregrine-Boussinesq equations for long waves over variable seabeds can be derived.Many efforts have been made for the development of Boussinesq equations.These efforts have been made with the aim of enhancing dispersion (ratio of water depth to wave length) characteristics of the equations in order to preserve their validity for deep waters applications.As one of the first attempts, Witting [2] used the momentum equations on the integrated depth in one dimension so he could present practical equations based on the velocity terms that were defined on the free surface.In the governing equations, Taylor series were used to represent velocity components.Coefficients of the expansion were derived in a way to achieve the best agreement with the linear dispersive properties.Later on, Madsen et al. [3] developed the dispersive properties of these equations to include deep water conditions by obtaining resulting terms based on the International Journal of Oceanography assumption of long waves in shallow waters and adding them to the classic equations of Abbott et al. [4].Other forms of the extended Boussinesq equations were derived by Nwogu [5] where the velocity components in an arbitrary distance from the calm free surface were used.These equations were derived directly from 2-dimensional continuity and Euler equations on a seabed with variable depth.Another type of extended Boussinesq equations was obtained by Beji and Nadaoka [6] using an algebraic manipulation of the classical Peregrine equations.Extended Boussinesq equations carry similar nonlinear and dispersive properties.This means that, while their algebraic representation differs from each other, the dispersive properties of these equations could be indicated to be equivalent after performing a respective change of variable.
Various physical phenomena have been simulated so far using the extended Boussinesq equations for many practical applications.Simulation of undular bore that was carried out by Peregrine [7] can be noted as the first physical modeling using Boussinesq equations in which a second-order accurate finite difference method with one corrector step for the free surface was employed.Almost twenty years later, Abbott et al. [4] introduced the first physical model for the regular wave shoaling on a sloping beach using the Boussinesq equations.They used a finite difference scheme for spatial terms of the equation in order to cancel out the nonphysical dispersion sources.Later, Schäffer et al. [8] and Nwogu [5] modeled the regular and irregular waves propagation on the sloping bottom using the extended Boussinesq equations.In the first attempt, Schäffer et al. [8] applied the breaking effects to the equations when the wave nonlinearity was being increased and the wave was getting closer to breaking up conditions.Moreover, experiments of Beji and Battjes [9] that were carried out for the regular wave propagation over a submerged trapezoidal bar are among the most applicable practical experiments in this field.They used the Boussinesq equations of Abbott et al. [4] while introducing a predictor-corrector scheme in order to numerically model their experiments.
Using a finite difference method, Dingemans [10] compared the classical and extended Boussinesq equations with each other in the case of regular wave propagation over a submerged breakwater.In the meantime, with the introduction of a fourth-order predictor-corrector method for time integration and spatial discretization of extended Boussinesq equations of Nwogu [5], Wei and Kirby [11] made efforts to decrease the truncation errors to achieve a fourth-order accurate scheme.They applied a multitude of physical models such as random wave shoaling on a slope, 2-dimensional sloshing wave evolution in a basin, and regular wave propagation over shoal in order to verify the accuracy of their numerical method.Ohyama et al. [12] introduced a finite difference method for Nwogu's Boussinesq equation [5] and modeled the regular wave propagation over a bar.In another attempt, Wei et al. [13] applied Nwogu's Boussinesq equation [5] for the simulation of solitary wave propagation on a sloping beach.Ambrosi and Quartapelle [14] applied the modified Taylor-Galerkin finite element method to the classical Peregrine equations in order to model the solitary wave propagation and also its interaction with a cylinder.Using Galerkin finite element method with a complex temporal scheme (known as Sprint), Walkley and Berzins [15] used the Nwogu's extended Boussinesq equations for the numerical simulation of regular wave propagation over a submerged breakwater and the 1-dimensional solitary wave propagation over a beach attached to the coastal wall.Wave propagation in portal regions over regular geometries was also modeled by Li et al. [16] using Beji and Nadaoka extended Boussinesq equations while applying Galerkin Finite Element Method.Walkley and Berzins [17] extended their numerical method to consider 2-dimensional case of wave propagation in portal regions over irregular geometries.In the past decade, Sørensen et al. [18] simulated the 2-dimensional physical phenomenon of surf zone using unstructured meshes.Nonlinear oscillations and harbor resonance problem were investigated by Woo and Liu [19] using a finite element method on Nwogu's extended Boussinesq equations.In another attempt, the problem of sloshing in a closed tank was modeled by Lin and Man [20] using Nwogu's extended Boussinesq equations with the application of a finite difference method.In the recent years, Ghadimi et al. [21] simulated solitary wave shoaling over coastal beds using Beji and Nadaoka's extended Boussinesq equations and with the application of a finite element method.They [21] presented a reasonable estimation of wave breaking using their numerical technique.More recently, Liu et al. [22] modeled solitary wave run-up in a cylinder group using a finite element method approach on Beji and Nadaoka's extended Boussinesq equations.
Although there seemed to be a wide range of applications for the extended Boussinesq equations with the above mentioned studies being only a brief overview on their implementations, the use of these equations for modeling submerged bars and breakwaters is rarely reported as in Ting et al. [23] and Yao et al. [24].
In the present study using the FEM approach of Ghadimi et al. [21] for the numerical solution of Beji and Nadaoka's extended Boussinesq equations, interaction of solitary wave with submerged bars for nondiffracting waves condition is investigated.Since the geometry of the submerged bars includes sharp vertical edges, it has been shown that these geometries cannot be directly implemented using the current numerical method.Therefore, a set of geometrical techniques are studied and then applied for the problem at hand.

Governing Equations
In the present work, Beji and Nadaoka's extended Boussinesq equations are employed.These equations include various terms of free surface, (, , ), integrated velocity vector in depth, u = (, V), depth of the seabed profile which is measured from the calm free-surface level, (, ), and gravity acceleration, , which are presented as follows: International Journal of Oceanography where  is a free variable for which the value of 1/5 is proposed by Beji and Nadaoka [6].Nonlinearity and dispersive properties are introduced by parameters  and , respectively.As the wave approaches the beach, the ratio of wave amplitude to water depth (/), which is the source of nonlinearity, would increase and the ratio of water depth to wave length (/) will decrease.For Beji and Nadaoka's extended Boussinesq equations with the consideration of the above mentioned value for , the error of the dispersion relation with  ≈ 0.5 would be less than 0.5 percent.It is worth mentioning that the accumulative error of the equations for  ≈ 0.17 would be less than 2 percent.

Numerical Method
The FEM technique of Ghadimi et al. [25] is employed for the numerical solution of the governing equations in 1dimensional form.The Galerkin FEM using linearly interpolated elements and functions is used for the discretization of the spatial terms.In order to couple the discretized form of the temporal and spatial terms, a fourth-order Adams-Bashforth-Moulton predictor-corrector method has been applied.Since the present method is 2nd-order accurate in terms of spatial discretization, Von-Neumann stability analysis has been employed.Therefore, the aforementioned scheme would be stable for   = √(Δ/Δ) < 0.5, where   is the courant number.

Initial and Boundary Conditions.
In the present study, the computational domain is a rectangular channel with the length and depth of  and , respectively.A solitary wave of amplitude of  is propagated from the left side of this channel at the time  = 0.This solitary wave is defined with profile velocity given as follows: in which  is the propagation velocity of the solitary wave.
A fully reflecting boundary condition is imposed on the exit boundary, and the normal velocity on the boundary is assumed to be zero: Here, n is the unit normal vector on the reflecting boundary.

Solitary Wave Propagation over the Submerged Bars
In the present study, the propagation of solitary wave over two types of submerged bars is simulated.The first simulation considers a, submerged bar.This geometry has vertical edges and is impossible to be implemented in the model directly using the present method.Therefore, a set of polygonal approximations with inclined sides have been employed instead of the vertical edges, and the results are presented.The second submerged bar is an attached step to the wall which also has vertical edges and, therefore, requires further modifications.
4.1.Solitary Wave Propagation over a Semicircular Bar.For the case of solitary wave propagation over a semicircular submerged bar, the experimental investigations of Cooker et al. [26] are employed.They performed a set of experiments to investigate the interaction of solitary waves with amplitude  over a semicircular submerged bar of radius  and introduced a relation between the reflected and transmitted wave for breaking and nonbreaking conditions.Since a breaking model is not employed in the present method, the submerged bar of radius  = 0.6 m and a solitary wave of amplitude  = 0.311 m are used for modeling the interaction of solitary wave with the seabed depicted in Figure 1 (in the deep side of a tank  = 1 m).
The scale of the wave nonlinearity just before reaching the submerged bar in the fixed-depth channel is  = / = 0.311 which is far from the breaking conditions.When the wave reaches the semicircular submerged bar due to the decrease in depth, it is expected for the wave amplitude to elevate.In the region close to the semicircular bar, the nonlinearity based on the experimental results would be equal to  = 0.4.Once again, this case is safely away from the breaking conditions, and during its propagation, the solitary wave will theoretically exhibit a nonbreaking procedure.Due to the decrease in water depth which occurs when the wave passes over the submerged bar crest, wave would be dispersed, and the wave crest, which was sharpened due to shoaling, requires strict considerations.As stated before, the semicircular submerged bar shown in Figure 1 cannot be numerically modeled using the present FEM approach.Since the computational domain is discretized using finite elements, the neighboring nodes of an element laying on the convex semicircular surface would have discrete depths and as the wave gets sharpened and the free surface gets accelerated, a discrete profile would be obtained over time from the numerical results (Figure 2).As a remedy to this numerical phenomenon, it is required to approximate the convex surface with a regular or irregular polygon (triangle or quadrilateral) in order to have more than two adjusting nodes on each side of the slope.In fact, there should exist a minimum of one intermediate node on each slope so that elements on that side would be discretized more than once (Figure 3).
In the present study, three gauges have been employed in order to record the time series of the propagating solitary wave with two of them being placed before the submerged bar with distances of 0.9 m and 2.5 m from the bar, while another one is placed 0.4 m behind the bar.Four numerical experiments have been arranged with four different polygonal approximations of the semicircular bar (Figure 4), and the resulting time series of the shoaling numerical solution are compared with the gauge results.In the present numerical method, the computational domain inside the channel has a length of  = 40 m and is discretized using 1000 linear finite elements.The time step for the numerical solution is set equal to Δ = 0.01. Figure 4 represents four different models employed for the approximation of semicircular submerged bar that are used in the present numerical simulation.The approximated geometries are set in a way that they are all inscribed in a semicircle with a radius of 0.6 m.Model 1 represents an inscribed triangle, and model 2 is a regular semioctagon in which the intersection of the middle sides forms a 45 degrees angle with respect to the center of the geometry.Other models include geometries with a 30-and 60-degree angle between the intersections of the middle sides.Therefore, the distributions of the resulting finite element nodes on each side are shown in Figure 4.The calculated time series for the models introduced earlier are shown in Figures 5-8.
As depicted in these figures, numerical and experimental results have a slight difference between them in terms of phase angle which could be attributed to the irrotational flow assumption for the Boussinesq governing equations and also the proposed approximation of the semicircular submerged bar.Numerical results recorded by the first gauge (from left) give two peaks in the time series for all four models.The first peak is due to the approaching wave crest while the second one denotes the reflecting wave crest.Based on the results of the experiments, a distinct crest for the reflecting wave cannot be distinguished, while both experimental and numerical results demonstrate a similar time range for the reflecting wave.The second gauge shows the moment at which the wave passes the first slope (left side) of the submerged bar and therefore wave shoaling occurs.An overestimation is observed in the results of the second gauge which shows its relevance to the results of the first gauge.The celerity of the reflecting wave using the Boussinesq model is lower than that of the empirical theories and this causes an overestimation in the free surface elevation that was recorded by the second gauge.In the meantime, a third gauge records wave dispersion and shows two peaks.After crossing the submerged bar, the wave starts to disperse and therefore the second reflection occurs; the idea that is, supported by the second peak recorded in the time series.Contrary to the first gauge, the second peak recorded by the third gauge cannot be observed.In fact, the second reflection occurs sooner in the Boussinesq model.After crossing the submerged bar, the wave amplitude using fourth approximation model (Model 4) is modeled more precisely.It is shown that only a triangular approximation of the semicircular submerged bar would give favorable results for the free surface elevation over the bar when the Beji and Nadaoka extended Boussinesq equations are applied.This is due to the fact that the wave-front would be accelerated while crossing over the bar which would cause instabilities for the case of the models 2, 3, and 4.This phenomenon can be observed as nonphysical dispersions that are shown on Figures 6, 7, and 8.That is why   some deviations are observed between the present results and the ones obtained from the experiments.Mentioned cases are the main reason of difference in numerical and experimental results.It is worth mentioning that the main goal of the present study was to have a suitable numerical modeling of the physical problem at hand using the Boussinesq equations due to the fact that their computational speed exceeds their rival alternatives in fluid motion equations.Time history of the propagating wave over Model 1 (Triangular Approximation) is in Figure 9 for time of 0.6 seconds.

Solitary Wave Propagation over a
Step Attached to the Coastal Wall.Solitary waves are stable during their propagation over a fixed depth, and this characteristic makes them more destructive.Therefore, a set of vertical concrete caissons is generally used in front of the coastal walls in order to act as a breakwater at times of low tide.These costal steps would   become unstable when collide with these kinds of solitary waves and would be destructive over time.On the other hand, their presence in the coastal line for the prevention of the wave impacts on the coastal wall is of high importance.Therefore in this part of the present study, propagating solitary wave over a step, that is, attached to the coastal wall, is modeled using the Beji and Nadaoka extended Boussinesq equations.Experimental results of Grilli et al. [27] are used for the validation purposes.In their experimental model [27], they introduced a solitary wave with nondimensional amplitude of  = / = 0.33 in the deep side of a tank ( = 1 m) while the wave would approach a step with a depth of ℎ 1 = 0.67 m.The wave profile is recorded over propagation using different gauges.Using the present numerical model, the introduced problem is simulated, and the free surface elevations are compared with the experimental results in Figure 10.Due to the fact that the vertical face of the step in this model causes a discontinuity between the two different depths, a gradual slope is used in a way that enables the required continuity of the finite elements in this region.The obtained numerical solution is in a good agreement with the experimental results for the respective gauges.

Conclusion
In the present paper, Beji and Nadaoka extended Boussinesq equations have been applied for numerical modeling of solitary waves over submerged bars with sharp vertical edges using finite element method.A semicircular submerged bar along with a stepped seabed was implemented.Since these bars have vertical a discontinuity occurs in depth which makes it impossible have direct implementation of the present numerical method.Thus, a gradual slope was employed instead of the vertical edges.In the first type of submerged bar, inscribed polygons in the semicircle were used instead of the semicircular submerged obstacle.Irregular geometries were chosen in such a way that the intersection of the middle sides makes a 30-, 45-and degree models, respectively.Numerical results for all four models demonstrate that in these experimental cases, the Boussinesq model simulates the collision causing the first reflection of solitary wave to occur a bit sooner and the second reflecting wave appears a bit later.However, due to being more efficient compared to other equations used in simulating fluid motions, the present Boussinesq model shows to be a suitable criterion for the present experimented phenomenon.Due to the fact that the wave front gets sharpened after crossing over the submerged bar, implementation of sequential polygonal slopes would cause instabilities for the dispersed wave.As a result, only the triangular approximation of the semicircular bar is proved to be acceptable in the present numerical model.Moreover, the solitary wave propagation over a step attached to the coastal wall was investigated, and favorable agreement of the numerical results with the experimental data was achieved.It is worth mentioning that, due to the discontinuity problem that was discussed earlier, the step was replaced by a gradual slope in the second test case.

Figure 1 :Figure 2 :
Figure 1: Physical modeling of the interaction of a solitary wave and a semicircular submerged bar.

Figure 3 :Figure 4 :
Figure3: Approximation of nodes height on a semicircular bar so that a minimum of two elements are defined on each slope.

Figure 5 :
Figure 5: Time series of the gauges compared to the experimental results for Model 1.

Figure 6 :
Figure 6: Time series of the gauges compared to the experimental results for Model 2.

Figure 7 :
Figure 7: Time series of the gauges compared to the experimental results for Model

Figure 8 :
Figure 8: Time series of the gauges compared to the experimental results for Model 4.