Numerical investigation of free surface elevation and celerity of solitary waves passing over submerged trapezoidal breakwaters

In this paper, evolution of both free surface elevation and celerity of solitary wave interacting with submerged breakwater is numerically investigated by solving one-dimensional extended Boussinesq equations derived by Madsen and Sorensen. Spatial discretization is done by Galerkin finite element method (FEM) and for time integration, predictor-corrector method of Adams-Bashforth-Moulton is used. Propagation of solitary waves has been simulated over four different seabed slopes and computed results, compared against published work of Grilli et al. [17], indicate favorable agreement. Furthermore, a solitary wave with two different amplitudes is propagated over trapezoidal breakwater and evolution of free surface elevation is studied and validated. Finally, effects of side slopes of trapezoidal breakwater on wave characteristics has been investigated for 6 different slopes. It is clearly indicated that damping of the free surface elevation and celerity of the wave are strongly affected by the steepness of the breakwaters.


INTRODUCTION
In most of the engineering problems, mathematical models cannot be solved analytically and require a numerical solution.With an increase in computational technologies and because of the advantages of numerical methods over analytical solutions, many numerical models and software programs have been developed for various engineering practices.This is in addition to the fact that, more realistic models of greater complexity, can be investigated using numerical techniques.This makes it possible to predict the course of an event before it actually occurs, or to study various aspects of an event mathematically without actually running expensive and time-consuming experiments.
Simulation of solitary waves using different type of water equations has become the object of many industrial and scientific projects [1,2].Numerical and experimental works related to propagation of solitary waves mainly consist of transformation of a solitary wave over variable seabed topography and sloping beaches or propagating over a shelf, obstacles, barriers and submerged dikes.As an example, Santfo et al [3] investigated motion of solitary waves over triangular shapes in the context of shallow water wave approximations.Cooker et al [4] studied propagation of solitary waves over submerged semicircular cylinders for different cylinder size and wave heights.Based on the latter, Forbesl et al [5] investigated free surface flow over semicircular obstructions.Similar to the investigation of Huang and Dong [6] of the interaction of a solitary wave and a submerged dike, Ohava et al [7] modeled propagation of nonlinear waves over submerged dikes.Huang and Dong [8] have studied vortex generation and wave deformation.Experimental propagation of solitary waves over a rectangular dike is carried out by Lin et al [9].
In recent years, there has been much attention given to different types of Boussinesq equations for simulations of free surface problems due to their simplicity and reduced computational time which is achieved by integrating equations over the depth and eliminating one dimension in computational domain.Therefore, Boussinesq type equations would be an efficient alternative to other types of water equations such as incompressible Navier-Stokes equations.Boussinesq equations are capable of accurately simulating the water surface evolution due to the presence of different phenomena including shoaling, diffraction, refraction, and reflection in marine environment.
There have been many attempts to investigate Boussinesq type equations.For example, Ghadimi et al [10] by using extended Boussinesq equations of Beji and Nadaoka simulated transformation of periodic wave over sloping beaches.Based on these equations, Liu et al [11] also studied run up of solitary wave in a cylinder group.In other works Ghadimi et al [12,13] studied solitary wave transmission and regular wave propagation over submerged breakwater using Madsen and Sorensen extended Boussinesq equations.Also, Lin and Man [14] using Nwogu's extended Boussinesq equations, modeled sloshing in a closed tank.
Among previous works, there have been extensive efforts in simulating propagation of water waves over submerged obstacles implementing different types of water equation.However, studying not only free surface evolution but also celerity of solitary waves over trapezoidal breakwater in addition to the effect of side slopes of breakwater on these important characteristics using Boussinesq equations is rarely investigated.Therefore, the aim of the present work is to conduct this numerical investigation using Boussinesq equations derived by Madsen and Sorensen.Quality of the numerical approach is verified by applying the numerical model to a problem for which published numerical solution is known.

MATHEMATICAL FORMULATION
Boussinesq equations of Madsen and Sorensen [15] consist of continuity and momentum equations as presented in (1) (2) Equations ( 1) and ( 2) are related to continuity and momentum equations, respectively, where h(x) is still water depth, g is gravity acceleration, d is defined as d = h + S, S being the water surface level.In this equation, Q is depth-integrated velocity component.Subscript x denotes partial differentiation with respect to space.B is a dispersion coefficient which, based on the work by Madsen and Sorensen [15], a value 1/15 is selected.

NUMERICAL METHOD 3.1 FINITE ELEMENT MODELING OF GALERKIN
As stated earlier, spatial discretization of equations is accomplised by Galerkin finite element method.Solution of a complicated problem can be obtained by dividing the region of interest into small regions and approximating the solution over each sub-region by a simple function.Approximation of dependent variables by finite element method is presented in equation 3. ( where f i is the nodal value of dependent variables and N i is the standard basis function.

SOLUTION FORMULATION
By implementing linear Galerkin method, final forms of the equations are presented in equations 4 to 6. ( A detailed description of numerical method dedicated to resolution of an auxiliary equation for modeling Madsen and Sorensen's Boussinesq equations, has been presented in a previous work by the present authors [13].

SPATIAL DISCRETIZATION
By applying Galerkin method, the element matrix for continuity, momentum and auxiliary equations becomes Here, G q denotes boundary of the element with unit normal n.For Dirichlet boundary conditions, the boundary integrals are eliminated.

TIME DISCRETIZATION
To discretize the derivatives in time, the predictor-corrector method of Adams-Bashforth-Moulton is implemented.The assembled continuity and momentum global equations are written in matrix form in equation 15.
[M]{ f .} = E (15) At the predictor step, the amount of dependent variables can be obtained using Eq.16. ( Subsequently, fourth-order Adams-Moulton scheme is used to calculate final values of the variables at the corrector step, which is calculated from the predictor step. (17)

WAVE ABSORBING BOUNDARY CONDITIONS
There are several methods for dealing with unbounded domains in numerical simulations.Boundary conditions are a required component of the mathematical model.Different boundary conditions may cause quite different simulation results.Improper sets of boundary conditions may introduce nonphysical influences on the simulation system, while a proper set of boundary conditions can avoid that.Therefore, arranging the boundary conditions for different problems becomes very important.Here, a wave absorbing boundary is needed to ensure that no unwanted reflection occurs at boundaries of the solution domain.Wave damping can be achieved using several numerical methods which allow waves to exit the solution domain without reflection.To achieve this goal, damping method of sponge layers proposed by Larsen and Dancy [16] is adopted.In this method, surface elevation and fluxes are divided by a coefficient m(x), after each time step.m(x) takes the following form: Here, d is the distance between boundary and the sponge layer, d s is often equal to one or two wave lengths, Dd is size of the elements, and a is a specified constant.In the present work, a = 4 is the prescribed value.

NUMERICAL SIMULATION 4.1 SOLITARY WAVE SHOALING ON PLANE BEACHES
In this section, shoaling of solitary waves on four different slopes (S = 1:8, 1:15, 1:35 and 1:100) using Madsen and Sorensen's equations is considered.Initial surface profile and velocity can be obtained using equations 19-21. ( Here, c is wave celerity and α is wave amplitude.
The obtained results in various time steps are compared with the results of Grilli et al. [17].Computational domain and its parameters are illustrated in Figure 1.where the slope starts.Here, coefficient d is defined as .
Figure 2 shows the comparison of obtained results for nonbreaking solitary wave with the results of Grilli et al [17] for four different seabed slopes at three different wave heights (d).Based on the computed results of propagation of the solitary wave, one can see the increase of nonlinearity and asymmetry of the solitary wave due to its shoaling.Although slight under prediction of free surface profile is observed, the reported results seem to have good agreement with the work of Grilli et al for all slopes in different time steps.

PROPAGATION OF SOLITARY WAVE OVER A SUBMERGED BREAKWATER
In this section, propagation of solitary waves over submerged trapezoidal breakwaters is considered using one dimensional Boussinesq equations derived by Madsen and Sorensen.
Free surface elevations are compared with those by Grilli et al. [18] which are based on fully nonlinear potential flow equations.Computational domain is illustrated in Figure 3.For this particular test case, numerical domain is 40.0 m long and is discretized into 400 elements with a mesh size of Δx = 0.1m.With respect to the work of Grilli et al. [18], the water depth (h) is 1.0m and two incident waves of 0.06m and 0.1m heights are considered.The time step is fixed and equal to 0.005s.Slope of both sides of the breakwater (with 0.8m height), is 1:2.Slopes are symmetrical toward 0 and the corners of left slope are at x' 1 = -2.0mand x' 2 = -0.4m.Also, sponge layer is considered at the right side of the computational domain.Results of computations are presented at at which time the wave crest is located at .
Figures 4 and 5 show the numerical results of propagation of solitary waves over a submerged trapezoidal breakwater using Madsen and Sorensen Boussinesq equations and the work by Grilli et al [18] for different time steps.By comparing the obtained results, no major differences can be observed in evolution of free surface profile and all figures show high agreement of the results.
In Figures 4 and 5, first curve shows the wave profile before reaching the breakwater in a constant depth, which maintains its initial shape.The second curve is the crest of solitary wave which is starts to increase due to shoaling effect and reaches the maximum height in the third curve.Over the breakwater, crest exchange takes place and wave profile decomposes into two waves (as shown by the fourth curve).Eventually, in the last three curves, the transmitted wave continues its propagation down the beach and the reflected wave propagates backward.At the right side of the computational domain, the transmitted waves are well damped by the sponge layer.

EFFECT OF BREAKWATER'S SIDE SLOPE ON DEFORMATION OF FREE SURFACE ELEVATION AND WAVE CELERITY
In this test case, effect of breakwater's side slope on deformation of the free surface elevation is investigated.All wave and computational domain characteristics are set as in test case (4.2).The analysed slopes of the four different breakwaters are S = 1:1, 1:2, 1:4, 1:8 and 1:16.The geometry of these breakwaters are displayed in Figure 6.
An initial wave amplitude of A = 0.06 has been chosen.For all cases, the maximum wave crest heights (h) and its position (x) in different time steps, shown in Figures 4 and 5, have been computed and presented in Table 1.
As seen in Table 1, in each time probe, not only the elevation but also the position of the maximum wave crest is affected by side slope of the breakwater.Therefore, one can easily deduce that celerity of the wave crest is changed with respect to the breakwater slope.In Table 2, the average wave crest celerity is calculated and presented for each time interval.For a better observation of the effect of breakwater slope on the wave, wave crest elevation is illustrated versus time in Figure 7.
It is easily observed from Figure 7 that, up to time probe C, where the crest of wave reaches the breakwater top, the effects of different slopes on the wave are nearly similar.However, as wave crest passes above the breakwater, elevation of the wave crest is clearly better damped by the breakwaters with steeper slopes.Therefore, best option for cancelling solitary waves elevation among the tested breakwaters is the steepest one (slope 1:1).The celerity of the crest of the waves can also be ploted against time to show the phase changes caused by different breakwater slopes.This plot is shown in Figure 8.
As observed, there is almost no celerity differences between slopes before the crest reaches the top of breakwater.When crest passes the breakwater, celerity of the wave is strongly affected by the slope of the breakwater, and as clearly observed in Figure 8, the steeper breakwater gives higher celerity to the wave crest, causing it to pass more rapidly over the breakwater.One can deduce the fact that amplitude of the wave is damped by a conversion of energy which results in a higher celerity.

CONCLUSIONS
In this study, numerical simulation of nonlinear and dispersive Boussinesq type equations of Madsen and Sorensen was considered.By using linear Galerkin finite element method, spatial domain has been discretized by subdivision of the continuum into non-overlapping elements.Time integration has been performed using Adams-Bashforth-Moulton predictorcorrector method.In order to absorb the outgoing waves, a sponge layer is considered at the right side of all computational domains.In the first test case, propagation of solitary waves over four different seabed slopes has been investigated and results have been compared with published results presented by Grilli et al.It was shown that the obtained results have good agreement for all slopes in different time steps compared to the mentioned numerical work.In another test case, a solitary wave with two different amplitudes is propagated over trapezoidal breakwater and evolution of free surface elevation is studied and validated.Finally, effects of breakwater's side slope on wave characteristics has been investigated for 6 different slopes.It has been demonstrated that steeper breakwater causes a better damping on the elevation of the wave.It is further shown that wave celerity is strongly affected by the slope of the breakwater and steeper breakwaters cause more increase in celerity of the wave.

Figure 2 :
Figure 2: Comparison of obtained wave profiles by the present work against the work by Grilli et al. [17] for slopes of 1:8, 1:15, 1:35 and 1:100 shown in plots a, b, c and d, respectively.Wave heights are =0.2 in a, c, and d and = 0.3 in b Int. Jnl. of Multiphysics Volume 9 • Number 1 • 2015 67

Figure 8 :
Figure 8: Wave crest celerity variation for different side slopes

Figure 7 :
Figure 7: Comparaison of crest elevation at different time probes for different side slopes

Table 2 :
Average wave crest celerity over different time intervals for various side slopes

Table 1 :
Maximum crest elevation and position for various side slopes at different probe times