# Nonlinear Numerical Investigation on Higher Harmonics at Lee Side of a Submerged Bar

Authors: D. Ning; X. Zhuo; L. Chen; B. Teng
Journal: Abstract and Applied Analysis
 (2012)

---

## Abstract

The decomposition of a monochromatic wave over a submerged object is investigated numerically in a flume, based on a fully nonlinear HOBEM (higher-order boundary element method) model. Bound and free higher-harmonic waves propagating downstream the structure are discriminated by means of a two-point method. The developed numerical model is verified very well by comparison with the available data. Further numerical experiments are carried out to study the relations between free higher harmonics and wave nonlinearity. It is found that thenth-harmonic wave amplitude is growing proportional to the nth power of the incoming wave amplitude for weakly nonlinear wave condition, but higher-harmonic free wave amplitudes tend to a constant value for strong nonlinear wave condition.


---

## Body

## 1. Introduction

Submerged bars are common marine structures widely used in coastal engineering for different applications. Besides wave diffraction, higher-harmonic waves are generated due to wave shoaling effect, when the incoming waves propagate over a submerged structure. The generation of the superharmonic waves changes the swell spectrum because a significant part of the incoming wave energy may be transferred to higher frequencies. The free higher-harmonic waves are prominent at the lee side of the submerged bars and generally show strong nonlinear effect, whose wavelengths are much shorter than that of incident wave, wave amplitudes grow with increasing incoming wave amplitude, even to comparable to the incoming wave amplitude. Such free higher harmonics may debase sailing conditions and do harm to coastal structures. Thus, it is quite necessary to take account of their nonlinear characteristics.Numerous researching works have been performed on the problem of submerged bar, experimentally and numerically. For example, Beji and Battjes [1], Luth et al. [2] made physical experiments of wave transformation over a submerged bar and accordingly stated the generation mechanism of higher-frequency wave energy. Jeng et al. [3] experimentally studied the wave propagation over a submerged breakwater in front of a vertical wall, and the corresponding wave surfaces are recorded and analyses. Cho et al. [4] reported their experimental study of wave reflection over submerged breakwaters in tandem, and concluded that reflection coefficients of impermeable submerged bars are bigger than those of permeable ones and trapezoid bars are advised to use in the view of applicability. Rambabu and Mani [5] numerically studied the effects of incident wave factors, water depth over submerged bar, and section width of submerged bar on the transmission coefficient, and they accordingly obtained the optimal section width of submerged bar. Madsen et al. [6], Madsen and Sφrensen [7], Engsig-Karup et al. [8], Fang and Zou [9] adopted various Boussinesq-type equations to model wave propagation and transformation over submerged bars or bottom topography. However, the analysis of the higher harmonic waves scatting by submerged bar is relatively scare. Grue [10] experimentally researched the decomposition of higher harmonic waves due to submerged rectangular step and horizontal cylinder in deep water and got the conclusion that superharmonic wave amplitudes grow up to a saturation values with the increasing of incoming wave amplitude. Brossard and Chagdali [11], Liu et al. [12] experimentally and numerically studied the harmonic generation by waves over a submerged plate, respectively.Among various numerical studies, the Boussinesq-type equation has been commonly used in the analyses of nonlinear wave evolutions in shallow water. However, due to assumptions of both weak wave nonlinearity and weak wave dispersivity, the Boussinesq equation may not be valid for predictions of free higher harmonics in deep-water region and of highly nonlinear interaction between bound and free waves in the shallow-water region. On the other hand, applications of the higher-order boundary element method (HOBEM) to obtain a fully nonlinear numerical solution have become extremely successful in the past two decades. This method is achieved by directly solving the boundary integral equation with fully nonlinear free surface boundary conditions. The essential advantages include reducing the originaln-dimensional partial differential equation to an (n-1)-dimensional boundary integral equation, which allows the numbers of introduced unknowns to be reduced significantly, and no limitation to water depth and wave height if there do not exist wave breaking.In this paper, the free higher harmonics due to a monochromatic wave propagation over a submerged bar are addressed by using the authors’ developed fully nonlinear HOBEM [13]. This paper is arranged as follows. A brief description of the numerical wave tank model is given in Section 2, in which a two-point decomposition technique of higher harmonics from a periodic wave train is also introduced. In Section 3, the computed wave profiles at various locations and the higher harmonic wave amplitudes are compared to experimental observations in order to confirm the validity of the numerical model. The characteristics of free higher harmonics for weakly and strong nonlinear wave conditions are further addressed. Finally, the conclusions are presented in Section 4.


## 2. Mathematical Formulation

### 2.1. Governing Equation and Boundary Conditions

For a monochromic wave propagating over a submerged bar as shown in Figure1, the numerical analysis are performed using a numerical wave tank model previously developed by the authors [13]. In the model, the incident wave is generated through a nonreflective wave generator which combines vertically distributed wave-making sources on SI, introduced by Brorsen and Larsen [14]. Two damping layers are located at both ends of the computational domain to dissipate the outgoing waves.Figure 1
Sketch of the numerical wave tank model.Fluid is assumed to be ideal, so that a velocity potentialφ(x,y,z,t) exists and satisfies the following Poisson equation inside the fluid domain Ω:(2.1)∇2ϕ(x,y,z,t)=q*(xs,y,z,t),inΩ,
and the boundary conditions are expressed as(2.2)DX→F(t)Dt=∇ϕ-μ(x)(X-X0)DϕDt=12∇ϕ⋅∇ϕ-gη-μ(x)ϕ},on⁡Sf∂ϕ∂n=0,on⁡Sc,Sd,SB,
in which Sf is corresponding to the instantaneous free surface, Sd the tank bottom, SB the surface of submerged body, Sc the tank lateral surfaces. q*=2vδ(x-xs) is the pulsating volume flux density of the internal source distribution and v is prescribed second-order Stokes wave velocity. The damping function is adopted as the following:(2.3)μ(x)={αω(x-x1(2)Lb)2,(x>x1orx<x2),0,otherwise,
where x1(2) is a reference value specifying the at-rest position of the fluid particle, α is damping coefficient, Lb is length of damping zone. In this paper, α=1, Lb is one and half times incident wave length.


### 2.2. HOBEM

By using the second Green’s theorem, the prescribed boundary value problem can be transformed the following boundary integral equation:(2.4)α(p)ϕ(p)=∫S(ϕ(q)∂G(p,q)∂n-G(p,q)∂ϕ(q)∂n)dS+∫Ωq*G(p,q)dΩ,
where p and q are source and field points, and α(p) is the solid angle. For the Green function G, an image Green function is utilized so that the integration surface only includes the incident boundary SI and the free surface boundary and submerged bar SB. The Green function can be obtained by the superposition of the image of the Rankine source about the sea bed and the infinite images about the two lateral walls and written as (see [15])(2.5)G(ξ⃗,x⃗)=-14π[GC(x-x0,y-y0,z-zo)+GC(x-x0,y+y0,z-zo)+GC(x-x0,y-y0,z+zo+2h)+GC(x-x0,y+y0,z+zo+2h)],
where(2.6)GC(X,Y,Z)=1X2+Y2+Z2+∑n=-∞∞(1X2+(Y+2nB)2+Z2-12nB),B is the tank width, and h is the static water depth.Then the boundary surface is discretized with a number of elements. The geometry of each element is represented by the shape functions, thus the entire curved boundary can be approximated by a number of higher-order elements. Within the boundary elements, physical variables are also interpolated by the same shape functions, that is, the elements are isoparametric. In the integration process, the solid angle, the single-layer and double-layer integration are directly resolved (see [16]). Since the discretized integral equation is always variant in time, all the boundary surfaces are remeshed and updated at each time step using the mixed Eulerian-Lagrangian scheme and 4th-order Runga-Kutta approach. Once (2.4) is solved, we can obtain the time history of surface evaluation at any position as described by Ning and Teng [17].


### 2.3. Separation of Free Waves from Phase-Locked Waves

After waves pass the submerged bar, higher harmonic waves generated by nonlinear wave-wave interactions in the shallow water over the bar will leave the bar leeward as free waves. So the surface elevation atx at the lee side of the bar can be written as(2.7)η(t,x)=∑n=1∞an(F)cos⁡(knx-nωt+ψn(x))+∑n=2∞an(L)cos⁡(n(kx-ωt+ψ1(x))),
where an(F) is the amplitudes of the free transmitted waves with frequencies of integer times of the incident wave frequency; an(L) are the amplitudes of the nth-harmonic phase-locked waves, ψ1(x) are the initial phase angles of the fundamental waves, and ψn(x)(n≥2) are the nth-harmonic free waves; k and kn are the wave number of the fundamental waves and the nth-harmonic free waves, and they satisfy the following dispersion relations:(2.8)ω2=gktanh⁡kh,(nω)2=gkntanh⁡knh,n=2,3,…,
respectively. The fundamental wave amplitude, as well as the higher-order free and locked wave amplitudes, is obtained from the time histories of the surface elevation. The Fourier transform is introduced as follows:(2.9)ηn(x)=1T∫0Tη(x,t)e-inwtdt=An(x)+iBn(x),
where An(x) and Bn(x) are the corresponding real and imaginary components, respectively.Substituting (2.7) into (2.9) and making use of orthogonality of trigonometric function, we may obtain(2.10)a1(F)=A12(x)+B12(x),ψ1(x)=actanA1(x)B1(x)-kx,(2.11)an(L)cos⁡(kx+ψ1(x))+an(F)cos⁡(knx+ψn(x))=An(x),an(L)sin⁡(kx+ψ1(x))+an(F)sin⁡(knx+ψn(x))=Bn(x).Then, we can obtain the amplitudes of thenth-harmonic locked waves and free waves after applying two fixed points’ (Δx apart) surface elevations into (2.11) through an easy iterative solving process.


## 2.1. Governing Equation and Boundary Conditions

For a monochromic wave propagating over a submerged bar as shown in Figure1, the numerical analysis are performed using a numerical wave tank model previously developed by the authors [13]. In the model, the incident wave is generated through a nonreflective wave generator which combines vertically distributed wave-making sources on SI, introduced by Brorsen and Larsen [14]. Two damping layers are located at both ends of the computational domain to dissipate the outgoing waves.Figure 1
Sketch of the numerical wave tank model.Fluid is assumed to be ideal, so that a velocity potentialφ(x,y,z,t) exists and satisfies the following Poisson equation inside the fluid domain Ω:(2.1)∇2ϕ(x,y,z,t)=q*(xs,y,z,t),inΩ,
and the boundary conditions are expressed as(2.2)DX→F(t)Dt=∇ϕ-μ(x)(X-X0)DϕDt=12∇ϕ⋅∇ϕ-gη-μ(x)ϕ},on⁡Sf∂ϕ∂n=0,on⁡Sc,Sd,SB,
in which Sf is corresponding to the instantaneous free surface, Sd the tank bottom, SB the surface of submerged body, Sc the tank lateral surfaces. q*=2vδ(x-xs) is the pulsating volume flux density of the internal source distribution and v is prescribed second-order Stokes wave velocity. The damping function is adopted as the following:(2.3)μ(x)={αω(x-x1(2)Lb)2,(x>x1orx<x2),0,otherwise,
where x1(2) is a reference value specifying the at-rest position of the fluid particle, α is damping coefficient, Lb is length of damping zone. In this paper, α=1, Lb is one and half times incident wave length.


## 2.2. HOBEM

By using the second Green’s theorem, the prescribed boundary value problem can be transformed the following boundary integral equation:(2.4)α(p)ϕ(p)=∫S(ϕ(q)∂G(p,q)∂n-G(p,q)∂ϕ(q)∂n)dS+∫Ωq*G(p,q)dΩ,
where p and q are source and field points, and α(p) is the solid angle. For the Green function G, an image Green function is utilized so that the integration surface only includes the incident boundary SI and the free surface boundary and submerged bar SB. The Green function can be obtained by the superposition of the image of the Rankine source about the sea bed and the infinite images about the two lateral walls and written as (see [15])(2.5)G(ξ⃗,x⃗)=-14π[GC(x-x0,y-y0,z-zo)+GC(x-x0,y+y0,z-zo)+GC(x-x0,y-y0,z+zo+2h)+GC(x-x0,y+y0,z+zo+2h)],
where(2.6)GC(X,Y,Z)=1X2+Y2+Z2+∑n=-∞∞(1X2+(Y+2nB)2+Z2-12nB),B is the tank width, and h is the static water depth.Then the boundary surface is discretized with a number of elements. The geometry of each element is represented by the shape functions, thus the entire curved boundary can be approximated by a number of higher-order elements. Within the boundary elements, physical variables are also interpolated by the same shape functions, that is, the elements are isoparametric. In the integration process, the solid angle, the single-layer and double-layer integration are directly resolved (see [16]). Since the discretized integral equation is always variant in time, all the boundary surfaces are remeshed and updated at each time step using the mixed Eulerian-Lagrangian scheme and 4th-order Runga-Kutta approach. Once (2.4) is solved, we can obtain the time history of surface evaluation at any position as described by Ning and Teng [17].


## 2.3. Separation of Free Waves from Phase-Locked Waves

After waves pass the submerged bar, higher harmonic waves generated by nonlinear wave-wave interactions in the shallow water over the bar will leave the bar leeward as free waves. So the surface elevation atx at the lee side of the bar can be written as(2.7)η(t,x)=∑n=1∞an(F)cos⁡(knx-nωt+ψn(x))+∑n=2∞an(L)cos⁡(n(kx-ωt+ψ1(x))),
where an(F) is the amplitudes of the free transmitted waves with frequencies of integer times of the incident wave frequency; an(L) are the amplitudes of the nth-harmonic phase-locked waves, ψ1(x) are the initial phase angles of the fundamental waves, and ψn(x)(n≥2) are the nth-harmonic free waves; k and kn are the wave number of the fundamental waves and the nth-harmonic free waves, and they satisfy the following dispersion relations:(2.8)ω2=gktanh⁡kh,(nω)2=gkntanh⁡knh,n=2,3,…,
respectively. The fundamental wave amplitude, as well as the higher-order free and locked wave amplitudes, is obtained from the time histories of the surface elevation. The Fourier transform is introduced as follows:(2.9)ηn(x)=1T∫0Tη(x,t)e-inwtdt=An(x)+iBn(x),
where An(x) and Bn(x) are the corresponding real and imaginary components, respectively.Substituting (2.7) into (2.9) and making use of orthogonality of trigonometric function, we may obtain(2.10)a1(F)=A12(x)+B12(x),ψ1(x)=actanA1(x)B1(x)-kx,(2.11)an(L)cos⁡(kx+ψ1(x))+an(F)cos⁡(knx+ψn(x))=An(x),an(L)sin⁡(kx+ψ1(x))+an(F)sin⁡(knx+ψn(x))=Bn(x).Then, we can obtain the amplitudes of thenth-harmonic locked waves and free waves after applying two fixed points’ (Δx apart) surface elevations into (2.11) through an easy iterative solving process.


## 3. Numerical Results and Discussions

As a validation of the present model, a case for monochromic wave propagating a submerged trapezoid bar as shown in Figure2 was considered firstly. The same problem was ever experimentally studied by Beji and Battjes [1] and Luth et al. [2]. In the present numerical simulation, some parameters for wave period T=2.02s, wave amplitude A=0.01m, and water depth h=0.4m are used. The corresponding computational domain is taken as 10λ×0.12λ (λ=2π/k denotes wave length), time step defined as T/80, and meshed with 200×2×10 cells in x, y and z directions after numerical convergent tests, in which the first and last 1.5λ are used as the damping layers.Figure 2
Definition sketch of wave propagating over a trapezoid bar.Figure3 presents the distribution of wave elevations along the plane of symmetry at t=18T and 20T. From the figure, it can be seen that there is a good agreement between these two times, which states that the numerical convergence has been obtained in the proposed model. It is also found that the wave profiles at the weather side of the submerge object are still regular and periodic, but they are quite skewed at the lee side. It means that a number of wave components are introduced the incident flow at the lee side.Figure 3
Wave profile along the plane of symmetry att=18T and 20T.The comparisons of wave elevation time series at specific points marked in Figure2 among the present results, experimental data [1], and the solutions of the extended Boussinesq model of Nwogu [9] are given in Figure 4. From the figures, it can be seen that there are good agreements for the three methods at points x=2.0m,10.5m and 12.5m. But for the points at the lee side such as x=19.0m, the Boussinesq model cannot give accurate description because the wave nonlinearity is stronger behind the shoal. Therefore, the fully nonlinear model, such as the present scheme, is advised to analyze such strong nonlinear wave problem.Comparisons of free surface displacement among the experiment, Boussinesq model, and present method.
(a)
x=2.0m(b)
x=10.5m(c)
x=12.5m(d)
x=15.7m(e)
x=17.3m(f)
x=19.0mThe results in Figures2 and 3 show that much higher-harmonic waves are involved at the lee side of the submerged bar. So it is necessary to analyze these higher harmonics. Equation (2.7) is used to separate various waves at point x=25.75m by truncating n to 7. The downstream neighboring point is also used here for solving (2.11). Figure 5 gives the comparison of free surface at point x=25.75m between the numerical solution and the fitting curve for weakly nonlinear condition (A=0.05m, T=2.02s, h=0.8m) and strong nonlinear condition (A=0.01m, h=0.4m), respectively. The other parameters are the same with the above. It proves the validity of the present two-point separation method from good agreement. Figure 5(a) shows a regular wave profile similar to incident wave. But Figure 5(b) presents a strong deformed wave profile, even secondary wave peaks can be apparently found.Comparison of free surface atx=25.75m between the numerical solution and the fitting curve.
(a)
Weakly nonlinear condition (A=0.05m, h=0.8m)(b)
Strong nonlinear condition (A=0.01m, h=0.4m)Once the parametersan(L), an(F) and initial phase angles are solved, the various order wave records can be obtained as shown in Figures 6 and 7. From the figure, it can be seen that the periods of second and third harmonic waves are corresponding to 1/2 and 1/3 times wave period, and phase-locked waves and free waves keep same periods but different phases. Locked waves are released to the corresponding free wave in a large degree after the wave propagating over the submerged bar, especially for strong nonlinear condition shown in Figure 7. With the increasing of wave nonlinearity, the higher harmonic free waves are found to be more prominent at the obstacle’s lee side. The second-harmonic free wave amplitude is comparable to the incoming wave amplitude as shown in Figure 7.Time histories of various order harmonic waves for weakly nonlinear condition (A=0.05m, h=0.8m).
(a)
1st-harmonic free wave(b)
2nd-harmonic waves(c)
3rd-harmonic wavesTime histories of various order harmonic waves for strong nonlinear condition (A=0.01m, h=0.4m).
(a)
1st-harmonic free wave(b)
2nd -harmonic waves(c)
3rd-harmonic wavesFrom above results, it shows that the incoming wave is mainly broken up into superharmonic free waves at the lee side of submerged structure due to nonlinear free surface effect. Thus, the higher harmonic free waves are especially researched in the following section. Figure8 gives the distribution of first-, second-, and third-harmonic free wave amplitudes with the increasing of incident wave amplitudes for static water depth h=0.8m. In the figure, the real lines are the linear, quadratic, and cubic fitting curves obtained from the computed data. It means that the nth-harmonic wave amplitude is growing proportional to the nth power of the incoming wave amplitude for weakly nonlinear wave condition. Figure 9 shows the distribution of first-, second-, and third-harmonic free wave amplitudes with the increasing of incident wave amplitudes for static water depth h=0.4m. With the water depth h decreased, the wave nonlinearity is increased. The relation between harmonic free wave amplitudes and incoming wave amplitudes does not obey the rules shown in Figure 8 again. It is observed that the generation of higher harmonic free waves is more powerful than that in Figure 8. The various harmonic free wave amplitudes tend to a constant value as the incoming wave amplitude A≥0.025m. Surprisingly, the portion of the second-harmonic free wave amplitude is even larger than that of first-harmonic wave amplitude since A≥0.01m.First-, second-, and third-harmonic free wave amplitudes at the lee side of submerged bar forh=0.8m.
(a)
Fundamental free wave(b)
2nd- and 3rd-harmonic free wavesFirst-, second-, and third-harmonic free wave amplitudes at the lee side of submerged bar forh=0.4m.
(a)
Fundamental wave(b)
2nd- and 3rd-harmonic free wavesTo further testify the conclusion obtained in Figure9, another case for monochromic wave propagating a submerged rectangular step is considered. The rectangular shelf has a cross-section 0.5 m long and 0.41 m high. The distance h0 between the uppermost point of the geometry and the mean free surface is chosen as 0.0375 m and 0.05 m, respectively. The incident wave frequency w is defined as 5.969 and 6.597 rad/s, respectively. Figure 10 shows the distribution of higher free harmonic wave amplitudes with different incident wave amplitudes, in which the experimental data from Grue [10] are also given and the longitudinal coordinates are dimensionless divided by the incident wave amplitude A. From the figure, it can be seen that there are better agreements between numerical solutions and experimental data. The similar conclusions with those in Figure 9 are also drawn.Distribution of higher free harmonic wave amplitudes with incident wave amplitude.
(a)
(h0=0.0375m, w=5.969rad/s)(b)
(h0=0.05m, w=6.597rad/s)


## 4. Conclusions

The phenomenon of wave propagation over a submerged bar is examined by a powerful numerical wave tank method. The transmitted waves are separated into the free waves and the locked waves with different harmonic frequencies by a two-point method. Numerical results show that, at the leeside of the bar, the portions of the free waves are much larger than those of the locked waves with the same frequency. For weakly nonlinear wave conditions, thenth-harmonic free wave amplitude is growing proportional to the nth power of the incoming wave amplitude and the fundamental free wave plays an important role. However, for strong nonlinear wave condition, the fundamental free wave amplitude is decreased and second- and third-harmonic free wave amplitudes are increased with the increasing of incoming wave amplitude, and all tend to a constant value, that is, a kind of saturation statues. Even the phenomena of higher-harmonic free wave amplitude exceeds the lower-harmonic free wave amplitude occur.



---

*Source: 214897-2012-04-03.xml*