A Unique Finite Element Modeling of the Periodic Wave Transformation over Sloping and Barred Beaches by Beji and Nadaoka's Extended Boussinesq Equations

Abstract

This paper presents a numerical model based on one-dimensional Beji and Nadaoka's Extended Boussinesq equations for simulation of periodic wave shoaling and its decomposition over morphological beaches. A unique Galerkin finite element and Adams-Bashforth-Moulton predictor-corrector methods are employed for spatial and temporal discretization, respectively. For direct application of linear finite element method in spatial discretization, an auxiliary variable is hereby introduced, and a particular numerical scheme is offered to rewrite the equations in lower-order form. Stability of the suggested numerical method is also analyzed. Subsequently, in order to display the ability of the presented model, four different test cases are considered. In these test cases, dispersive and nonlinearity effects of the periodic waves over sloping beaches and barred beaches, which are the common coastal profiles, are investigated. Outputs are compared with other existing numerical and experimental data. Finally, it is concluded that the current model can be further developed to model any morphological development of coastal profiles.

1. Introductions

Coastal profile can vary considerably during a year or even during a single storm event. As a result of this change in profile, long shore bars can be formed and existing bars can be shifted over barred beaches or destroyed and constant slope beaches are built. The cross-shore sediment transport plays a significant role in the development of the coastal profiles. In this case, long shore bars will be formed on a sandy constant slope profile as a result of the interaction between the onshore and offshore transport flows in the nearshore zone [1]. In the meantime, the behavior of periodic waves over the sloped and barred beaches is different. The nonlinearity of the periodic wave travelling over beaches with constant slope profile continually increases following the same trend over the barred beaches. However, in the case of waves travelling over the barred beaches, after the waves pass over the bars they get decomposed and their nonlinearity decreases. The breaking and nonbreaking conditions may occur in each of the mentioned cases, and studying them may be significant in development of the coastal profile.

Several governing equations exist for simulation of waves over the coastal zone. Fully nonlinear potential flow (FNPF) and incompressible Navier-Stokes equation are common governing equations that are used for simulation of the periodic wave interactions in the coastal zone. While two- and three-dimensional Navier-Stokes equations provide accurate mathematical models for wave motions, the numerical analysis of such an actual problem would involve a huge computational domain with an unknown water surface elevation. Therefore, as an alternative, the popular shallow water equations are used to simulate the periodic wave interactions which have been developed rapidly in the last decades. By integrating the Navier-Stokes equation over depth, the three-dimensional and two-dimensional problems can be reduced to two-dimensional and one-dimensional computational problems, respectively. Shallow water equations substantially reduce the computational cost associated with the wave modeling when diffraction, refraction, shoaling, and harmonic interaction are of particular importance. Boussinesq equation is one type of shallow water equations which is derived from the Euler or Navier-Stokes equations. In deriving the Boussinesq equations, free surface of the water and the averaged velocities are approximated by solving the continuity and momentum equations, respectively. These equations are ordinarily applied for determining the water surface elevation inside the harbors and wave interactions in near the shore.

The improvement of linearized dispersion property of the Boussinesq equations has been the basis for all developments of these types of equations from the past to the current time. Any improvement in the dispersion characteristic of these equations would bring about a better accuracy in results associated with the deeper waters. Validity of the linearized dispersion property of the original Boussinesq equations [2, 3] is limited to actual shallow waters, that is, for the range of 𝑘ℎ≤0.76 (𝑘 is the wave number and ℎ is the water depth). With derivation of the developed Boussinesq equations [4–6], linearized dispersion property of the Boussinesq equations has drastically improved and they can even be applied to the relatively deep water, that is, for the range of 𝑘ℎ≤3. However, in recent years, by derivation of high-order Boussinesq equations [7–10], the applicable range of the Boussinesq equations have increased (6≤kh≤40), but nonetheless, the complexities involved in sophisticated mathematical form of high-order equations has made the numerical modeling and/or finding the accurate numerical solution very difficult.

In the present work, the one-dimensional Beji and Nadaoka extended Boussinesq equations [6] have been used as governing equations to model regular wave shoaling and its decomposition in nearshore zone. Needless to say that, at nearshore zone, the dispersive properties of the extended Boussinesq equations and high-order Boussinesq equations are the same. The Beji and Nadaoka extended Boussinesq equations (BNBE) is derived by a simple mathematical manipulation of Peregrine’s equations.

From the mathematical viewpoint, the system of BNBE contains third-order free surface derivatives in momentum equations, and also the form of BNBE is not so straightforward when employed for variable depth. For this reason, the BNBE have not been considered as frequent as the other equations such as Nwogu Boussinesq equations by the researchers. Furthermore, finite difference scheme has been applied in most of the numerical models of the time-dependent Boussinesq-type equations [3, 6, 11]. Although, the usage of finite difference approach is easy, the finite element method is more efficient. While Finite element method can be applied to discretize complex domain, the finite difference method does not have the capability of handling the irregular boundaries which are encountered in several coastal engineering problems. Another important characteristic of the finite element scheme is its superiority in applying a finer mesh for some regions of interest, which would indeed yield in the minimization of the number of required grid points. Therefore, one can state with relative certainty that the use of finite element method for solving Boussinesq-type equations is preferred over the finite difference method by researchers. For example, Do Carmo et al. [12] and Do Carmo et al. [13] presented numerical solutions for the Boussinesq equations based on finite element method, excluding the dispersion modification terms. Ghadimi et al. [14] applied Galerkin finite element method for solving the extended Boussinesq equations and calculated the solitary wave shoaling on plane beaches. Li et al. [15] applied BNBE for simulation of regular and solitary wave’s dispersion in nearshore zone. Their numerical method was based on a finite elements approach with special treatment of third-order spatial derivatives.

The current paper presents a new and particular technique based on finite element scheme to investigate the propagation of wave over sloping and barred beaches, where both nonlinearity and dispersion properties of the waves are accurately modeled. Finite element method used here is based on the rewritten equations in lower-order differentials form by introducing an auxiliary variable [15, 16]. The auxiliary variable is derived by introducing a novel form of related equations. By showing the good modeling ability of the wave properties, it is also demonstrated how capable the suggested scheme can be for modeling the morphological development. Accordingly, the most typical beach profiles such as sloping beaches and barred beaches are considered in this study.

2. Boussinesq Equations System

The physical system depicting the entire seabed and wave profile can be represented by a water depth 𝐻, a wavelength 𝜆, and a wave amplitude 𝑎. The nonlinearity and dispersion existing in the system are denoted by the ratios 𝜖 and 𝜎 as follows: 𝑎𝜖=𝐻,𝐻𝜎=𝜆. (1) As discussed in the literature, BNBE are produced by adding extra dispersive terms to Peregrine’s momentum equations [2]. Mathematical form of these equations in one-dimensional domain is as follows: 𝜂_𝑡+𝜕()𝑢𝜕𝑥𝐻𝑢=0,_𝑡+𝑢𝜕𝑢𝜕𝑥+𝑔𝜕𝜂=ℎℎ𝜕𝑥1+𝛽3𝜕²𝑢_𝑡𝜕𝑥²+𝜕ℎ𝜕𝑥𝜕𝑢_𝑡ℎ𝜕𝑥+𝛽𝑔ℎ3𝜕³𝜂𝜕𝑥³+𝜕ℎ𝜕𝜕𝑥²𝜂𝜕𝑥², (2)(3) where 𝛽 is a free parameter and 𝑢(𝑥,𝑡) is the depth-averaged velocity vector in BNBE system.

The value of free parameters in extended boussinesq equation systems are obtained from the best agreement achieved between the linearized and exact dispersion relations. After introducing a new boussinesq equations system, the equations are linearized and phase velocities (𝐶) or group velocities will be compared to the exact dispersion relation or Stokes’ first-order theory (𝐶²=𝑔ℎ(tanh𝑘ℎ)/(𝑘ℎ)). The linearized dispersion relation properties of the BNBE that can be found [6] are given as in 𝐶²()=𝑔ℎ1+𝛽𝑘ℎ²/31+1+𝛽(𝑘ℎ)²./3 (4) In this paper, the value of 𝛽=0.177 has been considered. This value is the best choice because the best agreement between the linearized dispersion relations (4) and the exact dispersion relations is achieved at this value, with a maximum error of less than 6%.

3. Numerical Method

3.1. Finite Element Formulation

Application of a linear finite element method from the view point of simplicity in formulation is similar to the application of a central finite difference method, but with this difference the finite element can also be applied to the unstructured mesh. Because of the fact that (3) includes a third-order derivative in space, we cannot directly apply the linear Galerkin finite element scheme to the system of equations. Linear elements can only be applied, when the weighted residual form of the problem contains second-order spatial gradients (second-order free surface derivatives in BNBE) in which case the linear finite element approximation cannot be applied. Therefore, special treatment of the second-order derivatives becomes necessary.

Walkley and Berzins [16] presented a new technique, based on the linear Galerkin finite element method, to solve the Nwogu Boussinesq equations (NBE). To overcome the difficulty associated with the third-order spatial derivatives in NBE, they introduced an auxiliary variable and rewrote equations system in its lower-order form. Applying this auxiliary equation, the high-order terms of Nwogu Boussinesq equations are eliminated and thus the solution of NBE by linear elements becomes possible. The form of the suggested auxiliary variable is based on the nature of the differential form of NBE. Therefore, this model is only limited to special cases where this auxiliary term is inherently embodied in the governing equations like NBE, and not equations like BNBE (e.g., [6]) that do not explicitly possess such a term.

However, in this study, a flexible method for choosing the auxiliary variable is introduced which is not restricted by the mathematical form of the governing equations. For spatially discretizing the equations, the standard Galerkin finite element scheme is utilized. The finite element formulation of the dependent variables is as follows: 𝑓≈^𝑛𝑑_𝑖=1𝑁_𝑗𝑓_𝑖, (5) where 𝑓_𝑖, 𝑁_𝑗, and 𝑛𝑑 are the values of any independent variable at the nodal points, the standard basis functions, and the number of nodes, respectively. To discretize the spatial domain Ω, nonoverlapping elements are applied to subdivide the continuum and subsequently each element and node is numbered. The presence of third-order spatial derivatives in momentum equation (3) makes the application of a standard linear interpolation finite element method for conversion of the integral equation to allegorical equations difficult. A particular scheme for overcoming this difficulty is the usage of an auxiliary equation that reduces the order of the original equation. Based on this scheme, an extra equation for 𝑤 is introduced in the Boussinesq equations system as follows: 𝜂_𝑡+𝜕()𝑢𝜕𝑥𝐻𝑢=0,_𝑡+𝑢𝜕𝑢𝜕𝑥+𝑔𝜕𝜂=ℎℎ𝜕𝑥1+𝛽3𝜕²𝑢_𝑡𝜕𝑥²+ℎ_𝑥𝜕𝑢_𝑡+𝜕𝑥𝜕𝑤,ℎ𝜕𝑥𝑤=𝛽𝑔ℎ3𝜕²𝜂𝜕𝑥²+ℎ_𝑥𝜕𝜂.𝜕𝑥 (6)(7)(8) The basis of this auxiliary equation is the fact that all nodal values over their corresponding elements are assumed constant. Multiplying these equations by the test functions 𝑁_𝑘 and integrating over the entire domain, the equations related to the Galerkin method are obtained. As a result of these steps, the continuity equation transforms to 𝑁_𝑘𝜂_𝑡+𝜕(𝜕𝑥𝐻𝑢)𝑑Ω^𝑒=0. (9) Pursuant to discretization of all terms in (7), the following matrix equation for an element Ω^𝑒 is found: 𝑀_𝑘𝑖𝜂̇_𝑖+𝑄_𝑘𝑖𝑗𝐻_𝑖𝑢_𝑗+𝑄_𝑘𝑖𝑗𝑢_𝑖𝐻_𝑗=0, (10) where the overdot in the first term represents the temporal differentiation with respect to time. Also, 𝑀_𝑘𝑖 and 𝑄_𝑘𝑖𝑗 are defined as 𝑀_𝑘𝑖=𝑁_𝑘𝑁_𝑖𝑑Ω^𝑒,𝑄_𝑘𝑖𝑗=𝑁_𝑘𝑁_𝑖𝜕𝑁_𝑗𝜕𝑥𝑑Ω^𝑒. (11) Using Green’s theorem and applying the integration by parts, the momentum and auxiliary equations can also be written similar to (10) as in 𝑀_𝑘𝑖+1+𝛽ℎ𝐸_𝑘𝑖−ℎ_𝑥𝑇_𝑘𝑖𝑢̇_𝑖+𝑄_𝑘𝑖𝑗𝑢_𝑖𝑢_𝑗+𝑔𝑇_𝑘𝑖𝜂_𝑖−𝑇_𝑘𝑖𝑤_𝑖𝑀=0,_𝑘𝑖𝑤_𝑖=−𝛽𝑔ℎ𝐸_𝑘𝑖−ℎ_𝑥𝑇_𝑘𝑖𝜂_𝑖, (12)(13) where 𝑇_𝑘𝑖=𝑁_𝑘𝜕𝑁_𝑖𝜕𝑥𝑑Ω^𝑒,𝐸_𝑘𝑖=13ℎ𝜕𝑁_𝑘𝜕𝑥𝜕𝑁_𝑖𝜕𝑥𝑑Ω^𝑒−𝑆𝐴_𝑘𝑖. (14)𝑆𝐴_𝑘𝑖 in (14) is the boundary integral that is defined as follows: 𝑆𝐴_𝑘𝑖=13ℎ𝑁_𝑘𝜕𝑁_𝑖𝜕𝑥⋅𝑛_𝑥𝑑Γ^𝑒. (15) Here, Γ^𝑒 identifies the boundary of the element having unit normal 𝑛. Equations (14) and (15) become necessary only when natural boundary conditions are prescribed on portions of the boundaries of the computational domain.

Finally, after discretization, the assembled global forms of (11) and (12) related to BNBE system take the form [𝑀]𝑓̇=𝐄_𝟏. (16) Also, the assembled global form of auxiliary equation (13) is obtained as follows: [𝑀]{𝑤}=𝐄_𝟐. (17) In (16) and (17), 𝑀 is the coefficient matrix, 𝑓=𝜂  or  𝑢, while 𝐄_𝟏 and 𝐄_𝟐 are vectors that are determined by the known values of 𝑢,𝜂, and 𝑤. System of linear equations in (17) can be solved explicitly at each time step, and (16) must be integrated in time by a typical high-order method.

3.2. Predictor-Corrector Scheme

If the variable 𝑓 is discretized at 𝑡=𝑛Δ𝑡, the Adams-Bashforth-Moulton (ABM) method is adopted to integrate in time domain, in the following form [17]. (1)Right-hand side of (16) is evaluated at time levels 𝑛, 𝑛−1, and 𝑛−2.(2)Equation (16) is temporally integrated by applying the predictor stage of the ABM method: [𝑀]𝑓^𝑛+1=[𝑀]𝑓^𝑛+Δ𝑡×𝐸1223₁^𝑛−16𝐸₁^𝑛−1+5𝐸₁^𝑛−2. (18)(3)Right-hand side of (16) is evaluated at time level 𝑛+1.(4)Equation (16) is temporally integrated by applying the corrector stage of the ABM method: [𝑀]𝑓^𝑛+1=[𝑀]𝑓^𝑛+Δ𝑡×9𝐄24_𝟏𝑛+1+9𝐄_𝟏^𝑛−5𝐄_𝟏^𝑛−1+𝐄_𝟏^𝑛−2. (19)These steps are repeated until convergence is reached. In the present work, the system of equations is solved using generalized minimal residual (GMRES) method [18].

3.3. Boundary Conditions

Inflow and outflow boundaries are two types of general boundary conditions which are employed here.

3.3.1. Inflow Boundary

For this type of boundary condition, initial value for water surface elevation for the entire computational domain is assumed zero, and the elevation 𝜂 on the boundary for each time step is approximated as follows: 𝜂(𝑥,𝑡)=𝑎sin(𝑘𝑥−𝜔𝑡). (20) Linear theory [19] is applied to find the wave velocity for the incident wave elevation as in 𝜔𝑢=𝑘ℎ𝜂. (21) The system of discretized equations becomes a complete set by addition of a specific variable of 𝑤 on the inflow boundary for each model. Therefore, 𝑤_𝑖 for BNBE is determined as follows: 𝑤_𝑖−𝑎𝛽𝑔ℎ−ℎ𝑘²3sin(𝑘𝑥−𝜔𝑡)+𝑘ℎ_𝑥cos(𝑘𝑥−𝜔𝑡)=0. (22)

3.3.2. Outflow Boundary

Minimization of the reflecting wave back into the domain is an important issue at the outflow boundary. Therefore, in this study, a viscous damping layer, that is, a sponge layer, is situated right before the outflow region to absorb the reflecting wave energy. For more detail of this method one can refer to the work done by Larsen and Dancy [20] and Li et al. [15].

3.4. Stability Analysis

In this segment, von Neumann stability analysis [21] has been conducted to investigate the stability of the suggested numerical model. For the stability analysis of BNBE, wave with small amplitude (linear) is considered over a constant water depth. Thus, (6), (7), and (8) are linearized as follows: 𝜂̇+ℎ𝜕𝑢̇𝜕𝑥=0,𝑢+𝑔𝜕𝜂−ℎ𝜕𝑥1+𝛽²3𝜕²𝑢̇𝜕𝑥²−𝜕𝑤ℎ𝜕𝑥=0,𝑤−𝛽𝑔²3𝜕²𝜂𝜕𝑥²=0. (23) Matrix equations for each element Ω^𝑒 can also be written as 𝑀_𝑘𝑖𝜂̇_𝑖+ℎ𝑇_𝑘𝑖𝑢_𝑗𝑀=0,_𝑘𝑖+1+𝛽ℎ𝐸_𝑘𝑖𝑢̇_𝑖+𝑔𝑇_𝑘𝑖𝜂_𝑖−𝑇_𝑘𝑖𝑤_𝑖𝑀=0,_𝑘𝑖𝑤_𝑖+𝛽𝑔ℎ𝐸_𝑘𝑖𝜂_𝑖=0. (24)(25)(26) Here, it is noted that 𝜂^𝑛_𝑖=𝜂₀𝐺^𝑛𝑒^𝐼𝜃𝑖,𝑢^𝑛_𝑖=𝑢₀𝐺^𝑛𝑒^𝐼𝜃𝑖,𝑤^𝑛_𝑖=𝑤₀𝐺^𝑛𝑒^𝐼𝜃𝑖, (27)(28)(29) where 𝐼=√−1 and 𝜃=𝑘Δ𝑥 are the imaginary unit and the phase angle, respectively. Also, (𝜂₀,𝑢₀) is the eigenvector of the problem, and 𝐺 is the amplification factor. By substituting (24) into the Adams-Bashforth predictor method (18), we can obtain the following: Δ𝑥6𝜂^𝑛+1_𝑖−1+4𝜂_𝑖^𝑛+1+𝜂^𝑛+1_𝑖+1=Δ𝑥6𝜂^𝑛_𝑖−1+4𝜂^𝑛_𝑖+𝜂^𝑛_𝑖+1−ℎΔ𝑡𝑢2423^𝑛_𝑖+1−𝑢^𝑛_𝑖−1−16𝑢^𝑛−1_𝑖+1−𝑢^𝑛−1_𝑖−1+5𝑢^𝑛−2_𝑖+1−𝑢^𝑛−2_𝑖−1. (30) Introducing (27) and (28) into (30) will result in 𝛽₁𝐺²(𝐺−1)𝜂₀+𝛽₂23𝐺²−16𝐺+5𝑢₀=0, (31) where 𝛽₁ and 𝛽₂ are given as in 𝛽₁𝛽=2cos𝜃+4,₂=ℎΔ𝑡()4Δ𝑥2sin𝜃𝐼. (32) Also, similar to the above scheme for (24), (25) becomes 𝛽₃𝐺²(𝐺−1)𝑢₀+𝛽₄23𝐺²−16𝐺+5𝑔𝜂₀+𝑤₀=0, (33) where 𝛽₃ and 𝛽₄ are given as in 𝛽₃=()−2ℎ2cos𝜃+41+𝛽²(Δ𝑥)²(),𝛽2cos𝜃−2₄=Δ𝑡()4Δ𝑥2sin𝜃𝐼. (34) Expansion of (26) at every time step would lead to Δ𝑥6𝑤^𝑛_𝑖−1+4𝑤^𝑛_𝑖+𝑤^𝑛_𝑖+1ℎ−𝛽𝑔²𝜂3Δ𝑥^𝑛𝑖−1−2𝜂^𝑛_𝑖+𝜂^𝑛_𝑖+1=0. (35) By applying (27) and (29), it is concluded that 𝑤₀=𝛽₅𝜂₀, (36) where 𝛽₅ℎ=2𝛽𝑔²(Δ𝑥)²2cos𝜃−2.2cos𝜃+4 (37) By introducing (36) to (33), (33) can be rewritten as 𝛽₃𝐺²(𝐺−1)𝑢₀+𝛽₆23𝐺²−16𝐺+5𝜂₀=0, (38) where 𝛽₆=𝛽₄𝑔−𝛽₅. (39) To have nontrivial solutions for eigenvector (𝜂₀,𝑢₀) for (31) and (38), we must have 𝐺²(𝐺−1)²−𝛽₂𝛽₆𝛽₁𝛽₃23𝐺²−16𝐺+5²=0. (40) Equation (40) can be numerically solved for |𝐺| in terms of 𝜃, by assuming Δ𝑥=0.1ℎ and Cr=(0.5,1.0,1.5) where Cr is the Courant number. That is. Cr=𝑔ℎΔ𝑡.Δ𝑥 (41) Figure 1 demonstrates the largest modulus of the amplification factor in terms of 𝜃. It is obvious that the predictor scheme is stable when the Courant number is smaller than or equal to unity.

Figure 1 

Largest modulus of the amplification factor |𝐺|max  in terms of phase angle 𝜃 for the predictor method.

The outlined technique for studying the stability of the Adams-Moulton predictor scheme in (42) can also be extended to examine the stability of the Adams-Moulton corrector method in (19). This would result in 𝐺²(𝐺−1)²−𝛼₂𝛼₄𝛼₁𝛼₃9𝐺³+19𝐺²−5𝐺+1²=0, (42) where 𝛼₁𝛼=2cos𝜃+4,₂=ℎΔ()𝛼8Δ𝑥2sin𝜃𝐼,₃=()−2ℎ2cos𝜃+41+𝛽²(Δ𝑥)²(),𝛼2cos𝜃−2₄=ℎ𝑔−2𝛽𝑔²(Δ𝑥)²2cos𝜃−22cos𝜃+4Δ𝑡()8Δ𝑥2sin𝜃𝐼. (43) The largest modulus of amplification factor in terms of phase angle 𝜃 is illustrated in Figure 2.

Figure 2 

Largest modulus of the amplification factor |𝐺|max in terms of phase angle 𝜃 for the corrector method.

Based on this figure, it is obvious that the corrector scheme is stable when the Courant number is less than or equal to 0.5.

4. Numerical Experiments

Accuracy of the suggested model is examined by different simple and complex test cases. The first test considered is the regular wave shoaling on a constant slope beach where nonlinearity effect is examined. Time histories of the free surface elevation computed by the proposed model are compared against the existing results of experimental test conducted to predict the nonlinearity properties of the model. In continuation, in a more complex case, a time series of velocity in addition to free surface evaluation are also determined and compared with the available numerical and experimental results. Another test case considered which is more related to the dispersion property of the model is the simulation of periodic wave propagation over a bar. Here, free surface elevation determined numerically is compared against the experimental profile. Finally, the phenomenon of nonlinear wave propagation over a barred beach is modeled by the current numerical schemes. As a result, the wave decomposition is examined behind or beyond a bar, and the computed dispersed waves are compared against available numerical findings.

4.1. Regular Wave Shoaling on Constant Slope Beach

A beach with a constant slope is the simplest assumption that can be considered for the morphological purpose. To demonstrate the strength of the present numerical scheme in modeling the nonlinear terms in the system of BNBE, the model is applied to simulate the propagation of regular wave over a constant slope beach where nonlinearity effect increases (under the effect of the slopes). Computational domain is shown in Figure 3.

Figure 3 

Schematic of the computational domain for regular wave shoaling on a constant slope beach.

The computational domain with length 𝐿=19m is divided into 489 elements. The slope is 1 : 25, and the starting point of the slope with respect to the inflow boundary is located at 𝑥=4.6m. Water depth at the inflow boundary is equal to 0.056m and the regular wave is absorbed at the outflow boundary where the water depth is reduced to 0.1m. The considered wave period and wavelength are 1.0s and 1.55m, respectively. The wave is allowed to propagate until the time 𝑡=30.0s, with a time step Δ𝑡=0.025s. The velocity boundary condition at the incoming boundary for the simulation of regular wave is known from the linear wave theory.

In order to analyze the accuracy of the obtained results by the current numerical scheme, time histories of the free surface level at the first and second gages, which are located at ℎ=0.24m and at ℎ=0.1m, are determined. The synchronized results of the numerical method and the available experimental data [5] are presented in Figure 4 for the sake of comparison.

(a)

(b)

Figure 4 

Comparison of time series of regular wave shoaling on a constant slope beach: ℎ at (a) 0.24m and (b) 0.1m.

In both of the presented plots, the time histories of water surface elevation, computed at the first and second gages, using the suggested model, are in favorable agreement with the experimental data. It is quite clear that the presented numerical model has good ability to model the wave nonlinearity effects in the shoaling phenomena.

4.2. Long Periodic Wave Propagation over Sloping Beaches

In this case study, the ability of the present model for determining the horizontal particle velocity near the bed and the nonlinearity effect of the free surface are investigated. Good prediction of horizontal particle velocity may be important when the cross-shore sediment transport in nearshore is to be modeled, if the Boussinesq equation is to be applied. Gilbert et al. [22], in continuation of the research done by Grilli et al. [23], studied the propagation of the long periodic waves over semiburied cylindrical objects at the bottom of a wave tank. The aim of that study was to find a numerical model for simulation of nonlinear wave forced sediment transport, where near bed flow increases near the objects in shallow water coastal area. For this situation, they offered an experimental as well as a numerical work. Their experiment involved a laboratory wave tank with a sandy bottom, while their numerical model included solution of a two-dimensional fully nonlinear potential flow (FNPF). Some experiments were run without any bottom obstacles. Here, these experimental cases were conducted to calculate the free surface and the particle velocity by one-dimensional BNBE model at definite gages for the case where there was no obstacle at the bottom. It was shown that the effect of the barrier is quite low indicating a maximum of 6% error in the water surface level [23]. Schematic of their numerical wave tank is illustrated in Figure 5. Wave is generated at the left boundary (𝑥=0) at the initial water depth of ℎ₀=1.05m. Bed slope starts with a 1 : 6 slope at 𝑥=1.219m and is followed by a 1 : 24 slope at 𝑥=3.048m.

Figure 5 

Schematic of the computational domain for 1D-BNBE computations of a periodic wave deep with height 𝐻0 and period 𝑇.

As detailed in Figure 5, the surf zone is modeled by AB following the AB (ramp) which itself starts at 𝑥_𝑎=12.8m in the AB. The seabed slope diminishes at the point 𝑥_𝑏=16.5m, prior to the wave breaking point, and then (for 𝑥>𝑥_𝑏) the water depth gradually increases to ℎ₁=0.5m and rapidly ripples under the wave action. Experimental measurements are only made when the seabed shape levels out. The fully nonlinear potential flow model did not offer the precise shape of the ripples, as this would not strongly influence the shoaling waves [22, 23].

To dissipate the energy in order to prevent wave overturning, wave energy was absorbed for 𝑥>𝑥_𝑏 in the present model and also in the FNPF model.

Similar to the experimental test, a periodic wave with period 𝑇=2𝜋/𝜔=2.5s and height 𝐻₀=0.19m is simulated by the present model. Computational domain (shown in Figure 5) is discretized by 400 linear finite elements, and each time step is equal to Δ𝑡=0.025s. Free surface elevations obtained by BNBE model are compared with the available experimental data and FNPF model at gages G1 (𝑥=12.802m) and G2 (𝑥=14.021m) as shown in Figure 6. Here, consistent trends are displayed by model at both gages. The agreement is quite good for the results at gage G1, but the computed result shows a little less overestimation at G2. On the other hand, horizontal velocity was calculated by BNBE model and compared with the available experimental data and the results of FNPF model. Horizontal velocity data was measured at 0.1m above the bottom at gages G1 and G2 without the obstacle. As stated earlier, BNBE model computes the depth-integrated horizontal velocity. Figure 7 shows these comparisons. The results display overestimation of the maximum velocity magnitude by 25% at all gages compared to FNPF model and experimental data. The reason is obvious. The velocity near the bed is smaller than the averaged velocity that is predicted by the Boussinesq models. Also, the asymmetry of the horizontal particle velocity calculated by the boussinesq models, compared to that found by the FNPF model and the experimental data at these gages, is relatively good.

(a)

(b)

Figure 6 

Comparison of the results for free surface elevations by BNBE model against the existing results of FNPF model and the experimental data at gages: (a) G1; (b) G2.

(a)

(b)

Figure 7 

Comparison of the results for horizontal velocity by BNBE model against the existing results of FNPF model and the experimental data at gages: (a) G1; (b) G2.

Comparison of the results in Figures 6 and 7 indicates lager onshore velocities under crests and smaller offshore velocities under troughs. Such velocity patterns are directly related to the observed sediment transport [24]. The good agreement of the velocity displayed can be indicative of the fact that BNBE can be used to determine the cross-shore sediment which will be used in future work.

4.3. Periodic Wave Propagation over a Bar

Before applying the suggested model for simulation of the wave propagation over the barred beach, a simple bar over a constant bed is considered here with the hope that wave properties (wave height and phase celerity) can be modeled accurately after the passing of waves over the bars. This test is more important for showing the accurate numerical modeling of the dispersion properties of the models. The experimental data of this test case was firstly published by Beji and Battjes [25]. Figure 8 illustrates the computational domain for modeling the flume experiment.

Figure 8 

Bathymetry for the bar.

The input wave is of wavelength 𝜆=3.73m, wave period 𝑇=2.02s, and height 𝐻=0.02m. wave enters into a calm water of depth ℎ=0.40m before reaching the bar. To neutralize the effect of exiting wave, a sponge layer is utilized at the right boundary. The computational domain is discretized with grid size elements of Δ𝑥=0.05m in the wave propagation direction. The wave is allowed to propagate until the time 𝑡=40.0s, and a time step Δ𝑡=0.026s is used for the time integration. Comparison of the numerical free surface elevation with the available experimental data for the gage stations at 𝑥=13.5 and 17.3m is shown in Figure 9.

(a)

(b)

Figure 9 

Comparison of time series of free surface over a bar, (a) 𝑥=13.5m; (b) 𝑥=17.3m.

Results show no significant difference between the current numerical method and that of existing experimental data, and the model has the same behavior after passing the gage located at 𝑥=17.3m.

4.4. Shoaling of Regular Waves over Barred Beaches

The assumption of a constant slope beach is often not a good approximation as sandy coasts exposed to a wave climate will often build up barred beaches [1]. Schematic of this numerical test which was first used by Grilli and Horrillo (1999) is shown in Figure 10.

Figure 10 

Schematic of the computational domain for a periodic wave modeling.

The main profile considered for the beach has the same form as Dean’s [26] equilibrium beach profile. This profile has an average slope of 1 : 50 along with a bar of 0.2 m height with a 1 : 20 seaward slope and 1 : 10 shoreward slope located toward its top. Geometrical profile of this bar is similar to that used by Beji and Battjes in their experiments [25]. Similar to their test cases, wave in the current numerical study is decomposed after propagating over the bar. After passing the bar, the wave nonlinearity increases until the wave reaches the breaking zone (i.e., AB). A sponge layer is located at this zone to absorb the outgoing wave. The length of the domain is 𝐿=40.0m and is divided into 1000 linear elements along the 𝑥 axis. The constant depth in the primary part of domain is equal to 0.6m. A regular wave with a height of 0.06m is generated at the incoming boundary. In the first test conducted here, a regular wave with period 𝑇=2.394s is simulated. Free surface elevation computed by the present numerical model for 𝑡=3.592𝑇s and 𝑡=3.911𝑇s is compared against the existing fully nonlinear potential flow model (FNPF) [27] in Figure 11.

(a)

(b)

Figure 11 

Computed surface elevation for wave period 𝑇=2.394s: (a) 𝑡=3.592𝑇s and (b) 𝑡=3.911𝑇s.

In the second test, numerical results and the existing FNPF results at 𝑡=12.5𝑇s and 𝑡=13.0𝑇s are presented in Figure 12 with smaller wave period of 𝑇=2.075s.

(a)

(b)

Figure 12 

Computed surface elevation for wave period 𝑇=2.075s: (a) 𝑡=12.5𝑇s and (b) 𝑡=13.0𝑇s.

Comparison of the results produced by the present models with that of the existing numerical results obtained by the FNPF model indicates that the nonlinearity and dispersive properties of the wave are predicted well. Figure 12 shows that waves height predicted by BNBE model is underestimated just like the results obtained in the first test case. For 𝑥>26, the wave nonlinearity increases after passing over the slope, and consequently a slight phase difference is observed between the results of the present model and that of FNPF. Strong nonlinear wave decomposition in the region beyond the bar is predicted by BNBE which is similar to the results of the FNPF model. However, as the wave reaches the top of the slope, that is, the breaking zone, the nonlinearity increases and becomes stronger, but the wave height is overpredicted by the BNBE model (as shown in Figure 12(b)). Considering the fact that BNBE model is a one-dimensional approximation, the obtained numerical results may be considered remarkable.

5. Conclusions

In the current study, a particular numerical method is presented for the Beji and Nadaoka extended Boussinesq equations to determine the local wave properties (i.e., wave height and wave velocity) in nearshore zone by one-dimensional cases. Numerical forms of these equations are revisited, and, by introducing a novel form of the auxiliary variable, equations are rewritten in lower-order forms making the application of the linear finite element method possible.

In order to show how well the nonlinearity and dispersion properties of the waves are predicated, two benchmark tests have been conducted. Regular wave shoaling on a constant slope beach and periodic wave propagation over a bar have been modeled. In these test cases, accurate modeling of nonlinearity and dispersion of regular waves are of high importance, respectively. Results of these numerical case studies show good agreements for numerical model against the available experimental findings. In the case of long periodic wave propagation over the sloped beach, the nonlinearity of the wave was noticed to increase while reaching the top of the slope. However, even though the asymmetric nature of wave shape was not predicated well by Beji and Nadaoka’s model, the particles’ velocities were well predicted, compared against the experimental results. This is an indicative of the fact that the Beji and Nadaoka Boussinesq equations can be developed for computation of the cross-shore sediment transport under the influence of long periodic wave near the shore.

Modeling the periodic waves over the barred beaches was studied in another test, and the obtained results were compared with the numerical results available by the FNPF model. Comparison showed that the Beji and Nadaoka Boussinesq equations are almost capable to model the decomposition of the nonlinear waves beyond the bars. Contrary to the two-dimensional FNPF model, this was achieved only by a one-dimensional approximation.

Based on the results presented, it can be concluded that the current numerical scheme can be further developed to simulate any morphological development of coastal profiles under the influence of long periodic waves.