For an arbitrary-shape harbor, the domain is divided into two regions as shown in Fig. (1): Region A is the limit of the harbor and Region R is the infinite open-sea. ∂A is near island region as shown in Fig. (1), and also in the harbor basin region shown in Fig. (2), ∂C is the breakwater boundary shown in Fig. (2). The governing equations for the respective regions can be expressed as:

Fig. (1). Configuration of an Offshore Harbor.

(3)

and

(4)

Where, water depth is constant or deep.

The total velocity potential can be separated into incident (Φ^I), reflected (Φ^R), and scattered (Φ^S) wave potentials, such that Φ = Φ^I + Φ^S (for offshore harbor), and Φ = Φ^I + Φ^R + Φ^S (for coastal harbor).

The incident wave potential can be expressed as:

(5)

where ε₀ = 1, ε_n = 2, n = 1, 2…. are the Neumann factors, a₀ is the amplitude of incident waves, θ^I is the incident angle, J_n is the Bessel function of the first kind and order n, θ is the angle variable in polar coordinates, and r is the radial coordinate.

For a semi-infinite ocean with a straight coastline, the reflected wave potential can be written as:

(6)

The scattered wave potential satisfies the Helmholtz equation and Sommerfeld radiation condition, which can be written as:

(7)

where H_n¹ is the Hankel function of the first kind of order n, α_n and β_n are unknown coefficients yet to be determined. The boundary conditions fall into the following domains:

(1) The domain B (Fig. 1) represents the harbor, along whose boundaries appropriate boundary conditions must be specified. On a solid boundary, the normal velocity must vanish, i.e. ∂Φ/∂n = 0, where n is the outward normal vector to the boundary.

In practical problems, reflection may be partial where energy absorption occurs. Therefore, the general boundary condition can be represented as:

(8)

Where, α is a boundary absorbing coefficient that must be empirically prescribed.

(2) In the infinite sea region R, the boundary ∂R is taken as a circular or semicircular contour. A Sommerfeld radiation condition is imposed at infinity in all horizontal directions and the scattered wave Φ^S must be an outgoing wave at infinity:

(9)

(3) A matching continuity condition is imposed along the boundary that may be expressed as:

(10)

(11)

The present method derives a function from the principle of conservation of energy. The total wave energy in the inner Region A,E_A (total), is the sum of the total wave energy in Region A,(E_A) energy loss from boundary absorption dissipation (E_f), bottom friction energy loss (E_bf), and energy transmission through the permeable breakwater (E_Tf), which may be expressed as:

(12)

In the outer Region R:

So the total wave energy can be expressed as:

(13)

For a progressive long wave, the average total energy may be obtained from the contributions of both potential energy (PE) and kinetic energy (KE).

(1) The potential wave energy for per unit wave length is:

(14)

where ρ is the water density, g is the acceleration due to gravity, and . The mean free surface is located at z = 0 and the free surface elevation η can be derived from the linearized Bernoulli’s equation and the dispersion relation as:

(15)

Substituting Eq.(15) into (14), the simplified potential wave energy equation is:

(16)

(2) The kinetic wave energy for per unit wave length is:

(17)

where the velocity potential is given as , the velocity potential for monochromatic waves with angular frequency ω, is . Substituting into Eq.(17) results in:

(18)

Then the simplified kinetic wave energy equation is:

(19)

where .

Total wave energy is equal to potential (Eq.16) and kinetic (Eq.19) components, expressed as:

(20)

The instantaneous energy flux of a progressive long wave passing through a section of unit width, characterized by its horizontal normal n, can be expressed as:

(21)

where P is the complex excess pressure, , and U_c is the complex particle velocity in the normal direction, ; then,

(22)

A portion of wave energy dissipation is caused by bottom friction. Bottom energy flux per unit area (E_fb) can be expressed as:

(23)

where U_b is the instantaneous particle velocity at the bottom.

Conventionally the shear stress is defined as:

(23a)

where (τ_b,max) real and (U_b,max)real² and k_b are the real-valued maximum shear stress, potential theory particle velocity, and wave friction at the bottom, respectively. Jonsson [25]

For convenience, we define:

(24)

where and has dimensions of ft²/sec³

Particle velocity at the bottom (U_b) can be expressed as:

(25)

(26)

Substituting Eqs. (24) and (26) into Eq. (23) results in a bottom friction energy flux equation:

(27)

(27a)

The wave partial reflection due to wave energy dissipation at the absorbing boundaries may exist such as the instantaneous wave energy flux at the absorbing boundaries. E_f^A is the instantaneous wave energy flux, measured positive, from A out through ∂B to ∂A. One can integrate the instantaneous energy flux as Eq.(22) with respect to time then yield:

(28)

Because of the absorbing boundary condition Eq.(8), , then along the absorbing boundaries ∂B can be written as:

(29)

Permeable harbor structures, such as rubble-mound breakwaters, have been widely constructed to provide protection from wave attack by reflecting and dissipating incident wave energy. At a permeable breakwater only part of the incident wave energy is reflected back to the ocean, while the rest is absorbed and/or transmitted into the harbor through the structure. If this transmitted energy is exceedingly high, it may create undesirable and unsafe wave conditions inside the harbor.

Consider a two dimensional interaction of a monochromatic wave train with a single homogeneous, isotropic, porous structure of width located between two semi-infinite fluid regions, as shown in Fig. (2). The damping region occupies 0 < x < l and the flow domain extends to infinity (x → ± ∞). The wave field can be divided into three regions:

Region 1: x < 0, Energy dissipation can occur (incident wave region)

Region 2: 0 < x < l, Energy dissipation can occur (absorbed and transmitted wave region)

Region 3: x > l, No energy dissipation (transmitted wave region)

Since the solution in adjacent regions must be continuous at each interface, the appropriate boundary conditions are the continuity of pressure and horizontal mass flux at y = 0 (between region 1 and 2) and y = l (between region 2 and 3), expressed as,

(30)

(31)

where ε, ε^' are the porosity of the porous medium, the solution of the velocity potential in each region (similar to Liu, et. al. [26]) can be written as:

(32)

(33)

(34)

where a₀ is the incident wave amplitude, B and F are unknown constants to be determined, R is the reflection coefficient, and T represents the transmission coefficient.

The wave transmission through a porous structure, such as the instantaneous wave energy flux at Region 3 (Fig. 2), can be written as:

Fig. (2). Configuration of a coastline harbor with permeable breakwater.

(35)

From Eqs. (32) and (34), assuming that Φ₃ = k_TΦ_I, results in,

(36)

where Φ_I is the incident wave potential, and k_T is an empirical transmission coefficient, Eq.(35) can then be written as:

(37)

Combining total wave energy (Eq.20), bottom friction (Eq.27), absorbing boundary (Eq.29), and transmitting boundary (Eq.37) terms, a functional can be obtained from Eqs.(12) and (13) such that, the stationary functional is then:

(38)

The discretization of terms I₁ to I₈ will be discussed in detail in the following section.

The interpolation function is used to evaluate the functional, and the Gaussian quadrature four-point method of numerical integration is used to solve for phase velocity and group velocity. Using this approach the functional, Eq. (38), can be evaluated as follows:

(1) First term (energy in region A):

(39)

with

(40)

where Wn is the weighting function, and

(2) Second term (energy flux along ∂A):

(41)

with

(3) Third term (energy flux along ∂A):

(42)

with

(42a)

(4) Fourth term (energy flux along ∂A):

(43)

with

(43a)

(5) Fifth term (energy flux along ∂A):

(44)

with

(44a)

(6) Sixth term (boundary absorption energy flux along ∂B):

(45)

with

(45a)

where L is the total number of line elements along the inner boundary ∂B, L_B^e is the length of each segment, and α is a boundary absorption coefficient, dependent on empirical and laboratory data.

(7) Seventh term (bottom friction energy flux along A):

(46)

with

(46a)

where W_n is a weighting factor and Kb is an empirical bottom friction coefficient.

(8) Eighth term (energy transmitting through porous structures along ∂C):

(47)

with

(47a)

where L^' is the total number of line elements along the porous structures (breakwater) ∂C, L_C^'e is the length of each segment, and k_T is the transmission coefficient.

Assembling the local stiffness matrix into a global stiffness matrix results in the functional being written as:

(48)

By combining nodal unknowns and coefficient unknowns and defining the total unknown vector {ϕ} by

(49)

where N=M+M and

(50)

the functional can be written as:

(51)

where the total load vector

(52)

and the total stiffness matrix [K] is usually large and banded,

(53)

then, the stationary functional of F(Φ) implies that

(54)

Substitution of equation (54) into equation (51), leads to

(55)

This set of N linear algebraic equations for N unknowns can be solved by the LU elimination method. The total velocity potential at all nodal points and unknown coefficients can also be determined.

A fully-open rectangular harbor finite element grid layout is shown in Fig. (3). The present model results are compared with existing experimental data of Ippen and Goda [27], Lee [1], and Chen & Mei’s [2, 3] model, as shown in Fig. (4). The incident wave angle θ^I = -π/2 the present model with boundary absorption α = 0.07, bottom friction coefficient K_b = 0.3. The resulting curve is the most exact replica of the experimental results. It reveals an accurate agreement between the numerical results and experimental data.

Fig. (3). Finite element grid layout for rectangular harbor resonance model.

Fig. (4). Response of amplification factor (R) for different wave number (Kl) at the end wall of the fully-open rectangular harbor.

Fig. (5) shows the effect of energy dissipation due to bottom friction. For deep harbor (h = 10.128 in.) case as shown in top Fig. (5), the reduction in wave amplitude inside the harbor is quite small; for shallow harbor (h=0.1 in.) case as shown in bottom Fig. (5), the wave response is quite sensitive.

Fig. (5). Comparison of the maximum response curve for a fully-open rectangular harbor using 4 different bottom friction coefficients (fω).

In order to test the effect of energy dissipation due to fully absorbing boundaries, as shown in Fig. (6), when the boundary absorption coefficient (α) increases, it shows an increase in the highly sensitive response curves reactively. In order to express more effectively using three-dimensional plots, Figs. (7 and 8) show a three dimensional plot of the relative surface elevations for the harbor resonance response of the first and second resonance mode, for the inner boundary absorption coefficients α equal to 0.0, 0.1, 0.3, and 0.5 respectively.

Fig. (6). Comparison of the maximum response curve for a fully-open rectangular harbor using 4 different boundary absorption coefficients (α).

Fig. (7-1). The relative surface elevations for the harbor resonance response of the first resonance mode.

In a realistic harbor, which consists of breakwaters, quay walls, and wharfs, boundary reflection factors may be quite different. Fig. (8) shows the effect of partial boundary absorption (P.B.A.) on R_max for five different boundary absorption conditions: (1) No boundary absorption, (2) One end wall with P.B.A. (4 elements, α = 0.3), (3) One side wall with P.B.A. (12 elements, α = 0.3), (4) One side and one end walls with P.B.A. (16 elements, α = 0.3), (5) All three walls with P.B.A. (27 elements. α = 0.3). The results seem reasonable and point to the fact that as the total area of absorption (i.e. number of absorbing boundary elements) increases, the harbor response should be reduced.

A wider rectangular harbor with breakwater, with its finite element grid layout shown in Fig. (9) reveals the effect of permeable transmitting boundaries for three cases as shown in Fig. (10). The harbor with an impermeable breakwater creates a larger response as the energy inside the harbor can’t be reflected back to sea. Apparently a harbor with permeable breakwater has more energy flux transmitted through the porous structure into the harbor, and thus produces the largest response, even larger than that produced for a harbor with impermeable breakwater. Fig. (11) shows variation of maximum amplification (R_max) with a different transmission coefficient (K_T) for the first, second, and third resonant modes for a harbor with permeable breakwater. Larger amplification occurs for (K_T) and smaller wave number (Kl). When K_T is within the range 0.3 to 0.8, a larger amplification is produced in the first and second resonant modes.

Fig. (7-2). The relative surface elevations for the harbor resonance response of the second resonance mode.

Fig. (8). The test result of partial boundary absorption (P.B.A.) on Rmax for five different boundary absorption conditions.

Fig. (9). A finite element grid layout of a rectangular harbor with breakwater.

Fig. (10). Comparison of the maximum response curve for a rectangular harbor (13X8X10.128) using three cases: (1)Fully-open harbor, (2)With impermeable breakwater 2 inch entrance long, (3) With permeable breakwater KT = 0.3.

Fig. (11). Comparison of maximum amplification factor (Rmax) with different transmission coefficient (KT) for the first, second, and third resonant modes for a harbor with permeable breakwater.

Fig. (12) shows the effect of variable harbor entrance using three cases: (1) Fully open entrance, (2) Partial 4 inch entrance width, (3) Partial 2 inch entrance width. As the entrance width decreases, the response curve amplitude increases and Kl decreases. The model thus confirms the Harbor Paradox of Miles and Munk [28].

Fig. (12). Comparison of 3 different entrance width cases (D=8 in., 4 in., 2 in.).

Fig. (13) shows the effect of varying length and width for a rectangular harbor with a constant basin area and entrance width, using four different horizontal dimensions. As the basin length decreases and width increases, the response curve amplitude for the second resonant mode increases and has a larger Kl.

Fig. (13). Comparison of the maximum amplification factor (Rmax) and the harbor width to length ratio (W/L) for the first and second resonant modes.