Next Article in Journal
A Quantitative Flood-Related Building Damage Evaluation Method Using Airborne LiDAR Data and 2-D Hydraulic Model
Next Article in Special Issue
Wake Effect Assessment in Long- and Short-Crested Seas of Heaving-Point Absorber and Oscillating Wave Surge WEC Arrays
Previous Article in Journal
Calibrating a Hydrological Model by Stratifying Frozen Ground Types and Seasons in a Cold Alpine Basin
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Internal Wave Generation in a Non-Hydrostatic Wave Model

1
Department of Civil Engineering, Ghent University, Technologiepark 60, 9052 Ghent, Belgium
2
Flanders Hydraulics Research, Berchemlei 115, 2140 Antwerp, Belgium
3
Department of Hydraulic Engineering, Delft University of Technology, Stevinweg 1, 2628 CN Delft, The Netherlands
*
Author to whom correspondence should be addressed.
Water 2019, 11(5), 986; https://doi.org/10.3390/w11050986
Submission received: 18 April 2019 / Revised: 6 May 2019 / Accepted: 7 May 2019 / Published: 10 May 2019

Abstract

:
In this work, internal wave generation techniques are developed in an open source non-hydrostatic wave model (Simulating WAves till SHore, SWASH) for accurate generation of regular and irregular long-crested waves. Two different internal wave generation techniques are examined: a source term addition method where additional surface elevation is added to the calculated surface elevation in a specific location in the domain and a spatially distributed source function where a spatially distributed mass is added in the continuity equation. These internal wave generation techniques in combination with numerical wave absorbing sponge layers are proposed as an alternative to the weakly reflective wave generation boundary to avoid re-reflections in case of dispersive and directional waves. The implemented techniques are validated against analytical solutions and experimental data including water surface elevations, orbital velocities, frequency spectra and wave heights. The numerical results show a very good agreement with the analytical solution and the experimental data indicating that SWASH with the addition of the proposed internal wave generation technique can be used to study coastal areas and wave energy converter (WEC) farms even under highly dispersive and directional waves without any spurious reflection from the wave generator.

1. Introduction

Numerical wave propagation models are commonly used as engineering tools for the study of wave transformation in coastal areas. The number of numerical models based on the Navier-Stokes equation has recently increased remarkably, since they offer detailed and accurate predictions of the wave field but at a very high computational cost. As a result, numerical models solving approximated equations usually averaged over the vertical—like Boussinesq and nonlinear shallow water equations—are an essential tool, especially when large domains and sea states of several hours are considered.
Boussinesq-type models [1,2,3] account for both nonlinearity and frequency-dispersion of the waves by using high-order derivative terms in the equations. As an alternative, models based on the non-hydrostatic approach [4,5,6] can resolve the vertical flow structure and can improve their frequency dispersion by using additional layers rather than increasing the order of derivatives as in the case of Boussinesq-type models. A representative model of this latter category is SWASH [7], a phase-resolving wave propagation model based on the nonlinear shallow water equations with added non-hydrostatic effects. In the field of wave transformation in coastal areas, SWASH has already reached a fairly mature stage since it can incorporate nonlinear shallow-water effects, such as the generation of bound sub- and super-harmonics and near-resonant triad interactions [8,9,10,11].
In order to simulate waves in the nearshore zone correctly, the generation and absorption of waves at the boundary of models need to be modelled accurately. In the SWASH model, incident waves are generated by prescribing their horizontal velocity component normal to the boundary of the computational domain over the vertical direction. Additionally, to absorb and to prevent re-reflections in front of the numerical wave generator, a weakly reflective wave generation boundary condition [12] is applied in which the total velocity signal is a superposition of the incident velocity signal and a velocity signal of the reflected waves [13]. Verbrugghe et al. [14] applied a similar method to create open boundaries within a Smoothed Particle Hydrodynamics (SPH) solver. However, this method is based on the assumption that the reflected waves are small amplitude shallow water waves propagating perpendicular to the boundary of the computational domain and hence this method is weakly reflective for directional, dispersive waves. Furthermore, Wei and Kirby [15] found that this type of radiation condition at the wave generator boundary can lead to numerical errors when long time simulations are performed. Recently, a generating absorbing boundary condition (GABC) has been developed which is an enhanced type of a weakly reflective wave generation boundary condition and can partially absorb dispersive and directional waves [16,17]. However, in this method the level of re-reflection strongly depends on the initial approximations since the characteristics of the reflected waves (i.e., wave angle, wave celerity) inside the numerical domain cannot be estimated a priori. In general, it is not possible to find practical boundary conditions that do the above task perfectly.
On the other hand, models utilizing a sponge layer are very effective in absorbing reflected waves. However, this implies that the waves have to be generated inside the computational domain instead of on the boundary. This internal wave generation technique in combination with numerical wave absorbing sponge layers was firstly proposed by Larsen and Dancy [18] for Peregrine’s classical Boussinesq equations [19]. Later, Lee and Suh [20] and Lee et al. [21] achieved wave generation for the mild slope equations of Reference [22] and the extended Boussinesq equations of Nwogu [23], respectively, by applying the source term addition method. Lee et al. [21] have shown empirically that the velocity of disturbances caused by the incident wave can be properly obtained from the viewpoint of energy transport. Further, Schäffer and Sørensen [24] theoretically derived the energy velocity by adding the delta source function to the mass conservation type equation and integrated asymptotically the resulting equation at the generation point. However, Wei et al. [25] found that the source term addition method in a single source line may cause high frequency noise in case of non-staggered computational grid. To deal with this problem, Wei et al. [25] derived a spatially distributed (Gaussian shape) source function for internal wave generation, where a mass source is added in the continuity equation. Later, Choi and Yoon [26] and Ha et al. [27] used this technique in a Reynolds-averaged Navier-Stokes (RANS) equations model and a Navier-Stokes equations (NSE) model, respectively. However, in both cases they used directly the formula of [25], which was derived from the Boussinesq equations under the shallow water assumption. Thus, their model was not capable to accurately generating deep water waves and consequently high-frequency components in case of irregular waves.
In the present paper, a source term addition method and a spatially distributed source function for internal wave generation are implemented in the non-hydrostatic model SWASH, in order accurately generate regular and irregular long-crested waves. In addition, the energy velocity will be derived for the governing equations of SWASH in case multiple layers are implemented. To the present authors’ knowledge, these internal wave generation techniques which are commonly used in Boussinesq models [3] and mild-slope wave models [28], have not been derived and used in a non-hydrostatic wave model before, due to the complexity of the governing equations. So, the main objective of the present work concerns the implementation of the internal wave generation in SWASH to accurately generate even highly dispersive and directional waves, while at the same time re-reflections at the position of the wave generator are vanishing. Moreover, SWASH has already been used in the field of marine energy by simulating the wave-induced response of a submerged wave-energy converter [29]. This kind of application requires a homogeneous wave field in the whole numerical domain. Thus, the implementation of internal wave generation is important as it introduces noteworthy improvements in the model, which can then make full use of its benefits for the study of WEC farms.
The paper is structured as follows. The numerical model SWASH is described in Section 2 where the governing equations are presented. The implemented internal wave generation techniques are derived in Section 3. Section 4 provides a detailed overview of the model results for regular and irregular long-crested waves. In addition, validation results are presented where the accuracy of the model is compared with experimental data. The last sections provide conclusions and a summary discussion of the present study.

2. SWASH Model

2.1. Mathematical Formulation

SWASH is an open source non-hydrostatic wave model [7]. The governing equations of the model are based on the nonlinear shallow water equations with added non-hydrostatic effects, which are derived from the incompressible Navier-Stokes equations. The numerical domain is bounded vertically by the free-surface z = η ( x , t ) and the bottom z = d ( x ) , where t is time, d is the still water depth and x and z are the Cartesian coordinates, where z is directed positive upwards. The equations for the 2D vertical domain are
u x + w z = 0
u t + uu x + wu z = 1 ρ 0 ( p h + p nh ) x + τ xz z + τ xx x
w t + uw x + ww z = 1 ρ 0 p nh z + τ zz z + τ zx x
where u ( x , z , t ) and w ( x , z , t ) are the horizontal and vertical velocities, p h and p nh are the hydrostatic and non-hydrostatic pressures, ρ 0 is density and τ xx , τ xz , τ zx , τ zz are the turbulent stresses. The kinematic boundary conditions are
w = η t + u η x ( z = η )
w = u d x ( z = d )
The dynamic boundary conditions at the surface are a constant pressure and no surface stresses. The free surface elevation η is determined by considering the continuity for the entire water column
η t + x d η udz = 0
At the bottom boundary, a bottom stress is included based on a quadratic friction law
τ b = c f U | U | h
where h = η + d is the total water depth, U is the depth averaged velocity and c f is a dimensionless friction coefficient.
At the outlet of the domain a sponge layer can be employed to minimize reflections. The sponge layer formula as described in Reference [30] is used, where the free surface elevation η and the velocity component u are relaxed at every time step.
In SWASH the full process of wave breaking is not simulated but instead a breaking wave is considered analogous with a hydraulic bore. However, a high vertical resolution is needed to resolve accurately this process. Hence, Smit et al. [31] proposed a breaking formulation which can reproduce wave breaking with a coarse vertical resolution. A hydrostatic pressure distribution is assumed when η / t > a gh , where g is the gravitational acceleration and a is the maximum wave steepness before breaking.
The numerical implementation is based on an explicit, second order finite difference method for staggered grids, where the mass and momentum are strictly conserved at a discrete level. In the horizontal direction rectilinear or orthogonal curvilinear grid can be applied, while in the vertical direction the computational domain is divided into a fixed number of layers. A more detailed overview of the numerical methods and equations are given in Reference [7].

2.2. Wave Generation in SWASH

In the SWASH model, incident waves are generated by prescribing their horizontal velocity component normal to the boundary of the computational domain over the vertical direction. To prevent re-reflections in front of the numerical wave generator, a weakly reflective wave generation boundary condition is adopted [12]. For the case of one layer the imposed depth-averaged horizontal velocity component is given by
u ( t ) = ω kd η t + g d ( η t η i )
where ω is the angular frequency, k is the wave number, d is the still water depth, g is the gravitational acceleration, η t is the target surface elevation and η i is the instantaneous surface elevation. In case of multiple layers, a hyperbolic cosine distribution is assumed for the horizontal velocity component. As it can be noticed from the second term of Equation (8), this method is valid only in case that the reflected waves are shallow water waves propagating perpendicular to the boundary of the computational domain, since a shallow water phase velocity, C = gd is used. Hence, it is assumed that the high frequency energy gets dissipated inside the numerical domain before arriving to the wave generation boundary [31]. In addition, a linear superposition of the incident and reflected wave is assumed and thus can only be applied to small amplitude waves.

3. Internal Wave Generation

In cases of directional and dispersive waves an internal wave generation technique can be applied to avoid re-reflections due to the weakly reflective wave generation boundary. Two different internal wave generation techniques are proposed here and are implemented in the SWASH model. The first one is a source term addition method proposed by Lee et al. [21], while the second one is a spatially distributed source function proposed by Wei et al. [25].

3.1. Energy Velocity

At first, a new energy velocity for the system of SWASH equations is mathematically derived following the methodology of Reference [24] who derived the energy velocity for Nwogu’s Boussinesq equations. In this paper, the derivation for the case of two equidistant vertical layers is presented, while the extension to more layers is straightforward.
Bai and Cheung [32] showed that the governing Equations (1)–(3) can be converted into a Boussinesq form by expressing the non-hydrostatic pressure and vertical velocity component in terms of surface elevation and horizontal velocity component. Adding a point source to the continuity equation, the linearized governing equations on a constant depth for the case of two equidistant vertical layers read
η t + d u x = Λ ( t ) δ ( x )
u t + g η x 5 16 d 2 3 u x 2 t 1 8 d 2 3 u ^ x 2 t = 0
u ^ t 1 8 d 2 3 u x 2 t 1 16 d 2 3 u ^ x 2 t = 0
where u is the depth integrated horizontal velocity component, u ^ is the inter-layer velocity variations, δ is the Dirac delta function and Λ is the source function applied to a single point. From Equation (9) it can be noticed that the u / x term must be the one balancing the Dirac delta function and thus u and u ^ have a discontinuity.
Taking the time derivative of Equation (10) and eliminating η from Equation (9) gives
2 u t 2 gd 2 u x 2 5 16 d 2 4 u x 2 t 2 1 8 d 2 4 u ^ x 2 t 2 = g   Λ δ x
Eliminating the second and third term of Equation (12) using Equation (11) to get
2 u t 2 8 g d u ^ + 1 2 gd 2 u ^ x 2 5 2 2 u ^ t 2 + 1 32 d 2 4 u ^ x 2 t 2 = g   Λ δ x
Taking the x derivative two times and eliminating the first term using Equation (11) gives
8 1 d 2 2 u ^ t 2 8 g d 2 u ^ x 2 + 1 2 gd 4 u ^ x 4 3 4 u ^ x 2 t 2 + 1 32 d 2 6 u ^ x 4 t 2 = g   Λ 3 δ x 3
Integrating 4 times Equation (14) from x = ε to x = + ε , with the limit ε 0 and requiring the integrals of u ^ to be continuous at x = 0 , we get
1 2 gd ( u ^ + u ^ ) + 1 32 d 2 ( 2 u ^ + t 2 2 u ^ t 2 ) = g   Λ
In Equation (15), u ^ is an odd function and thus the left and right contributions are identical.
In SWASH an approximation of the exact linear dispersion relation, which depends on the number of vertical layers, is used (see Appendix A). For the case of two equidistant layers it is given by
ω 2 = gk 2 d 1 + 1 16 ( kd ) 2 1 + 3 8 ( kd ) 2 + 1 256 ( kd ) 4
In addition, at the wave generation boundary, apart from the progressive waves, evanescent modes are included as well. These evanescent modes are a general characteristic of the Equations (9)–(11). Hence, Equation (16) can be rewritten as
ω 2 d g = ω e 2 d g ( kd ) 2 1 + 1 16 ( kd ) 2 1 + 3 8 ( kd ) 2 + 1 256 ( kd ) 4 = ( k e d ) 2 1 + 1 16 ( k e d ) 2 1 + 3 8 ( k e d ) 2 + 1 256 ( k e d ) 4
where k and k e are the wave numbers for progressive waves and evanescent modes, respectively. Equation (17) yields two solutions
( k e d ) 2 ( kd ) 2 , ( k e d ) 2 1 + 1 16 ( kd ) 2 1 16 + 5 256 ( kd ) 2
Furthermore, away from the source point, x = 0, the solutions of Equations (9)–(11) can be written as
η = η p 0 exp [ i ( ω t kx ) ] + η e 0 exp [ i ( ω t k e x ) ]
u = u p 0 exp [ i ( ω t kx ) ] + u e 0 exp [ i ( ω t k e x ) ]
u ^ = u ^ p 0 exp [ i ( ω t kx ) ] + u ^ e 0 exp [ i ( ω t k e x ) ]
where the subscripts p and e stand for progressive and evanescent modes, respectively. By using Equations (10) and (11), the inter-layer velocity amplitude u ^ p 0 can be expressed in terms of surface elevation amplitude η p 0
u ^ p 0 = 1 8 gk ω ( kd ) 2 1 + 3 8 ( kd ) 2 + 1 256 ( kd ) 4 η p 0
In addition, the double integral of u ^ with respect to x gives an odd function and it must vanish at the source point, x = 0 and thus from Equation (21) we can write
u ^ e 0 = k e 2 k 2 u ^ p 0
Finally, Equation (21) at x 0 + is imported in Equation (15) and eventually u ^ e 0 and u ^ p 0 are going to be replaced using Equations (22) and (23) respectively. As a result, Equation (15) becomes
Λ = 2 C e η p
where
η p = η p 0 exp [ i ω t ]
C e = 1 8 gk ω d ( ( kd ) 2 1 + 3 8 ( kd ) 2 + 1 256 ( kd ) 4 ) ( 1 k e 2 k 2 ) ( 1 2 1 32 ω 2 d g )
By eliminating ω and k e using Equations (16) and (18), respectively, it can be noticed that the energy velocity C e is equal to the group velocity C g as it can be calculated from the approximated dispersion relation by taking the derivative of ω with respect to k
C e = 64 dg ( 256 + 32 ( kd ) 2 + 5 ( kd ) 4 ) ( 16 + ( kd ) 2 ) ( 256 + 96 ( kd ) 2 + ( kd ) 4 ) 3 = ω k = C g
Similarly, the energy velocity (group velocity) C e can be derived for any number of vertical layers. In Figure 1, the normalised energy velocities C e / C g   Airy for one, two and three vertical layers are plotted as a function of dimensionless depth kd . In addition, the range of dimensionless depth kd as a function of number of vertical layers where the relative error in the normalised energy velocities C e / C g   Airy stays below 3% is given in Table 1. It can be observed that by increasing the number of layers, a better fit is achieved with the exact linear solution C e / C g   Airy = 1, fact that makes the model applicable for higher kd values.

3.2. Source Term Addition Method

In the source term addition method proposed by Reference [21], additional surface elevation η * is added with the desired energy to the calculated surface elevation η at the wave generation line for each time step and is given by
η * = 2 η I C e Δ t Δ x cos θ
where Δ x is the grid size in the x-axis, Δ t is the time step, θ is the angle of the incident wave ray from x-axis, η I is the water surface elevation of incident waves and C e is the energy velocity.

3.3. Spatially Distributed Source Function

In the spatially distributed source function proposed by Reference [25], a spatially distributed mass is added in the continuity equation. Hence, the spatially distributed source function is applied on an area in contrast to the source term addition method which is applied on a line as it is demonstrated in Figure 2. For a single wave component, the source function can be defined as follows
f ( x , y , t ) = g ( x ) Dcos ( ω t kysin θ )
where D is the source function amplitude and g ( x ) is the shape of the source function which can be arbitrarily chosen. Here, a smooth Gaussian shape has been applied. Wei et al. [25] derived the source function amplitude D for the extended Boussinesq equations of [23] and is given by
D = 2 η 0 ( ω 2 α 1 gk 4 d 3 ) cos θ ω kI ( 1 α ( kd ) 2 )
where α = 0.390 , α 1 =   α + 1 / 3 and I is the integral defined by
I = π β exp ( ( kcos θ ) 2 4 β )
where β = 20 / W 2 and W is the width of the source area.
Equation (30) is derived for the extended Boussinesq equations of Nwogu and thus cannot be applied directly to the SWASH equations for the case of multiple layers. However, Equation (30) can be written in terms of the energy velocity as
D = 2 C e η 0 cos θ I
Equation (32) can be implemented in SWASH by using the energy velocity that corresponds to the number of vertical layers used in the model. The energy velocity for the case of two vertical layers has already been derived in Section 3.1. In this way, we overcome the limitation that [26,27] observed, where their models, using Equation (30), were not applicable to deep water conditions.

4. Model Tests

The internal wave generation techniques, implemented in SWASH as described in Section 3, have been used to generate regular and irregular long-crested waves. The obtained results are validated against analytical solutions and experimental data including water surface elevations, orbital velocities, frequency spectra and wave heights.

4.1. Regular Waves

Firstly, the applicability and accuracy of the internal wave generation techniques to generate linear regular waves in transitional and deep water conditions are examined. The numerical basin is similar to the one in Figure 2 with sponge layers at the left and right boundaries with a width of 3 L   (where L is the wave length) and a flat bottom. The internal wave generator is positioned at a distance of 3 L from the sponge layer. Two different wave conditions are examined: one in transitional water ( kd = π/6) and one in deep water ( kd = π, highly dispersive wave). For dimensionless depth kd = π/6 one vertical layer is applied while for dimensionless depth kd = π two equidistant vertical layers are applied in order to keep the relative error below 3% (Figure 1, Table 1). For both cases, the grid cell size is chosen so that Δ x = L / 50 , while the time step is equal to Δ t T / 350 , where T is the wave period.
In Figure 3, the normalised water surface elevation η / η 0 calculated using the source term addition method (blue dashed lines) and the spatially distributed source function (red dashed lines) is presented. The computed results at a distance of 1 L from the internal generator and at a time period of t = 10 T to t = 18 T are compared to the results of the analytical solution (solid black lines) for the same conditions. The agreement is very good while the behavior of the two internal wave generation techniques is similar, with a small deviation for the case of the highly dispersive wave (Figure 3c,d). For the latter case, the spatially distributed source function provides a slightly better fit with the analytical solution.
In order to check the capability of the model for handling reflected dispersive waves a computational domain with a sponge layer at the left boundary and a fully reflective wall at the right boundary is used. The internal wave generator is positioned at the middle of the computational domain which has a total length of 12L. The wave height is H = 0.02 m, the wave period is T = 3 s and the water depth is d = 10 m. These wave conditions give a dimensionless depth of kd = 4.5 and thus two equidistant vertical layers are applied. The model is applied with a grid cell size of Δ x = 0.3 m and a time step of Δ t = 0.0125 s.
Figure 4 shows a snapshot of normalised water surface elevation η / η 0 at t = 50 T generated using the source term addition method (blue solid line) and the spatially distributed source function (red dashed line). The two computed profiles are identical, while their agreement with the analytical solution (black markers) is excellent.
The waves that are generated at the internal wave generator propagate towards both ends of the domain. The waves are fully reflected at the right boundary and then are fully absorbed at the sponge layer. This, together with the fact that the target wave is a linear wave and that the length of the numerical domain is an integer multiple of the wave length, results in a profile that is a standing wave with perfect nodal points. In addition, the results show that the reflected waves pass through the generation area without distortion and that the sponge layer is capable of absorbing the highly dispersive wave. It has to be mentioned that the weakly reflective wave generation (Section 2.2) cannot be applied in this case since the assumption that the reflected waves are shallow water waves propagating with a phase velocity of   C = gd is violated.
In the proposed internal wave generation techniques, the horizontal velocity component is not prescribed over the vertical direction in contrast to the weakly reflective wave generation (Section 2.2). However, as it can be observed from Figure 5 the horizontal and vertical velocity profiles are correctly calculated for both internal wave generation techniques at a distance of 1 L from the generation point. The wave conditions are the same as in the previous test, while ten equidistant vertical layers have been applied in order to achieve a good agreement with the analytical hyperbolic profile.
Here, it has to be mentioned that in case of deep water waves, when coarse vertical resolution is applied, the weakly reflective wave generation technique needs calibration to generate the target wave height, since the hyperbolic profile of the horizontal velocity component cannot be accurately described by a small number of vertical layers. On the other hand, the internal wave generation techniques that are presented in this paper do not need calibration, since the source is directly linked with the surface elevation instead of the velocity component.
In all the previous tests in this section, linear waves have been examined. In order to study the behaviour of the internal wave generation techniques under non-linear waves, waves with different wave heights are tested. The internal wave generator is positioned at the middle of the computational domain which has a total length of 32 L with sponge layers at the left and right boundaries and a flat bottom. The wave heights are H = 0.1 m, 1 m, 2 m, 3 m, while the wave period is T = 6 s and the water depth is d = 10 m. These wave conditions give a dimensionless depth of kd = 1.3 and thus two equidistant vertical layers are applied. The model is applied with a grid cell size of Δ x = 2.0 m and a time step of Δ t = 0.025 s.
Figure 6 shows snapshots of normalised water surface elevation η / η 0 at t = 50 T for the different wave heights. Profiles generated using the spatially distributed source function are only presented here since the source term addition method becomes unstable for high wave heights. From Figure 6, it can be observed that by increasing the ratio H / d , the waves are transforming from linear (Figure 6a) to non-linear where the wave crests are becoming sharper and the wave troughs flatter.

4.2. Irregular Waves

We consider two test cases of uni-directional irregular waves fitting a JONSWAP (Joint North Sea Wave Observation Project) spectrum, with a significant wave height, H s = 0.5 m, a water depth, d = 10 m, a peak enhancement factor, γ = 3.3 and two different peak wave periods, T p = 12 s and T p = 8 s. The frequency range is confined between 0.5 f p and 3 f p . The internal wave generator is positioned at the middle of the computational domain and sponge layers are placed at the left and right boundaries. For both test cases, the model is applied with a grid cell size of Δ x = 1.0 m and a time step of Δ t = 0.025 s. However, for the case of T p = 8 s two equidistant vertical layers have been applied to correctly describe the high-frequency part of the spectrum in contrast to the case of T p = 12 s where one layer is enough.
In Figure 7, a comparison is made between the target frequency spectrum ( S t ) and the simulated frequency spectra generated using the two internal wave generation techniques for the cases of T p = 12 s (Figure 7a) and T p = 8 s (Figure 7b). The surface elevations η at the electronic wave gauges, which are positioned at a distance of 3 L p (where L p is the wave length corresponding to the peak period) from the internal wave generator, are recorded from t =   25 T p to t =   300 T p with a sampling interval of 0.2 s. The recorded data are processed in segments of 2048 points per segment. A taper window and an overlap of 20% are used for smoother and statistically more significant spectral estimates. The resulting frequency spectra agree very well with S t apart from the one that corresponds to the source term addition method for the case of T p = 12 s, where high-frequency noise can be observed. In addition, it is found that the magnitude of this noise strongly depends on the grid resolution since a coarser resolution leads to less noise. Moreover, the small difference around the peak of the frequency spectrum in Figure 7b could be caused by numerical dissipation since the electronic wave gauges are positioned at a distance of 3 L p from the internal wave generator.

4.3. Oblique Waves in a Basin with Constant Depth

In this section, oblique waves are generated in a basin with constant depth by using the spatially distributed source function wave generation technique. The numerical basin is 190 m long, 90 m wide and 1 m deep. Sponge layers are placed at the left and right boundaries with a width of 50 m while periodic boundaries are applied at the top and bottom of the domain. The internal wave generator is parallel to the y-axis and is positioned at a distance of 100 m (x = 0) from the left boundary. The wave height is H = 0.01 m, the wave period is T = 4 s and the wave propagation angle is θ =   15°. The grid cell size is chosen so that Δ x   =   Δ y   =   0.15 m. In order to obtain a steady state wave field, waves are generated for a duration of 180 s with a time step Δ t =   0.0125 s.
Figure 8 and Figure 9 present comparisons between the computed normalised water surface elevation η / η 0 at t = 40 T and the corresponding analytical solution. It can be observed that the computed solution coincides with the analytical solution except for the region of the sponge layers. The excellent agreement indicates that the spatially distributed source function implemented in SWASH is able to accurately generate directional waves.

4.4. Wave Propagation over a Shoal in a Three-Dimensional Numerical Basin

Finally, a three-dimensional version of the developed model with the spatially distributed source function wave generation technique is applied to study regular waves propagating over a shoal. The experiment conducted by Reference [33] has served as a standard test case for validating phase resolving wave models.
The bathymetry of the experimental setup consist of an elliptic shoal resting on a plane sloping seabed and the entire slope is turned at an angle of 20° with respect to the x-axis. Detailed information on the geometry can be found in Reference [33]. The numerical basin is 45 m long (−20 < y < 25) and 20 m wide (−10 < x < 10). The internal wave generator is parallel to the x-axis and is positioned at a distance of 10 m (y = −10 m) from the bottom boundary. The remaining numerical domain includes two sidewalls at x = −10 m and x = 10 m and two sponge layers at y = −15 m and y = 20 m with a width of 5 m.
The wave height is H = 0.0464 m, the wave period is T = 1 s and the water depth at the position of the internal wave generator is d = 0.45 m. These wave conditions give a dimensionless depth of kd = 1.9 and thus two equidistant vertical layers are applied. The model runs for 50 s without any stability issues since the reflected waves that are reaching the offshore boundary are absorbed by the sponge layer, which is positioned behind the internal wave generator. The model is applied with a grid cell size of Δ x = Δ y = 0.05 m, while the time step is automatically adjusted during the simulation depending on the CFL (Courant–Friedrichs–Lewy) condition where a maximum CFL value of 0.5 is used.
Figure 10 shows the plan view of the normalised wave height H / H 0 (where H is the local wave height and H 0 is the wave height at the wave generation boundary) in the whole computational domain where the diffraction and refraction patterns of the waves are visible. The wave heights of the model are obtained by averaging those of the last ten wave periods of the simulation (from t = 40 s to t = 50 s).
Figure 11 shows the comparison of normalised wave heights H / H 0   between numerical model results (red lines) and experimental data (black circles) along eight measurement transects. Additionally, in order to evaluate the model, the root mean square error (RMSE) and the Skill factor for the normalised wave heights of each section are calculated as:
RMSE =   i = 1 N ( P i O i ) 2 N   Skill = 1 i = 1 N ( P i O i ) 2 i = 1 N O i 2
where O and P indicate the observed and predicted values, respectively. Very good agreement between the numerical model and the experimental model is observed (Figure 11 and Table 2).

5. Conclusions

In the present paper new internal wave generation techniques have been developed in an open source non-hydrostatic wave model, SWASH, for accurate generation of regular and irregular long-crested waves.
Initially, two different internal wave generation techniques have been developed and implemented in SWASH model. The first one is a source term addition method based on the method proposed by Lee et al. [21], while the second one is a spatially distributed source function based on the method proposed by Wei et al. [25]. These techniques need an extension of the numerical domain to accommodate the sponge layer and the source area contrary to the weakly reflective wave generation boundary. Thus, the computational cost increases 15–30%, where the lower values stand for larger computational domains. However, the main advantage of the internal wave generation is that in cases of directional and dispersive waves the reflected waves are absorbed by the sponge layer that is positioned behind the internal wave generator in contrast with the weakly reflective wave generation boundary, which is not valid for these wave conditions due to the limitations as described in Section 2.2. Hence, the internal wave generation has a big advantage when wave energy converter (WEC) farms and man-made structures (e.g., breakwaters, artificial reefs, artificial islands) are examined, since the radiated and reflected waves, respectively, cannot be estimated a priori. In the source term addition method additional surface elevation is added to the calculated surface elevation while in the spatially distributed source function a spatially distributed mass is added in the continuity equation. In both cases, the source term propagates with a velocity which is called the energy velocity which for the system of SWASH equations is mathematically derived in Section 3.1.
Then, these wave generation techniques has been used to generate regular and irregular long-crested waves. The results indicate that the developed model is capable of reproducing the target surface elevations as well as the target frequency spectrum. The overall performance of the spatially distributed source function is better than the source term addition method since the latter becomes unstable for large wave heights and may cause high frequency noise. Finally, the developed model is also used to study wave transformation over an elliptic shoal (Berkhoff shoal experiment) where very good agreement is observed between the numerical model and the experimental results. The aforementioned observations reveal that the spatially distributed source function developed here in SWASH can be successfully used to study coastal areas and wave energy converter (WEC) farms even under highly dispersive and directional waves, while at the same time re-reflections at the position of the wave generator are vanishing.

Author Contributions

P.V. developed the fundaments of the internal wave generation implementation. P.V. set up the numerical experiments. P.V. compared the numerical experiments with experimental data. V.S., T.S., M.Z. and P.T. supervised, proofread the text and helped in structuring the publication.

Funding

This research was funded by Research Foundation—Flanders (FWO), Belgium, grant number of Ph.D. fellowship 11D9618N.

Conflicts of Interest

The authors declare no conflict of interest.

Appendix A. Derivation of Linear Dispersion Relation

The approximation of the exact linear dispersion relation for the case of two vertical layers is derived here. For the derivation, the linearized equations of [32] are used
η t + d u x = 0
u t + g η x 5 16 d 2 3 u x 2 t 1 8 d 2 3 u ^ x 2 t = 0
u ^ t 1 8 d 2 3 u x 2 t 1 16 d 2 3 u ^ x 2 t = 0
The plane wave solutions of Equations (A1)–(A3) can be written as
η = η 0 exp [ i ( ω t kx ) ]
u = u 0 exp [ i ( ω t kx ) ]
u ^ = u ^ 0 exp [ i ( ω t kx ) ]
Substituting Equations (A4)–(A6) into Equations (A1)–(A3) we get
η 0 = kd ω u 0
η 0 = ω gk ( u 0 + 5 16 ( kd ) 2 u 0 + 1 8 ( kd ) 2 u ^ 0 )
u ^ 0 = 1 8 ( kd ) 2 1 + 1 16 ( kd ) 2 u 0
Substituting Equation (A7) into Equation (A8) and eliminating the inter-layer velocity amplitude u ^ 0 using Equation (A9) the linear dispersion relationship can be obtained
ω 2 = gk 2 d 1 + 1 16 ( kd ) 2 1 + 3 8 ( kd ) 2 + 1 256 ( kd ) 4
In the same way, the linear dispersion relationship for one layer and three equidistant layers is given by
ω 2 = gk 2 d 1 1 + 1 4 ( kd ) 2
ω 2 = gk 2 d 1 + 5 54 ( kd ) 2 + 1 1296 ( kd ) 4 1 + 5 12 ( kd ) 2 + 5 432 ( kd ) 4 + 1 46656 ( kd ) 6

References

  1. Lynett, P.J.; Liu, P.L.F. A two-dimensional, depth-integrated model for internal wave propagation over variable bathymetry. Wave Motion 2002, 36, 221–240. [Google Scholar] [CrossRef]
  2. Sørensen, O.R.; Schäffer, H.A.; Sørensen, L.S. Boussinesq-type modelling using an unstructured finite element technique. Coast. Eng. 2004, 50, 181–198. [Google Scholar] [CrossRef]
  3. Shi, F.; Kirby, J.T.; Harris, J.C.; Geiman, J.D.; Grilli, S.T. A high-order adaptive time-stepping TVD solver for Boussinesq modeling of breaking waves and coastal inundation. Ocean Model. 2012, 43–44, 36–51. [Google Scholar] [CrossRef]
  4. Stelling, G.; Zijlema, M. An accurate and efficient finite-difference algorithm for non-hydrostatic free-surface flow with application to wave propagation. Int. J. Numer. Methods Fluids 2003, 43, 1–23. [Google Scholar] [CrossRef]
  5. Ma, G.; Shi, F.; Kirby, J.T. Shock-capturing non-hydrostatic model for fully dispersive surface wave processes. Ocean Model. 2012, 43–44, 22–35. [Google Scholar] [CrossRef]
  6. Bai, Y.; Cheung, K.F. Depth-integrated free-surface flow with a two-layer non-hydrostatic formulation. Int. J. Numer. Methods Fluids 2012, 69, 411–429. [Google Scholar] [CrossRef]
  7. Zijlema, M.; Stelling, G.; Smit, P. SWASH: An operational public domain code for simulating wave fields and rapidly varied flows in coastal waters. Coast. Eng. 2011, 58, 992–1012. [Google Scholar] [CrossRef]
  8. Nicolae Lerma, A.; Pedreros, R.; Robinet, A.; Sénéchal, N. Simulating wave setup and runup during storm conditions on a complex barred beach. Coast. Eng. 2017, 123, 29–41. [Google Scholar] [CrossRef]
  9. Suzuki, T.; Altomare, C.; Veale, W.; Verwaest, T.; Trouw, K.; Troch, P.; Zijlema, M. Efficient and robust wave overtopping estimation for impermeable coastal structures in shallow foreshores using SWASH. Coast. Eng. 2017, 122, 108–123. [Google Scholar] [CrossRef]
  10. Rijnsdorp, D.P.; Ruessink, G.; Zijlema, M. Infragravity-wave dynamics in a barred coastal region, a numerical study. J. Geophys. Res. Ocean. 2015, 120, 4068–4089. [Google Scholar] [CrossRef] [Green Version]
  11. Suzuki, T.; Hu, Z.; Kumada, K.; Phan, L.K.; Zijlema, M. Non-hydrostatic modeling of drag, inertia and porous effects in wave propagation over dense vegetation fields. Coast. Eng. 2019, 149, 49–64. [Google Scholar] [CrossRef]
  12. Blayo, E.; Debreu, L. Revisiting open boundary conditions from the point of view of characteristic variables. Ocean Model. 2005, 9, 231–252. [Google Scholar] [CrossRef] [Green Version]
  13. Rijnsdorp, D.P.; Smit, P.B.; Zijlema, M. Non-hydrostatic modelling of infragravity waves under laboratory conditions. Coast. Eng. 2014, 85, 30–42. [Google Scholar] [CrossRef]
  14. Verbrugghe, T.; Stratigaki, V.; Altomare, C.; Domínguez, J.; Troch, P.; Kortenhaus, A. Implementation of Open Boundaries within a Two-Way Coupled SPH Model to Simulate Nonlinear Wave–Structure Interactions. Energies 2019, 12, 697. [Google Scholar] [CrossRef]
  15. Wei, G.; Kirby, J.T. Time-Dependent Numerical Code for Extended Boussinesq Equations. J. Waterw. Port Coast. Ocean Eng. 1995, 121, 251–261. [Google Scholar] [CrossRef]
  16. Düz, B.; Borsboom, M.J.A.; Veldman, A.E.P.; Wellens, P.R.; Huijsmans, R.H.M. An absorbing boundary condition for free surface water waves. Comput. Fluids 2017, 156, 562–578. [Google Scholar] [CrossRef]
  17. Wellens, R. Wave Simulation in Truncated Domains for Offshore Applications; Delft University of Technology: Delft, The Netherlands, 2012. [Google Scholar]
  18. Larsen, J.; Dancy, H. Open boundaries in short wave simulations—A new approach. Coast. Eng. 1983, 7, 285–297. [Google Scholar] [CrossRef]
  19. Peregrine, D.H. Long waves on a beach. J. Fluid Mech. 1967, 27, 815–827. [Google Scholar] [CrossRef]
  20. Lee, C.; Suh, K.D. Internal generation of waves for time-dependent mild-slope equations. Coast. Eng. 1998, 34, 35–57. [Google Scholar] [CrossRef]
  21. Lee, C.; Cho, Y.S.; Yum, K. Internal generation of waves for extended Boussinesq equations. Coast. Eng. 2001, 42, 155–162. [Google Scholar] [CrossRef]
  22. Radder, A.C.; Dingemans, M.W. Canonical equations for almost periodic, weakly nonlinear gravity waves. Wave Motion 1985, 7, 473–485. [Google Scholar] [CrossRef]
  23. Nwogu, O. Alternative Form of Boussinesq Equations for Nearshore Wave Propagation. J. Waterw. Port Coast. Ocean Eng. 1993, 119, 618–638. [Google Scholar] [CrossRef]
  24. Schäffer, H.A.; Sørensen, O.R. On the internal wave generation in Boussinesq and mild-slope equations. Coast. Eng. 2006, 53, 319–323. [Google Scholar] [CrossRef]
  25. Wei, G.; Kirby, J.T.; Sinha, A. Generation of waves in Boussinesq models using a source function method. Coast. Eng. 1999, 36, 271–299. [Google Scholar] [CrossRef] [Green Version]
  26. Choi, J.; Yoon, S.B. Numerical simulations using momentum source wave-maker applied to RANS equation model. Coast. Eng. 2009, 56, 1043–1060. [Google Scholar] [CrossRef]
  27. Ha, T.; Lin, P.; Cho, Y.S. Generation of 3D regular and irregular waves using Navier-Stokes equations model with an internal wave maker. Coast. Eng. 2013, 76, 55–67. [Google Scholar] [CrossRef]
  28. Vasarmidis, P.; Stratigaki, V.; Troch, P. Accurate and Fast Generation of Irregular Short Crested Waves by Using Periodic Boundaries in a Mild-Slope Wave Model. Energies 2019, 12, 785. [Google Scholar] [CrossRef]
  29. Rijnsdorp, D.P.; Hansen, J.E.; Lowe, R.J. Simulating the wave-induced response of a submerged wave-energy converter using a non-hydrostatic wave-flow model. Coast. Eng. 2018, 140, 189–204. [Google Scholar] [CrossRef]
  30. Mayer, S.; Garapon, A.; Sørensen, L.S. A fractional step method for unsteady free-surface flow with applications to non-linear wave dynamics. Int. J. Numer. Methods Fluids 1998, 28, 293–315. [Google Scholar] [CrossRef]
  31. Smit, P.; Zijlema, M.; Stelling, G. Depth-induced wave breaking in a non-hydrostatic, near-shore wave model. Coast. Eng. 2013, 76, 1–16. [Google Scholar] [CrossRef]
  32. Bai, Y.; Cheung, K.F. Dispersion and nonlinearity of multi-layer non-hydrostatic free-surface flow. J. Fluid Mech. 2013, 726, 226–260. [Google Scholar] [CrossRef]
  33. Berkhoff, J.C.W.; Booy, N.; Radder, A.C. Verification of Numerical Wave Propagation Models for Simple Harmonic Linear Water Waves. Coast. Eng. 1982, 6, 255–279. [Google Scholar] [CrossRef]
  34. Berkhoff, J.C.W. Refraction and diffraction of water waves; wave deformation by a shoal, comparison between computations and measurements. In Report on Mathematical Investigation, Report W 154 part VIII; Delft Hydraulics Laboratory: Delft, The Netherlands, 1982. [Google Scholar]
Figure 1. Comparison of normalised energy velocities C e / C g   Airy for different number of vertical layers.
Figure 1. Comparison of normalised energy velocities C e / C g   Airy for different number of vertical layers.
Water 11 00986 g001
Figure 2. Definition sketch of the implemented internal wave generation methods (a) source term addition method and (b) spatially distributed source function.
Figure 2. Definition sketch of the implemented internal wave generation methods (a) source term addition method and (b) spatially distributed source function.
Water 11 00986 g002
Figure 3. Comparison between computed (dashed lines) and target (solid lines) normalised water surface elevation η / η 0 for the case of linear waves with dimensionless depth (a,b) kd = π/6 (one layer) and (c,d) kd = π (two equidistant layers).
Figure 3. Comparison between computed (dashed lines) and target (solid lines) normalised water surface elevation η / η 0 for the case of linear waves with dimensionless depth (a,b) kd = π/6 (one layer) and (c,d) kd = π (two equidistant layers).
Water 11 00986 g003
Figure 4. Snapshot of normalised water surface elevation η / η 0 at t = 50 T using two internal wave generation techniques and the analytical solution (black markers) for the case of deep water linear waves ( H = 0.02 m, T = 3.0 s, d = 10 m, two equidistant layers) with a sponge layer at the left boundary and a fully reflective wall at the right boundary.
Figure 4. Snapshot of normalised water surface elevation η / η 0 at t = 50 T using two internal wave generation techniques and the analytical solution (black markers) for the case of deep water linear waves ( H = 0.02 m, T = 3.0 s, d = 10 m, two equidistant layers) with a sponge layer at the left boundary and a fully reflective wall at the right boundary.
Water 11 00986 g004
Figure 5. Comparison between computed (markers) and target (solid lines) maximum horizontal and vertical velocity profiles using two internal wave generation techniques for the case of deep water linear waves ( H = 0.02 m, T = 3.0 s, d = 10 m, ten equidistant layers).
Figure 5. Comparison between computed (markers) and target (solid lines) maximum horizontal and vertical velocity profiles using two internal wave generation techniques for the case of deep water linear waves ( H = 0.02 m, T = 3.0 s, d = 10 m, ten equidistant layers).
Water 11 00986 g005
Figure 6. Snapshot of normalised water surface elevation η / η 0 at t = 50 T using the spatially distributed source function wave generation technique for (a) H / d = 0.01 (linear), (b) H / d = 0.1, (c) H / d = 0.2 and (d) H / d = 0.3.
Figure 6. Snapshot of normalised water surface elevation η / η 0 at t = 50 T using the spatially distributed source function wave generation technique for (a) H / d = 0.01 (linear), (b) H / d = 0.1, (c) H / d = 0.2 and (d) H / d = 0.3.
Water 11 00986 g006
Figure 7. Comparison between the frequency spectra resulting from using two internal wave generation techniques and the target frequency spectrum S t for irregular waves with H s = 0.5 m, d = 10 m (a) T p = 12 s (one layer) and (b) T p = 8 s (two equidistant layers).
Figure 7. Comparison between the frequency spectra resulting from using two internal wave generation techniques and the target frequency spectrum S t for irregular waves with H s = 0.5 m, d = 10 m (a) T p = 12 s (one layer) and (b) T p = 8 s (two equidistant layers).
Water 11 00986 g007
Figure 8. Comparison between (a) computed and (b) target normalised water surface elevation η / η 0 at t = 40 T for the case of oblique linear waves ( H = 0.01 m, T = 4.0 s, d = 1.0 m and θ = 15°) with a sponge layer at the left and right boundaries and periodic boundaries at y–direction.
Figure 8. Comparison between (a) computed and (b) target normalised water surface elevation η / η 0 at t = 40 T for the case of oblique linear waves ( H = 0.01 m, T = 4.0 s, d = 1.0 m and θ = 15°) with a sponge layer at the left and right boundaries and periodic boundaries at y–direction.
Water 11 00986 g008
Figure 9. Comparison between computed (dashed line) and target (solid line) normalised water surface elevation η / η 0 at t = 40 T for the case of oblique linear waves ( H = 0.01 m, T = 4.0 s, d = 1.0 m and θ = 15°): (a) cross section at y = 45 m, (b) cross section at x = 0 m.
Figure 9. Comparison between computed (dashed line) and target (solid line) normalised water surface elevation η / η 0 at t = 40 T for the case of oblique linear waves ( H = 0.01 m, T = 4.0 s, d = 1.0 m and θ = 15°): (a) cross section at y = 45 m, (b) cross section at x = 0 m.
Water 11 00986 g009
Figure 10. Normalised wave height H / H 0 in the whole computational domain for the experiment of Berkhoff [34].
Figure 10. Normalised wave height H / H 0 in the whole computational domain for the experiment of Berkhoff [34].
Water 11 00986 g010
Figure 11. Comparison of normalised wave heights H / H 0 between numerical model results (red solid lines) and experimental data (black circles) along different measurement transects.
Figure 11. Comparison of normalised wave heights H / H 0 between numerical model results (red solid lines) and experimental data (black circles) along different measurement transects.
Water 11 00986 g011
Table 1. Range of dimensionless depth kd as a function of number of vertical layers and the corresponding relative error in the normalised energy velocities C e / C g   Airy .
Table 1. Range of dimensionless depth kd as a function of number of vertical layers and the corresponding relative error in the normalised energy velocities C e / C g   Airy .
Number of Layerskd (-)Error (%)
1≤0.55≤3
2≤6.00≤3
3≤12.50≤3
Table 2. Root mean square error (RMSE) and Skill factor of the normalised wave heights for each section.
Table 2. Root mean square error (RMSE) and Skill factor of the normalised wave heights for each section.
Section12345678
RMSE0.0560.1160.1020.0740.1030.0820.1570.108
Skill0.9440.8840.9010.9290.9020.8890.9010.865

Share and Cite

MDPI and ACS Style

Vasarmidis, P.; Stratigaki, V.; Suzuki, T.; Zijlema, M.; Troch, P. Internal Wave Generation in a Non-Hydrostatic Wave Model. Water 2019, 11, 986. https://doi.org/10.3390/w11050986

AMA Style

Vasarmidis P, Stratigaki V, Suzuki T, Zijlema M, Troch P. Internal Wave Generation in a Non-Hydrostatic Wave Model. Water. 2019; 11(5):986. https://doi.org/10.3390/w11050986

Chicago/Turabian Style

Vasarmidis, Panagiotis, Vasiliki Stratigaki, Tomohiro Suzuki, Marcel Zijlema, and Peter Troch. 2019. "Internal Wave Generation in a Non-Hydrostatic Wave Model" Water 11, no. 5: 986. https://doi.org/10.3390/w11050986

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop