FINITE ELEMENT MODELING OF NONLINEAR WAVE TRANSFORMATION USING ELLIPTIC MILD SLOPE EQUATION

Accurate modeling of nonlinear wave transformation is important for studies related to harbor and nearshore design. While sophisticated finite-element models based on the elliptic mild-slope equation are often used for wave prediction, they do not include wave-wave interactions. These interactions, in general, involve transfer of energy and wave phase coupling among spectral components and are known to be quite significant especially in shoaling regions and inside harbors. To overcome this limitation, and to provide a basis for the effective modeling of nonlinear wave transformation in complex coastal and harbor environments, the development of a finite-element model based on the second-order nonlinear mild-slope equation of Kaihatu and Kirby (1995) is considered. The model uses an iterative procedure for solving the second-order boundary value problem. To ensure effective boundary treatment, a combination of frequently-used boundary conditions and the method of internal wave generation is used. Two cases involving nonlinear shoaling and harbor resonance are considered for model validation. Modeled results are compared with experimental data, and good agreement is observed in most cases. The methodology described in this paper can improve the applicability of existing finite-element based mild-slope models.


INTRODUCTION
Modeling of nonlinear wave-wave interactions in the presence of reflection, diffraction, and harbor resonance is critical for studies related to harbor design and tranquility.Many finite-element models (e.g.PHAROS, CGWAVE, MIKE21-EMS) based on the standard elliptic mild-slope equation of Berkhoff (1972Berkhoff ( , 1982) ) have been developed.Such models are well known for efficient handling of complex harbor and nearshore processes for design applications.Extensions of these models that incorporate nonlinear mechanisms such as breaking, wave-current interactions, amplitude dispersion and wave-direction-dependent boundary conditions are also available.For harbor design problems, the finite-element models, which allow the use of unstructured grids, are preferred over finite-difference models.The use of unstructured grids allows accurate delineation of domain boundaries (e.g.harbor walls, arbitrary-shaped coastlines, structural boundaries) which ensures proper imposition of boundary conditions.This is important for appropriate handling of wave reflection and diffraction effects in the vicinity of domain boundaries.In addition, for such models, the mesh density can easily be refined (in a desired region) for improved accuracy.Owing to these advantages, such models have been successfully applied to simulate wave transformation in complex harbor domains (e.g.Los Angeles/Long Beach harbor, Panchang and Demirbilek 2001;Douglas Harbor, Li et al. 2005;Venice Harbor, Kostense et al. 1988) and around coastal structures present in nearshore regions.
Yet, it must be noted that the governing equation for the above models is based on linear wave theory, and wave-wave interactions are completely ignored.However, these nonlinear interactions, which, in general, involve transfer of energy among spectral components, are known to be quite significant especially in shoaling regions.In harbors, the interactions between short period waves may generate long period waves, which can cause harbor resonance (Rogers and Mei 1978;Woo and Liu 2004) with large natural periods.
To incorporate nonlinear wave-wave interactions in mild-slope models, extensions of the standard mildslope equation were proposed by Kaihatu and Kirby (1992) and Tang and Ouelett (1994).The extension includes quadratic nonlinear terms related to interactions between triads of waves.In a recent study, Sharma et al. (2014) considered the second-order extension of Kaihatu and Kirby (1992) for the development of a finite-difference model.Sharma et al. (2014) also proposed an alternative equation (with minor adjustments) which was shown to have improved convergence properties.Their finite-difference model uses an effective method for boundary treatment and provides stable solutions using an efficient iteration procedure.They investigated the performance of their model by simulating nonlinear wave transformation for a variety of cases involving reflection, diffraction, breaking and triad interactions.Their model showed, in most cases, good agreement with measurements and the results from other relevant models.However, for complex domains involving arbitrary shaped boundaries, Sharma et al. (2014) suggested developing a model which uses unstructured grids.Here we discuss the development of a finite-element model which solves the same boundary value problem considered by Sharma et al. (2014).The methodology discussed in this paper can improve the applicability of the existing finite-element models mentioned above.
In this paper, we present the finite-element formulation for this complex problem and its implementation.Two cases are then used for model validation.

THE BOUNDARY VALUE PROBLEM
The second-order extension of the elliptic mild-slope equation, which includes quadratic nonlinear terms related to triad wave interactions, was introduced by Kaihatu and Kirby (1992) in the form where g = gravitational acceleration.Eq. ( 1) represents a set of N coupled elliptic equations where the subscript n(= 1, 2, ..., N) represents the nth frequency component.The right hand side in Eq. ( 1) contains quadratic nonlinear terms (hereafter referred to as QNLs).C n and C gn are phase and group velocity, respectively, for the nth frequency component; the complex time-dependent wave potential function φn is related to the resultant velocity potential Φ using where φn (x, y, t, t) = φn (x, y, t) z is the vertical coordinate starting at the SWL; and the wave angular frequency ω n and the wave number k n for the nth frequency component are related by the linear dispersion relation.The function f n (z) = cosh(k n (h + z))/cosh(k n h), where h denotes water depth.
In a recent study, Sharma et al. (2014) proposed an alternative form of Eq. ( 1): The above equation, in comparison to Eq. ( 1), contains one additional (time-dependent) term on the left hand side.Sharma et al. (2014) showed that Eq. ( 4) exhibits improved convergence properties when solved using the iterative procedure commonly used for handling nonlinearities (e.g.breaking, nonlinear dispersion) in elliptic models.They further extended Eq. 4) to include breaking, friction and source function associated with the internal wave generation: We refer the reader to Sharma et al. (2014) for details regarding the breaking factor α n , friction factor D n and source function {S } n .The development of a more sophisticated finite-element model based on the same boundary value problem discussed above is presented in the next section.

FINITE ELEMENT MODEL
To develop a finite element model over a two-dimensional computational region Ω enclosed by boundary Γ, we first redefine the boundary value problem (discussed above) in a concise form: the governing equation and the generalized boundary condition Note that the function f n on the right hand side in Eq. ( 6) has contributions due to nonlinear interaction effects in {QNL} n and source function {S } n ; b includes frictional and breaking dissipation effects.Eq. ( 6) is the generalized form of the standard linear boundary conditions on the total domain boundary Γ which is comprised of closed boundary (Γ c ) and open boundary (Γ o ) segments.
The development of a Galerkin finite-element model of time-dependent problems, in general, involves two main steps: (1) semidiscretization, and (2) time approximation.Details regarding such problems can be found in Reddy (1993).In the semidiscretization step, a weak formulation of the boundary value problem is first established (not shown here for brevity), and the spatial approximation of the dependent variable over finite elements is defined in the form where N e j for j = 1, 2, ..m is shape function.In Eq. ( 8), it is assumed that the time dependence of dependent variable φn is separable from its spatial variation.The substitution of Eq. ( 8) in the weak form results in a set of ordinary differential equations (ODEs) in time: where the element matrices [M e ], [K e ], { f e } and {Q e } are defined using f n N i dxdy. and q n N i dxdy.
Next, the element matrices are assembled over all elements to get the following matrix form: where [M], [K], F n and Q are assembled matrices.
In the second "time-approximation" step, the time discretization of the ODEs in Eq. ( 12) is obtained by using a finite-difference scheme, generating a set of algebraic equations which is solved to obtain the solution at all time steps.Here, for time-approximation, we apply backward-difference method which is an unconditionally stable scheme for such problems.It results in a system of equations (at each step) in the form Finally, the above system is solved using the standard Conjugate Gradient method to obtain solutions at all time steps.A simulation continues until the steady-state is achieved.

MODEL VALIDATION
For preliminary testing, the model was verified (with QNLs = 0) for linear problems for which analytical/other solutions are known.Good agreement is found between modeled results and analytical solutions, and the numerical solution scheme discussed earlier is found to perfrom satisfactorily for these linear problems.(Results for these problems are not shown for brevity.)Two additional tests are considered from the literature for data-model comparisons, and also to verify the performance of the present model in the presence of nonlinear interactions.The tests involve: (1) two-dimensional nonlinear shoaling over a topographic lens; and (2) nonlinear harbor resonant interactions in a rectangular harbor.

Wave Propagation over Semicircular Shoal
Many researchers (e.g.Liu et al. 1985, Madsen and Sørensen 1992, Kaihatu and Kirby 1992) have considered experimental data of Whalin (1971) to validate their model's capability to simulate nonlinear wave shoaling.The set-up for Whalin's experiments is shown in Fig. 2. The wave tank had a length of 25.60 m and a width of 6.10 m, and the topography consists of a semicircular shoal in the central portion of the tank (Fig. 2).Three sets of experiments were conducted by Whalin for incident periods of T = 1, 2 and 3 s.The present model is applied to the T = 1 s case with incident amplitude of 1.95 cm.Two harmonics are found to be sufficient for this case.For numerical simulation, the internal wave generation method is used with the generation line located near X = 0.0 m, and the sponge layers are placed to absorb waves near the offshore and the coastal boundaries.The higher-harmonics are assumed to have zero incident wave amplitude at the generation line.
Modeled results along the centerline Y = 3.078 m are compared with experimental data in Fig. 3.It can be seen in the plot that the wave-field is roughly linear near the wave-maker boundary; however, due to wave focusing over the shoaling region, the nonlinear interactions become significant and the second harmonic starts to grow.Note that the oscillating wavefield pattern near the shoaling region for the T = 1 s case should not be confused with the spurious oscillations caused by the approximate boundary conditions.According to Madsen and Sørensen (1992), these oscillations are due to the interactions between bound and free harmonics inside the shoaling region.
For an unstructured grid with 7100 triangular elements, and a "time" step of ∆t/T min =0.50, a total of 130 steps are found sufficient to achieve steady-state on a computer with 2 GB RAM and a 3 GHz processor.The computational time is less than 5 mins.

Nonlinear Harbor Interactions
Harbor resonance in the presence of wave-wave interactions may have a significant impact on harbor design applications.Rogers and Mei (1978) investigated the effects of nonlinear energy transfer on resonant excitations inside idealized rectangular harbors with narrow harbor region (see Fig. 4).Three harbor lengths in their experiments corresponded to the first three resonant peaks of the fundamental frequency.In a more recent study, Woo and Liu (2004) used a Boussinesq-type model with improved dispersive and nonlinear characteristics to simulate harbor interactions.
Here we use the experimental data of Rogers and Mei (1978) to examine the present model's capabilities to simulate nonlinear harbor interactions in the presence of strong reflection caused by harbor walls and coastlines.The computational domain for a 1:10 prototype harbor model (similar to Woo and Liu (2004)) of the experiments with the longest bay configuration is shown in Fig. 4.Only half the domain is considered for computations, for reasons of symmetry.The prototype model has a uniform water depth of 1.53 m all over the domain.All harbor walls, including the coastline, are assumed to be fully reflecting.Input conditions for the first three wave harmonics are similar to those used by Woo and Liu (2004).The primary and secondary harmonics have incident amplitudes of 0.035 m and 0.01 m respectively.Wwave periods for primary and secondary harmonics are 4.90 s and 2.45 s respectively.The third harmonic is assumed to have zero incident amplitue.The wave generation line is located at x = −18 m (shown in Fig. ??), and the scattered waves are absorbed at the upwave boundary using a two wavelength-wide sponge layer placed near the generation line.
It is observed that the wavefield inside the bay is mostly one-dimensional (not shown) which is consistent with the experiments.The modeled results along the bay axis (y = 0.0) are compared with the experimental data of Rogers and Mei (1978) in Figs. 5 for all wave harmonics.The existence of nodes and antinodes inside the bay is due to the standing wave formation.Although the model behaves satisfactorily for both cases, some mismatch between the model results and the experimental data is observed in the vicinity of bay entrance.A similar discrepancy was observed by Woo and Liu Woo and Liu (2004) with their Boussinesq-type model.As shown in Figs. 5, the modeled harmonics are in good agreement with the experiments for x ≥ 6 m.
Using an unstructured grid with 12,310 triangular elements, the steady-state solution is obtained in 300 steps.High-resolution mesh is used near the bay entrance.

CONCLUSIONS AND RECOMMENDATIONS
This study describes the development of a finite-element model based on the second-order extension of the standard elliptic mild-slope equation.The proposed model is capable of simulating nonlinear wave transformation in complex harbor and coastal domains.The model, using unstructured grids, allows accurate imposition of boundary conditions and handling of reflection/diffraction.For effective boundary treatment, the commonly used linear boundary conditions are supplemented with internal wave generation method and dissipative sponge layers.Model performance is verified for two different cases involving nonlinear shoaling and harbor resonance, and satisfactory model performance suggests that the proposed approach will enhance the applicability of the existing finite-element models based on the elliptic mildslope equation.In future, the model will be applied to real-world applications involving arbitrary shaped boundaries and complex bathymetric features.Rogers and Mei Rogers and Mei (1978).(Top) First harmonic; (middle) second harmonic; (bottom) third harmonic.

Figure 4 :
Figure 4: Computational domain for nonlinear harbor problem. Figure not drawn to scale.

Figure 5 :
Figure 5: Comparison of wave amplitudes along the longitudinal axis of the bay.Modeled harmonics (-) and experimental data (o) of Rogers and MeiRogers and Mei (1978).(Top) First harmonic; (middle) second harmonic; (bottom) third harmonic.