2.1 Mathematical Formulation
In this study, a rectangular fluid domain with a horizontal length
L and constant depth
h was considered in a Cartesian coordinate system. In this coordinate system, the mean water surface was on the x-axis while the vertical z-axis was directed upwards. Assuming that the fluid is incompressible and inviscid and the flow is irrotational, the velocity of the fluid defined as the gradient of the velocity potential
ϕ and the continuity equation, which is the governing equation, is transformed into the Laplace equation as expressed in
Eq. (1).
where, ∇ denotes a horizontal gradient operator. According to the study conducted by
Zakharov (1968), the surface potential
φs at the free surface
ζ can be defined by
Eq. (2).
Using the surface potential
φs, the kinematic boundary conditions and dynamic boundary conditions of the free surface can be expressed by
Eqs. (3)–
(4).
where,
g and
W denote the acceleration due to the gravitational acceleration and vertical velocity of the free surface expressed in
Eq. (5), respectively.
The vertical velocity
W can be calculated using the HOS model presented by
Dommermuth and Yue (1987) and
West et al. (1987). The free surface
ζ and surface potential
φs based on to time progress can be obtained by integrating
Eqs. (3)–
(4), which substitutes the calculated vertical velocity
W. It is assumed that the horizontal region is infinite in the HOS model and the periodic boundary condition is adopted as expressed in
Eq. (6).
Eq. (7) was adopted in terms of the bottom boundary condition, which is the no-penetration condition condition.
The velocity potential, which is the solution of the Laplace equation that satisfies the aforementioned free surface boundary condition (
Eqs. (3)–
(4)), lateral periodic boundary condition (
Eq. (6)), and bottom boundary condition (
Eq. (7)), can be defined by
Eq. (8) using a Fourier pseudo-spectral basis.
where,
Cm denotes the time-dependent Fourier coefficient of the velocity potential, and
km denotes wave number defined by
Eq. (9), respectively.
Similar to the velocity potential, the surface potential
φs and free surface
ζ can be defined by a Fourier pseudo-spectral basis as expressed in
Eqs. (10)–
(11).
Where,
Am and
Bm respectively denote the time-dependent Fourier coefficients of the surface potential and free surface. The expansions of the surface potential and free surface in
Eqs. (10)–
(11) were limited to the number of points in the free surface that define the physical quantity. Moreover, the Fourier expansion in
Eqs. (10)–
(11) provided a horizontal differential to the time integrations of the surface potential and free surface using
Eqs. (3)–
(4). However, because the vertical velocity
W on the free surface is unknown, the vertical velocity
W of the free surface can be measured using the HOS model according to the expansion of
Dommermuth and Yue (1987). In the HOS model, the surface potential
φs, which is the physical quantity of the surface, and the free surface
ζ were used to evaluate the vertical velocity W of the free surface. First, the velocity potential
ϕ can be expressed up to the given order of
M with the perturbation expansion for wave steepness as expressed in
Eq. (12).
where, the perturbation potential
ϕ(m) has a perturbation order of
m for the wave steepness. From the mean depth of
z = 0, the perturbation potential
ϕ(m) for the free surface can be expressed by
Eq. (13) using the expansion of Taylor series.
The vertical velocity
W can also be expressed with perturbation expansion as in
Eq. (14), as presented in
Eq. (12).
where, the vertical perturbation velocity
W(m) also has a perturbation order of
m for the wave steepness and the vertical perturbation velocity
W(m) for the free surface
ζ from the mean depth of
z = 0 can be expressed with Taylor expansion as expressed in
Eq. (15).
The perturbation vertical velocity
W(m) in
Eq. (15) was calculated via the perturbation velocity potential
ϕ(m), which can be used to measure the vertical velocity
W at
ζ using
Eq. (14).
2.2 Numerical Technique
In this study, the free surface
ζ and surface potential
φs can be predicted depending on time integration using
Eqs. (3)–
(4), which define the boundary conditions of the free surface. As for time integration, the fourth order Runge–Kutta method, which ensured stability, was adopted as presented by
Dommermuth and Yue (1987). The physical quantity and vertical velocity W of the spatial differential included in the boundary condition of the free surface were calculated by applying the FFT algorithm. In this study, a numerical technique was implemented using MATLAB (Matrix laboratory). Because it supports the FFT algorithm and multi-thread, MATLAB provides fast computation. Regarding the HOS method, the multiplication of the nonlinear free surface boundary conditions was performed in a physical space rather than a spectra space. Although such transformation is efficiently performed using FFT, it leads to an aliasing phenomenon. It is therefore imperative to clearly eliminate the aliasing phenomenon to ensure the stability of the analysis owing to its influence on numerical stability. In this study, the zero-padding method was used to eliminate the aliasing phenomenon. Zero padding is a method for performing multiplication operations without aliasing in the physical space by filling in zero to expand the spectral space prior to multiplication. This is accompanied by the expansion of the spectral space, leading to an increase in computational cost as the consequence of an increase in the calculation domain. Therefore, appropriate spatial expansion was achieved via the half-rule as expressed in
Eq. (16) (Canuto et al., 1988;
Gatin et al., 2017).
where,
NExtend and
N represent an extended space for zero padding and the number of free modes used in HOS simulation. The parameter
p is completely independent of the aliasing phenomenon when the order
M is selected. An increase in computational cost owing to the expansion of the aforementioned spectral space can be expected as the order
M increases, as expressed in
Eq. (16). In this study, the full de-aliasing(
p =
M) was adopted because only two-dimensional calculations were performed.
The analysis of the HOS simulation is generally performed from the initial conditions of free surface
ζ and surface potential
φs at the free surface, in which the initialization is known to have significant influence on the numerical stability. According to the study by
Dommermuth (2000), the simulation of the HOS method under the linear initial condition causes numerical instability because it does not satisfy the nonlinear free surface boundary conditions in
Eqs. (3) and
(4). Therefore,
Dommermuth (2000) proposed an adjustment scheme to initiate the linear initial condition to the nonlinear free surface condition. This scheme enables the smooth transformation of the nonlinear solution via the relaxation of the nonlinear term of the free surface against time.
Eqs. (17)–
(18) express the free surface condition to which the relaxation of the nonlinear term was applied by introducing the relaxation time
Ta and relaxation index
n proposed by
Dommermuth (2000).
In this study, the linear free surface
ζ1st and linear potential
ϕ1st were adopted as the initial conditions of the regular wave under deep sea conditions, as expressed in
Eqs. (19)–
(20).
where,
A0,
k, and
ω denote the wave amplitude, wave number, and wave frequency of the linear solution. Concerning the initial condition of the irregular wave simulation, the free surface and linear potential of the irregular wave for the space can be defined by using the superposition of linear components. The wave amplitude for the wave frequency constituting the irregular wave from the defined wave spectrum
S(
ω) can be calculated using
Eq. (21).
By substituting the calculated wave amplitude
Ai into
Eqs. (22)–
(23), the initial conditions for the irregular wave simulation,
ζirregular and
φirregular can be calculated.
where
∊idenotes a uniform random number between zero and 2
π, and the wave number
ki and wave frequency
ωi are calculated via the dispersion relation. In this study, the joint North Sea wave project (JONSWAP) spectrum (
DNV, 2010) was used to define the wave spectrum, as expressed in
Eq. (24)
where
Hs,
Tp = 2
π/
ωp,
γ, and
Aγ denote a significant height, peak period, peak shape parameter, and normalizing factor, respectively. Furthermore,
σ is a spectral width parameter, which is determined depending on the peak wave frequency
ωp.