Simulation and Modelling for a Three-Dimensional Ocean Sur- face Wave using and Inverse Fourier Transform

Ocean surface waves have been utilized as fundamental information in various fields of oceanic research. In this paper, we suggest a simulation and modelling technique for generating an ocean surface wave using an inverse Fast Fourier Transform (iFFT), and we subsequently verify the wave’s accuracy. The conventional method, linear superposition, requires recursive calculation because of the double summation and the time variable; to circumvent this issue, the new algorithm is presented. The Joint North Sea Wave Project (JONSWAP) spectrum is utilized for the ocean surface wave simulation example, and the parameters are the significant wave height (H ) and the zerocrossing wave period (T ). A coordinate transform for the wavenumber domain was used to apply the inverse FFT algorithm. To verify the accuracy of the simulation result, the relative error between the input condition and the analysis result was calculated. The result for T is below 4% relative error, and the maximum relative error for H is 7%. To avoid the Nyquist frequency for wave-field analysis and simulation, the minimum grid size was calculated by twice applying the maximum wavenumber.


Introduction
Ocean surface waves act as external forces on ships, offshore structures, and other floating bodies in the ocean. To prevent damage to these floating structures, ocean surface waves have been studied for a long time. Analyzing the physical phenomena that arise from ocean surface waves is crucial for investigating the dynamic stability of and structural problems in floating structures. Specifically, for hydro-structure interaction problems resulting from the enlargement of ships, ocean surface waves provide fundamental data for performing time series analysis.
Simulation of ocean surface waves is conducted in many research fields. The rapid development of computer hardware accelerates simulation technologies for ocean surface waves. In previous decades, when computer hardware was insufficient to support such calculations, the ocean's surface was modelled with limited help from statistical analysis [1]. With advancements in computer hardware, the process of modelling ocean waves has evolved to include graphic processing, such as realistic surface rendering [2].
Numerical ocean surface waves have been generated to provide essential data for the fields of naval architecture and ocean engineering. In a previous study [3], an ocean surface wave was used to calculate ship motions in a time series.
In another study [4], an ocean surface wave was utilized to serve as input data for simulating ships and buoys' motion. Moreover, using a wave spectrum, a time series ocean wave was generated to perform quantitative analysis on the structural behavior of ships [5].
In previous studies of ocean surface simulation, the ocean surface simulations have been performed using random numbers based on the Monte Carlo method [6]. In addition, a directional wave spectrum is computed to create a random ocean wave field based on a linear superposition of results from the conventional method, in which sequential computation with a repetitive subroutine is used to simulate the wave surface [7]. Furthermore, numerical ocean surface waves can be applied in x-band marine radar image simulation. The results of numerical ocean surface wave simulations are a fundamental factor that directly affects the accuracy of marine radar image simulation. Representative papers for this application include [8], [9], and [10].
In this study, based on the generation of a directional wave spectrum, the ocean surface wave is simulated, and the results are analyzed. In the efficiency of its calculation time, the proposed method, which utilizes an inverse Fourier transform, has an advantage for large-scale simulations. In particular, the various parameters and coordinate transforms are graphically described to facilitate the simulation procedure in the article. The organization of this paper is as follows. The directional wave spectrum is presented in Section 2. In Section 3, based on the currently used linear superposition method and the inverse Fourier transform method, the simulation of the ocean surface wave is presented [11]. Section 4 contains the results of previous sections. Finally, Section 5 is followed by conclusions.

Ocean Wave Spectrum
Because the wave spectrum acquired by observing long-term sea state information includes the statistical characteristics of observed ocean phenomena, the wave spectrum is useful for providing features of the actual sea state. Although the ocean wave spectrum is comprised of various components, the representative directional wave spectrum is applied in this study, and the directional wave spectrum can be expressed by the following: where ( ) is a wave spectrum that includes frequency characteristics. ( , ) is a directional spreading function given by The directional spreading function directionally spreads the wave spectrum without distorting the energy's magnitude.
A two-parameter JONSWAP spectrum is applied in this study, and the parameters are significant wave height Hs and zero-crossing wave period Tz [12]. The peak enhancement parameter γ is 3.3. Normally, the parameter's value ranges from 1 to 20. For the directional wave spectrum, the half-cosine 2s-power-type spreading function and the Mitsuyasu-type function are used [13].

Wave Simulation using Linear Superposition
To obtain wave elevation ( , ), the linear superposition method is used. Using the Riemann sense, the continuous integration is modified [7], and the equation can be expressed by where = 2 ∆ , ∈ {1, ⋯ , }.
is the function of frequency that contains number of components, and is the amplitude parameter. in Equation (3) is an equivalent random number that ranges from 0 to 2π. Equation (1) is used to compute the three-dimensional irregular wave elevation, and additional directional components are included.
The three-dimensional irregular random wave is , denotes a function of P values in the frequency domain and Q values in the direction domain. Using the above equation, subroutines to compute temporal and spatial dimensions are required. The computation time is proportionally increased by the third power of the increment for each domain. Therefore, the following inverse Fourier transform algorithm is suggested to improve computation time. Figure 1 is data flow diagram on properties and requirements of two-dimensional and three-dimensional irregular ocean surface wave depending on wave spectrum. The figure also shows existing method for ocean surface wave simulation using linear superposition.

Ocean Wave Simulation using Inverse Fourier Transform
To simulate a three-dimensional irregular wave with temporal and spatial dimensions, a set of procedures is required. The set of procedure is characterized by the wavenumber domain setting, the domain transform, and the inverse Fourier transform in the directional spectrum. 'Step 1' involves calculating the wavenumber domain. To establish the domain, the domain is derived from the geometrical information in Cartesian coordinates, which are computationally unique, regardless of the time domain. Equation (18) was proposed to account for the divided form of the terms for space and time. Consequently, the variables for time and space can be considered for the phase of the exponential term and the amplitude term, respectively. 'Step 2' introduces the function value calculation for the wavenumber domain, which belongs to the directional wave spectrum. For the computational procedure, the chain rule and the linear dispersion relation are necessary for applying the coordinate transform. 'Step 3' calculates the temporal information for the wave simulation. This step contains the method for generating a three-dimensional wave surface using an inverse Fourier transform. Figure 2 is directional spreading function with a Mitsuyasu-type spreading function which is frequency dependent function as Smax value. The figure also shows simulation technique using inverse Fourier transform.

Step 1 : Generation of the wavenumber domain ,
To apply the inverse Fourier transform algorithm, the wavenumber domain and the function value must be specified. The desired surface wave is assumed to exist in the Cartesian coordinate grid. The required variables and transform information for generating the surface wave in the grid are shown in Figure 3 and Table 1.   Table 1: along (a) space domain and (b) wave-number domain. The centre of grids in each domain denotes (0,0).
In Figure 3, the center of the generated grid is located at point (0,0), and the size of the grid is uniform. The wavenumber of the domain is generated using equivalent variables.
× is the number of the grid set, which indicates the number of grids in the xdirection and the y-direction. The number of variables in the wavenumber domain is selected to equal the number of grids. ∆ , ∆ , ∆ , and , ∆ indicate the size of the grid, and the variable for each coordinate and each value is constant.
, , , and represent the lengths of the generated area.
Using the variables in Table 1, the vector notation in the space domain and in its corresponding wavenumber domain can be expressed by and ⃗ = ̂+ .
The procedure in Step 1 is independent of time.

3.2.
Step 2 : Coordinate transform of the directional wave spectrum The function value in the generated ⃗ domain is obtained from the directional wave spectrum. The calculated directional wave spectrum is transformed to the function value in the generated ⃗ domain. The function value represents the period and directional energy value and is directly related to the amplitude of the simulated ocean surface wave field. Specifically, the function value has statistical characteristics equivalent to the period of the main ocean surface wave direction. The energy function value is the physical ocean surface elevation. To apply the inverse Fourier transform algorithm, the main point of the coordinate transform in the generated domain from Equation (8) is the energy distribution of the function value, which must be constant. Therefore, the coordinate transformation is conducted using the following integral process: where J is the Jacobian, which is produced using the transformation between polar coordinates and Cartesian coordinates. The transformation between the wavenumber domain and the wave frequency domain can be expressed by is the group velocity that represents the energy propagation velocity, which results from the linear dispersion relation [14].
Equation (13) represents the final transformation between functions of , , and ⃗ . Using the above process, function ⃗ is applied to Step 3.

3.3.
Step 3 : Three-dimensional wave generation using the inverse Fourier transform The inverse Fourier transform commonly used in the engineering field was applied at this step. Specifically, the function ( , ) can be decomposed into the sum of a linear function and its pair, the function ( , ). The relationship between ( , ) and ( , ) is the so-called Fourier transform [15]. The two dimensional Fourier transform is and its pair, the inverse Fourier transform, can be defined as Moreover, based on Table 1, the two-dimensional discrete inverse Fourier transform is defined as where m, n, m′, and n′ are sets of grids with sizes equal to each axis.
To generate the two-dimensional wave surface, the amplitude matrix can be composed with the wavenumber spectrum from Equation (13) to produce the following: The equation includes random-phase temporal information, and Equation (17) is necessary for applying the complex matrix for the inverse Fourier transform. The formulation is as follows: where is the random phase matrix, and its size must equal the dimension of the wavenumber domain. The variable, t, that belongs to the time domain is a fixed number in the subroutine. Moreover, a two-dimensional Fourier shift must be applied to the complex matrix defined by Equation (18). The next equation shows the two-dimensional inverse Fourier transform, which is defined as The real part calculated from Equation (19) can be regarded as the simulated wave field.

Results
This paper focused on reducing the computation time for generating ocean surface waves. Figures 4 through 6 show the results of the generated spectrum, which is essential for simulating the ocean wave surface. The specifications of the workstation are i7-3.4GHz Intel CPU, 16GB random access memory and MATLAB 2018b. Table 2 compares the results for computation time from the suggested algorithm to the linear superposition method. A set of illustrations in Figure 7 shows the three-dimensional ocean wave surface, and Figures 8 through 10 show the results of frequency analysis.    Figure 4 shows the spectral density for various values of gamma in the JONSWAP spectrum. The value of gamma γ has a mean of 3.3 in this paper, while this value typically ranges between 1 and 20. This value means that the peak variation in the spectral density function affects the ocean surface's shape [13]. Moreover, wave conditions consisting of significant wave height and a zero-crossing wave period in the two-parameter-type JON-SWAP spectrum are H = 4m and T = 9.5s.

Ocean Surface wave
A half-cosine 2s-power-type function with a Mituyasu-type subfunction of S was used in the spreading function, as shown in Figure 5, and the main direction of the wave field was 120 degrees. S is the spreading parameter, which was proposed by [16] for engineering application purposes. The values of the wind waves, the swell with a short decay distance, and the swell with a long decay distance are 10, 25, and 75, respectively. The spreading function contains the function of direction and peak wave frequency.  A time-based ocean surface wave can be generated, as shown in Figure 7, and this wave was required to calculate Equation (19) based on Equation (1). The sequential 16 images for the simulated wave surface are shown in Figure 7, and the illustration at the top is the first ocean surface wave field of the series at t = 1. The wave direction is 120 degree incoming, which is defined based on the x-axis, and the time step, ∆t is 1 s for the simulation.

Comparison of computation time for each algorithm
To realize the algorithm, the solver in MATLAB was utilized. The algorithm can be divided into two major types, which are linear superposition and the proposed algorithm, an inverse Fourier transform. Linear superposition can be separated into the fully computational type and the detailed matrix computation type.
In the fully computational case, the double summation in Equation (5) was used in a loop function to perform the calculation. Moreover, to calculate x, y, and t, loop functions were required for simulating variables in the spatial and temporal domains. In total, five loop functions were used, and these five functions represented one main loop with four sub-loop functions. For a greater grid size (space), a larger time domain, and a larger spectral domain, the increased computation time requirements can be augmented in proportion to the size of the domain. The simulation occupied about 833 s, with 128 by 128 pixels for one image. If the image is increased by a factor of four, computation time approximately quadrupled over the case of a single image, as presented in Table 2 For matrix computation, instead of using a double summation for P and Q to compute Equation (5) in the fully computational case, the double-loop function can be calculated in a matrix. For the matrix computation, only the main loop of the time domain and two sub-loops of the space domain are used. As shown in Table 2the time required increases in proportion to the grid size as the spatial domain and the temporal domain become large.
Using the inverse Fourier transform from Equation (19) (which represents the proposed algorithm), the computation time was reduced to a minimum of 1/100 or less than the previous two algorithms in Table 2 All results in Table 2 show the computation time with ten repetitions and include the average time.

Verification of the simulation
To verify the simulation, three different verification methods were attempted based on wave direction, wave period, and wave height. These representative parameters are typically used as input conditions and for result verification. Furthermore, the simulation was considered with both a specific wave condition and also a combination of stochastic possibilities. The relative error that was proposed applies the zero-crossing wave period and the significant wave height as the input and the analyzed result. For the first step, spectral analysis of the A-A' section was performed, and this section represents the main wave direction for the simulated ocean surface wave shown in Figure  8. The thick blue line and the solid green line represent the main direction of the input wave spectrum and the evaluated result, respectively. According to Figure 8, the result indicates that the peaks are reasonably matched for both the input wave spectrum and the analyzed spectrum. As shown in Figure 8, compared with the blue line, the energy distribution of the green line along the wave frequency axis is more widely distributed. In this case, the integration areas were obtained along with the wave frequency, which agreed well with both the input spectrum and the analyzed result. This finding is necessarily limited because the directional spreading function was applied to the wave spectrum. Spectral analysis result of its surface wave field. The radial direction is wave frequency, and the angular direction is the incoming wave direction.  For the next step, a two-dimensional Fourier transform was used for the directionality and period of the ocean surface wave. Figure 9(a) is the circular input directional wave spectrum, and Figure 9(b) is the result of the two-dimensional Fourier transform of the ocean surface wave. The comparative result in Figure 9 shows that the location peak of the analyzed result is equivalent to the directional wave spectrum, even if this result is slightly distorted because of the spreading effect. Figure 10 shows the analyzed image of a two-dimensional surface wave with several combinations of wave directions and wave periods as the input conditions. Table 4 presents the results, including the main direction and the zero-crossing wave period, and this table compares these results with the input conditions. To perform the analysis, the zerocrossing wave period from the analysis result must be calculated to the zeroth-order spec-tral moment, which must then be divided by the second-order spectral moment. The directionality results shown in Table 3 were obtained from measuring the peak point using the expression tan K /K . Compared with the inputs, the measured values are almost in agreement. Moreover, for the two-dimensional Fourier transform, the spectral density in each case is asymmetric around the origin because of the directional ambiguity resulting from one ocean surface wave image.

Wave height
The statistical analysis for wave height (performed with ten repetitions) is shown in Figure 11. The wave conditions are H = 4m and T = 9.5s, and the spectrum parameters are The wave height Rayleigh distribution with wave conditions H = 4m and T = 9.5s is shown in Figure 10, and based on the distribution, H / was calculated to be 3.742 m.

Relative error of the simulation
To validate the accuracy of the simulation, the relative error between the input condition and the analysis result was calculated. Significant wave height H and zero-crossing wave period T are the main input conditions of the simulation. H and T were calculated using spectral analysis from the generated ocean wave surface.
The relative error for the zero-crossing wave period at different significant wave heights is shown in Figure 12. It is notable that all the results are below 4% relative error for T . In Figure 13, the relative error for the significant wave height in a different zerocrossing wave period is presented. When T is less than 16 s, the relative error is about 2%. The maximum relative error for HS when T is larger than 16 s is 7%.

Avoiding the Nyquist frequency
When the sampling frequency is less than twice of the peak wave frequency T , aliasing can occur because of folding. Half of the peak frequency is called the folding frequency or the Nyquist frequency. To avoid this issue for wave field analysis and simulation, the minimum grid size was calculated by applying twice the maximum wavenumber using the dispersion relation is the peak wave period, which represents the maximum value of the two-dimensional wave field from the FFT image in Figure 9, and is gravitational acceleration. In the deep-water condition, the term tanh ℎ can be approximated as 1, and Equation (20) can be simplified by Figure 14. Example of Nyquist frequency with T = 9.3s. Figure 14 shows the analysis results for the mean T values after fifty calculations with T = 9.3s. The peak wavelength, λ , which is indicated with the blue vertical solid line in Figure 14, is about 134.8 m based on Equation (22). The red vertical solid line in Figure 14 is the minimum value of the Nyquist frequency, which represents half of the peak wavelength. The x-axis and y-axis are the grid size and the peak wave period, respectively. The width and height of the grids are equal. The black solid line, the red circle, and the blue cross markers represent the peak wave period for the input condition, the peak wave period for the wave spectrum, and the peak wave period for the two-dimensional wave field generated using FFT analysis, respectively.
In figure 15, the resolution of ocean wave field is increased as the grid size is decreased. Moreover, distortion can occur due to the relationship between grid size and peak wave period in simulation. That is, the grid size should be smaller than the peak wave period, and a result without distortion can be obtained as shown in Figure 5(a). As the grid size is close to the Nyquist frequency, the ocean surface image becomes distorted Figure 15(b). Or, if the grid size is larger, you will get unintended results as shown in Figure 15(c).

Conclusion
Using the inverse Fourier transform method, the proposed algorithm has an advantage over typical linear superposition methods because of its computation time efficiency. The simulated results produced matching statistical characteristics.
Generating a three-dimensional ocean wave surface in space and time based on the linear superposition method increased the computation time because of the sub-loop functions for space, time, and double summation. By using the proposed algorithm, the computation time was substantially reduced because the proposed method only computed the subroutine functions in the time domain, and the other spatial information was pre-calculated using the inverse Fourier transform.
The generated three-dimensional ocean wave surface was verified using three different techniques. First, for wave frequency, the output frequency was compared with the input frequency using spectral analysis and a two-dimensional Fourier transform. Based on the comparison result for the input and output frequencies, the peaks of the frequencies' spectral density were substantially matched. Second, wave directionality was analyzed using a two-dimensional Fourier transform, and the angular components were considerably matched. Finally, the statistical characteristics of wave height were compared with the Rayleigh distribution in a ten-time repeatability test. The relative error between input and output for significant wave height and the zero-crossing wave period was calculated to validate the simulation results. The relative error was below 4% in TZ. In HS the maximum relative error was 7%. The error be-tween input and output significant wave heights was 5%. The input peak wave period and the grid size were optimized to avoid the Nyquist frequency using a dispersion relation. The minimum grid size should be smaller than half of the maximum wavelength. Data Availability Statement: All data presented in this study are available.