### 1. Introduction

### 2. Mathematical Formulation and Numerical Technique

### 2.1 Mathematical Formulation

*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

*φ*at the free surface

^{s}*ζ*can be defined by Eq. (2).

*φ*, the kinematic boundary conditions and dynamic boundary conditions of the free surface can be expressed by Eqs. (3)–(4). where,

^{s}*g*and

*W*denote the acceleration due to the gravitational acceleration and vertical velocity of the free surface expressed in Eq. (5), respectively.

*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

*φ*based on to time progress can be obtained by integrating Eqs. (3)–(4), which substitutes the calculated vertical velocity

^{s}*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.

*C*denotes the time-dependent Fourier coefficient of the velocity potential, and

_{m}*k*denotes wave number defined by Eq. (9), respectively.

_{m}*φ*and free surface

^{s}*ζ*can be defined by a Fourier pseudo-spectral basis as expressed in Eqs. (10)–(11). Where,

*A*and

_{m}*B*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

_{m}*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

*φ*, which is the physical quantity of the surface, and the free surface

^{s}*ζ*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.

*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).

### 2.2 Numerical Technique

*ζ*and surface potential

*φ*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,

^{s}*N*and

_{Extend}*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.

*ζ*and surface potential

*φ*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

^{s}*T*and relaxation index

_{a}*n*proposed by Dommermuth (2000).

*ζ*

_{1}

*and linear potential*

_{st}*ϕ*

_{1}

*were adopted as the initial conditions of the regular wave under deep sea conditions, as expressed in Eqs. (19)–(20). where,*

_{st}*A*

_{0},

*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).

*A*into Eqs. (22)–(23), the initial conditions for the irregular wave simulation,

_{i}*ζ*and

_{irregular}*φ*can be calculated. where

_{irregular}*∊*denotes a uniform random number between zero and 2

_{i}*π*, and the wave number

*k*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

_{i}*H*,

_{s}*T*= 2

_{p}*π*/

*ω*,

_{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}### 3. Numerical Simulation

### 3.1 Adjusted Simulation of Single Wave

*kA*= 0.1 and presented the convergence result of the eighth harmonic component of the Stokes wave. In this study, the developed code was verified by performing calculation conditions, as presented in Table 1, which are the same conditions as those of Dommermuth (2000); after that, the results were compared.

*t*at 0.025 intervals. In the HOS simulation, the linear solution converged to a nonlinear Stokes wave, as illustrated in Fig. 1.

### 3.2 Irregular Wave Simulation

*M*to 1, 3, 5, and 8 for the same initial condition. Fig. 4 presents the normalized maximum wave height and kurtosis (

*κ*) based on the change in time. When the normalized maximum wave height exceeds 2 m, it is called an abnormality index (AI). It is known that the abnormality index and kurtosis are highly correlated with the occurrence of Freak wave (Kim, 2019). The abnormality index was calculated using

*H*

_{max}and

*H*

_{1/3}was determined via the zero crossing analysis of the wave distribution at every time step. The kurtosis was a measure adopted to statistically establish the difference with the normal distribution, as defined in Eq. (25). where

*σ*and

_{ζ}*ζ*denote the root mean square (RMS) of the spatial wave distribution of the wave height and the mean of the spatial wave distribution, respectively. It was inferred that as time progresses, a significant difference exists between the abnormality index and kurtosis of orders 3, 5, and 8 that reflect the nonlinear effect of the free surface when the order is one, as illustrated in Fig. 4. This phenomenon occurs because of the nonlinear interaction between the wave components. Although the abnormality index and kurtosis of orders 3, 5, and 8 developed similarly up to 260

_{mean}*T*, the abnormality index and kurtosis of orders 5 and 8 exhibited different results from that of order 3 after 260

_{P}*T*. This indirectly showed that the order needs to be 5 or above when simulating 260

_{P}*T*or above of the corresponding condition. Moreover, it was inferred that the normalized maximum wave height and kurtosis developed in a similar pattern, thus indicating that the two physical quantities are highly correlated.

_{P}### 4. Numerical Technique and Numerical Simulation for Wave Generation

### 4.1 Numerical Techniques for Generating Waves

*L*was required as the wave steepness increased.

_{NAZ}*x*and

_{a}*x*,

_{b}*γ*, and

_{damping}*λ*denote the boundary points of the damping zone, wave damping coefficient, and length of the damping zone, respectively. In this study,

_{damping}*γ*was set to one and

_{damping}*λ*was set at a wavelength longer than the generated wavelength. If the nonlinear adjustment and damping zones are applied to the free surface boundary condition of Eqs. (3)–(4), they are defined as expressed in Eqs. (28)–(29).

_{damping}*x*was equal to zero, and the wave transfer function was calculated using the harmonic simulation results of white noise under the boundary condition of the linear free surface, as illustrated in Fig. 6.

### 4.2 Wave Generation Numerical Simulation

*ζ*denotes the disturbance amplitude of the wave at the point where

_{influx}*x*= 0.

*k A*was 0.025 and 0.126. The simulation that did not apply the nonlinear adjustment zone was additionally performed to understand the characteristics of the technique.

*kA*is 0.025, it indicates the condition for generating waves with weak nonlinearity, where the amplitude difference between orders 1 and 5 is insignificant, as demonstrated in Fig. 9(a). In contrast, the condition when

*kA*is 0.126 indicates the generation of waves with strong relative nonlinearity, which shows nonlinear characteristics with a high crest and steady trough, as demonstrated in Fig. 9(b). To strictly observe these nonlinear characteristics, it was compared with the 5th order Stokes wave (Fenton, 1985), as illustrated in Fig. 10. It was inferred that the shape of the generated wave did not only agree well with the 5th order Stokes wave, but also the waves of desired periods and amplitudes can be generated via the wave transfer function determined in advance, as demonstrates in Fig. 8. Regarding the application of nonlinear free surface boundary conditions without the nonlinear adjustment zone, the amplitude and phase of the wave differed from the linear result when

*kA*was 0.025, and the simulation of regular wave generation diverged when

*kA*was 0.126, as illustrated in Fig. 9(a). Consequently, it was concluded that it is necessary to apply a nonlinear adjustment zone when simulating wave generation via wave disturbance. It was determined that the HOS-based numerical harmonic tank can be implemented using a numerical technique such as nonlinear adjustment and damping zones based on the above results.

*x*= 1252 m, 2505 m, 3750 m, and 4994 m).

*x*was equal to 1252 m, which is the closest point after the relaxation distance of the nonlinear adjustment zone during wave generation; however, the difference in the input wave spectrum increased as the wave propagated away from the origin of wave generation. Consequently, it can be inferred that the nonlinear interaction of the wave increased compared to the previously described wave steepness condition. Based on the comparison with the mean wave steepness of 0.05, it was deduced that the energy redistribution between wave components was more significant owing to the nonlinear interaction between wave components during wave propagation. Some differences between the results of orders 3 and 5 were determined, thus indicating that for the mean wave steepness of 0.1, nonlinear interactions can be fully reflected only when the order is considered to be 5 or above similar to the results of the irregular wave simulation in Section 3.2.

### 5. Conclusion

*W*was estimated based on the HOS order according to the expansion presented by Dommermuth and Yue (1987). Furthermore, the time integration for the nonlinear free surface boundary condition was stably performed using the Runge-Kutta 4th order method accordingly. The aliasing phenomenon caused by the multiplication of the nonlinear free surface boundary condition was eliminated using the zero padding method, thus securing numerical stability. The conventional HOS method is limited in its simulation of wave generation in a confined space, such as in a tank, with lateral periodic conditions. Therefore, in this study, a nonlinear adjustment and damping zones were introduced to implement a numerical tank capable of generating waves via surface disturbance. Several numerical simulations were performed using the developed code, therefore, the following conclusions were derived based on the results of performed simulations:

By performing the same nonlinear relaxation simulation as in the study by Dommermuth (2000), it was inferred that the linear solution converged to the eighth order Stoke wave solution. The HOS code developed in this study was verified by demonstrating that the convergence pattern was the same as that of Dommermuth (2000) and the convergence result was highly consistent with the analysis solution.

The abnormal index and kurtosis over time were observed via the simulation of an irregular wave field with a mean wave steepness of 0.1. A high correlation between the two physical quantities was determined based on the similarity of the time series development pattern between the abnormality index and kurtosis, which are known to exhibit high correlation with the occurrence of freak waves. It was inferred that to fully reflect the nonlinear interaction between waves in this mean wave steepness 0.1, order 5 or above needs to be considered.

In the simulation of an irregular wave field, it was determined that the difference between the wave spectrum given as the initial condition and that of the final time step increased as the order increased. It was therefore concluded that energy was redistributed between wave components via the change in the wave spectrum, which can be deduced to have resulted from the nonlinear interaction between wave components.

It was reconfirmed that the HOS method can efficiently simulate the development of irregular waves over a long period of time within a limited grid because the wave that outflowed from the ends was re-introduced through the starting point owing to the periodic boundary conditions at both ends.

Using the simulation of wave generation based on the presence of a nonlinear adjustment zone, the difference in amplitude and phase of the wave was determined, as well as its propensity to trigger numerical divergence in the absence of a nonlinear adjustment zone. Therefore, it is necessary to apply a non-linear adjustment zone when simulating wave generation via wave disturbance. In addition, it was determined that desired nonlinear regular and irregular waves can be generated by the transfer function and nonlinear adjustment zone for wave disturbance at a specific point under the boundary condition of the linear free surface.

The simulation of irregular wave generation based on several changes in mean wave steepness was performed and it was determined that the difference in the input spectrum increased as the wave steepness increased. This difference indicates an increase in energy redistribution between wave components owing to nonlinear interactions between wave components. From these results, the HOS method is inferred to be an excellent numerical tool for observing nonlinear interactions between wave components.