The fluid region is described by using the 3D Cartesian coordinates (x, y, z), with x the longitudinal coordinate, y the transverse coordinate, and z the vertical coordinate (positive upwards). Time is denoted by t. The location of the free surface is at z = η(x, y, t) at a specific time of t.

For real fluid flows in the sea surface, it is reasonable to neglect the effects of viscosity. Meanwhile, sea water has a low compressibility, and therefore normally sea water flows are incompressible. Considering the aforementioned facts, it is reasonably accurate to assume that the fluids on the sea surface are ideal (i.e., incompressible and inviscid). Furthermore, rotation of a fluid particle can be caused only by a torque applied by shear forces on the sides of the particle. Since shear forces are absent in an ideal fluid, the flow of ideal fluids is essentially irrotational. Then, the velocity potential $\varPhi \left( {x,y,z,t} \right)$ exists. If the water depth d at the sea bottom is constant, then for constant water depth d, the velocity potential $\varPhi \left( {x,y,z,t} \right)$ and the free surface elevation η(x, y, t) can be determined by solving the following boundary value problem:

In this study the following expansion is utilized to solve the system (1)–(4):

where ε is a small parameter which is typically proportional to the wave steepness. For an irregular sea state characterized by a specific wave spectrum ${S_{\eta \eta }}\left( \omega \right)$ in which $\omega $ denotes the angular frequency, it can be shown that a first order linear solution of the system (1)–(4) can be expressed as follows:

where $N$ is usually chosen to be a sufficiently large positive integer, i stands for the imaginary unit, ${c_n}$ denotes the sinusoidal wave component amplitude, ${\varepsilon _n}$ is the phase angle uniformly distributed in the interval [0, 2${\text π} $], ${\omega _n}$ denotes the radian frequency (${\omega _n} = 2{\text π} n/T$, T is the time interval) and ${k_n}$ denotes the wave number. ${\omega _n}$ and ${k_n}$ are related through the following linear dispersion relation:

where $d$ and $g$ are respectively the water depth and the gravitational acceleration.

Equation (7) describes an ideal linear irregular sea model and this model has been derived under the assumption that the wave heights are small compared to the wave length. Then, the surface elevation resulting from irregular ocean waves can be approximated as a superposition of multiple harmonic waves with different amplitudes and phases. Furthermore, in this ideal linear irregular sea model the random phase angles, ${\varepsilon _n}$ are assumed to be uniformly distributed between 0 and 2π. However, for modelling shallow water nonlinear waves, the above linear irregular sea model should be corrected by including second order terms as follows:

The terms $P({\omega _n},{\omega _m})$, ${r_{mn}}$ and ${q_{mn}}$ in Eqs (9) and (10) are called quadratic transfer functions, and there expressions are given as follows:

where the Kroenecker delta (${\text{δ} _{ - n,m}}$=1, if n+m = 0, 0 otherwise) is introduced to avoid a singular $P({\omega _n},{\omega _m})$. The wave surface elevations for the second order nonlinear waves can finally be obtained by combining Eq. (7) and Eq. (10) as follows:

The proposed second order random wave simulation method starts with taking a multi-directional spectrum $S\left( {\omega ,\theta } \right)$ and integrating the energy over all directions to give the total energy at each frequency. The obtained equivalent frequency spectrum $S\left( \omega \right)$ is then utilized to generate a nonlinear wave elevation time series by applying Eqs (8) and (14) at a specific site of the sea. Finally it should be pointed out that the coefficients ${c_n}$(or ${c_m}$) in Eq. (14) are determined by using a given spectrum $S\left( \omega \right)$ as follows:

where $\Delta \omega = {\omega _C}/N$, and ${\omega _C}$ is the upper cut off frequency beyond which the wave spectrum $S\left( \omega \right)$ may be assumed to be 0 for either mathematical or physical reasons.

As we know, a wave spectrum is actually obtained from the observed time series of water elevation. The transformation from water elevation to wave spectrum is based on the assumption that irregular waves can be treated as a combination of linear waves with different amplitudes, different frequencies and different phases. It would be better to give an explanation for how to reasonably consider the nonlinearity of waves in the inverse transformation (i.e., the transformation from the wave spectrum to the nonlinear water elevation). In the following, we will raise a specific calculation example of this inverse transformation.

This study shows as an example to use Eq. (7) to simulate the linear Gaussian wave time histories of a sea state with a JONSWAP spectrum with a significant wave height ${H_S}$= 6 m, a spectral peak period ${T_p}$=8 s, a peakedness factor $\gamma$=1 and corresponding to a water depth of 10 m.

Figure 1 shows our simulation of the linear part wave elevation time series ${\eta ^{\left( 1 \right)}}(x,t)$ which contains 150 wave elevation points (with an equal time distant of 0.1 s between two successive points). A wave elevation point has two coordinates (the first coordinate is time, and the second coordinate is the value of the wave elevation above the mean water level). The computation is based on Eq. (7) and the aforementioned JONSWAP spectrum.

Figure 2 shows our simulation of the second order correction part wave elevation time series ${\eta ^{\left( 2 \right)}}(x,t)$ which contains 150 wave elevation points (with an equal time distant of 0.1 s between two successive points). The computation is based on Eq. (10) and the aforementioned JONSWAP spectrum. This is precisely an example for how to reasonably consider the nonlinear correction part of waves in the inverse transformation (i.e., the transformation from the wave spectrum to the nonlinear correction part water elevation).

In Fig. 3 the blue curve shows our simulation of the entire $\eta (x,t)$ (the linear part ${\eta ^{\left( 1 \right)}}(x,t)$ plus the second order correction part ${\eta ^{\left( 2 \right)}}(x,t)$) wave elevation time series which contains 150 wave elevation points (with an equal time distant of 0.1 s between two successive points). The computation is based on Eq. (14) and the aforementioned JONSWAP spectrum. In Fig. 3, the red curve shows our simulation of the linear part wave elevation time series ${\eta ^{\left( 1 \right)}}(x,t)$ which contains 150 wave elevation points based on the aforementioned JONSWAP spectrum. In Fig. 3, the green curve shows our simulation of the second order correction part wave elevation time series ${\eta ^{\left( 2 \right)}}(x,t)$ which contains 150 wave elevation points based on the aforementioned JONSWAP spectrum. We can see that the wave crests of the entire nonlinear waves $\eta (x,t)$ have become steeper and higher. This is a precise example for how to reasonably consider the nonlinearity of the entire waves in the inverse transformation (i.e., the transformation from the wave spectrum to the entire nonlinear waves).

In this sub-section, the accuracy of the proposed second order random wave simulation method will be validated by a calculation example. Specifically, the proposed second order random wave simulation method will be applied for calculating the wave crest amplitude exceedance probabilities of a sea state with a multi-directional wave spectrum based on the measured surface elevation data at the coast of Yura. The measured surface elevation data at the coast of Yura were obtained at a location 3 km off Yura fishing harbor facing the Sea of Japan. The observations were carried out during the period from 11:10 am to 2:08 pm on November 24, 1987 by the Ship Research Institute, Ministry of Transport of Japan. Temporal sea surface elevations were measured with ultrasonic-type wave gages installed at three points in 42 m water depth. The sampling time interval during the measurement was 1 s. Based on these measured Yura coast surface elevation data, a multi-directional wave spectrum was estimated by using the Maximum Likelihood Method and is shown in Fig. 4. Fig. 5 shows a part of the measured wave elevation time series by the mid-point wave gage at this site.

In the field of ocean engineering, based on a specific wave spectrum, an empirical (or theoretical) model or numerical simulation methods can be used to calculate the exceedance probabilities of the wave crest amplitudes (which is directly related to the occurrence probability of the extreme waves). In the following, taking the multi-directional wave spectrum in Fig. 4 as a calculation example, this study will test the proposed second order random wave simulation method to calculate the exceedance probabilities of the wave crest amplitudes, and compare its accuracy with those of a theoretical model or a linear simulation method.

Figure 6 shows the calculated wave crest amplitudes exceedance probabilities based on the multi-directional spectrum shown in Fig. 4 and the measured Yura coast wave data. The continuous green curve in Fig. 6 represents the calculation results of the wave crest amplitudes exceedance probabilities directly obtained from the measured Yura coast wave data that contains 10700 wave elevation points. This continuous green curves based on the measured Yura coast wave data are used as the benchmark against the accuracy of the results from the various numerical simulation methods and from an existing theoretical wave crest amplitudes model which is checked.

The continuous red curve in Fig. 6 represents the results of the wave crest amplitudes exceedance probabilities obtained from using the theoretical Rayleigh distribution model expressed as follows:

where h represents the wave crest amplitude, ${H_{\rm{{S}}}}$ represents the significant wave height, A_c represents crest amplitude. The wave crest amplitudes exceedance probabilities obtained from using the theoretical Rayleigh distribution model deviate a lot from the corresponding benchmark results obtained directly from the measured Yura coast wave data. This is not a surprise because the theoretical Rayleigh distribution is an ideal linear model and the measured Yura coast wave elevation series are obviously nonlinear. The continuous black curve in Fig. 6 represents the results of the wave crest amplitude exceedance probabilities obtained from using the linear simulation method based on the multi-directional spectrum shown in Fig. 4. The linear simulation process started with taking the multi-directional spectrum shown in Fig. 4 and integrating the energy over all directions to give the total energy at each frequency. The obtained equivalent frequency spectrum was then used to generate a wave elevation time series of 800000 points by applying Eqs (7) and (8). Next, the wave crest amplitudes time series were extracted from these 800000 wave elevation points. Finally, the wave crests amplitudes exceedance probabilities were obtained by statistical and mathematical processing the extracted wave crest amplitudes time series. Obviously, the wave crest amplitudes exceedance probabilities obtained from using the linear simulation method also deviate substantially from the corresponding benchmark results obtained directly from the measured Yura coast wave data.

In order to more accurately calculate the wave crest amplitudes exceedance probabilities, the proposed second order random wave simulation method was tried to use. In Fig. 6, a pink * represents the result of the wave crest amplitudes exceedance probability obtained from using the proposed second order random wave simulation method based on the multi-directional spectrum shown in Fig. 4. The proposed second order random wave simulation process started with taking the multi-directional spectrum shown in Fig. 4 and integrating the energy over all directions to give the total energy at each frequency. The obtained equivalent frequency spectrum was then utilized to generate a nonlinear wave elevation time series of 800000 points by applying Eqs (8) and (14). Next, the wave crest amplitudes time series were extracted from these 800000 wave elevation points. Finally, the wave crests amplitudes exceedance probabilities were obtained by statistical and mathematical processing the extracted wave crest amplitudes time series. Note that in the figure legend the phrase “from nonlinear simulation” exactly means “from second order random wave simulation”. The wave crest amplitudes exceedance probabilities obtained from using the proposed second order random wave simulation method fit quite well with the corresponding benchmark results obtained directly from the measured Yura coast wave data. The accuracy of the proposed second order random wave simulation method is therefore convincingly validated.

Our specific calculation examples will be carried out regarding a heaving two-body point absorber wave energy converter. Figure 7 shows this specific wave energy converter modeled in WEC-Sim (http://wec-sim.github.io/WEC-Sim/), an open source computer software for simulating WEC performances. The dimensions of this specific wave energy converter are shown in Fig. 8. This two-body point absorber WEC consists of a float and a spar/plate. The float has a diameter of 20 m, a thickness of 5 m and a mass of 727.01 t. The spar has a height of 38 m and a diameter of 6 m. The spar/plate has a mass of 878.3 t.

It should be emphasized that this two-body point absorber WEC is installed in a coastal sea area with a water depth of only 49.5 m. Obviously, the random waves in this shallow water area will be nonlinear according to the nonlinear random wave theory in Ochi (1998).

In most wave energy exploitation engineering projects that met in the real world, the only information knew beforehand regarding a sea state usually will be a specific wave spectrum. In Section 2, the accuracy of the proposed second order random wave simulation method was already demonstrated in generating irregular waves based on a measured multi-directional wave spectrum. In the following calculation examples, the proposed second order random wave simulation method was utilized for generating nonlinear irregular waves for calculating the wave excitation load ${{\boldsymbol{P}}_{{\rm{ext}}}}\left( t \right)$ in Eq. (17). The absorbed power of the aforementioned WEC Power-Take-Off system was predicted for some specific sea states characterized with multi-directional JONSWAP spectra. Figure 9 shows one of these multi-directional JONSWAP spectra with a significant wave height ${H_{\rm{{S}}}}$=1 m, a spectral peak period ${T_p}$=12 s, a peakedness factor $\gamma $=1, and a cosine squared spreading function with the spreading parameter equal to 15. The mathematical expression of this multi-directional JONSWAP spectrum is as follows:

where $S\left( \omega \right)$ is a frequency spectrum and $D\left( \theta \right)$ is a spreading function.

The mathematical expression for the frequency spectrum $S\left( \omega \right)$ is as follows (Wang, 2014):

where $\text{δ}= 0.07,\;{\rm{if}}\; \;\omega < 2{\pi} /{T_p};$ $\text{δ} = 0.09, \; {\rm{if}}\; \;\omega \geqslant 2{\pi} /{T_p}.$ $\omega $ is the wave angular frequency.

The mathematical expression for the spreading function $D\left( \theta \right)$ in Eq. (18) is as follows:

where the spreading parameter s=15.

Figure 10 is a corresponding polar plot of the multi-directional JONSWAP spectrum as shown in Fig. 9.

Next, the calculation results of the absorbed power of the aforementioned wave energy converter were placed in an ideal linear sea versus in a multi-directional nonlinear random sea. Figure 11 shows the predicted WEC absorbed power time series under the sea state of linear irregular waves based on a multi-directional JONSWAP wave spectrum with ${H_{\rm{{S}}}}$=12 m, ${T_p}$=12 s, $\gamma $=1 and a cosine squared spreading function with the spreading parameter equal to 15. Figure 12 shows the predicted WEC absorbed power time series under the sea state of nonlinear irregular waves based on the same multi-directional JONSWAP wave spectrum with ${H_{\rm{{S}}}}$=12 m, ${T_p}$=12 s, $\gamma $=1 and a cosine squared spreading function with the spreading parameter equal to 15. The calculation results in Figs 11 and 12 were obtained by solving the WEC nonlinear dynamic filter Eq. (17) in WEC-Sim. However, WEC-Sim does not have built-in functions for generating nonlinear irregular waves for calculating the wave excitation load ${{\boldsymbol{P}}_{{\rm{ext}}}}\left( t \right)$ in Eq. (17). Therefore, this study externally generated nonlinear irregular waves by using the proposed second order random wave simulation method which imported into WEC-Sim. The nonlinear irregular wave simulation was started by taking the multi-directional JONSWAP wave spectrum (with ${H_{\rm{{S}}}}$=12 m, ${T_p}$=12 s, $\gamma $=1) and integrating the energy over all directions to give the total energy at each frequency. The obtained equivalent frequency spectrum was then utilized to generate a nonlinear wave elevation time series of 1200 points by applying Eqs (8) and (14). The generated nonlinear irregular waves time series were saved as a.mat file and then imported into WEC-Sim for calculating the wave excitation load ${{\boldsymbol{P}}_{{\rm{ext}}}}\left( t \right)$ in Eq. (17) and subsequently obtaining the WEC absorbed power time series in Fig. 12. For comparison purpose, in this study, externally generated linear irregular waves by using Eqs (7) and (8) were also imported into WEC-Sim. The linear irregular wave simulation was started by taking the multi-directional JONSWAP wave spectrum (with ${H_{\rm{S}}}$=12 m, ${T_p}$=12 s, $\gamma $=1) and integrating the energy over all directions to give the total energy at each frequency. The obtained equivalent frequency spectrum was then utilized to generate a linear wave elevation time series of 1200 points by applying Eqs (7) and (8). The generated linear irregular waves time series were saved as a .mat file and then imported into WEC-Sim for calculating the wave excitation load ${{\boldsymbol{P}}_{{\rm{ext}}}}\left( t \right)$ in Eq. (17) and subsequently obtaining the WEC absorbed power time series in Fig. 11. Having obtained the prediction results as shown in Figs 11 and 12, the statistical characteristic values were subsequently calculated based on the time series shown in these two figures. Our calculation results are summarized in the last row in Table 1.

By carefully comparing and analyzing the calculation results, the mean value of the 1 200 s WEC absorbed power under the sea state of ideal linear irregular waves is smaller than the corresponding mean absorbed power value when inputting nonlinear irregular waves. In order to investigate the influences of choosing different significant wave height values on the power performances of wave energy converters, the aforementioned WEC absorbed power values were sequently calculated under 11 other sea states characterized with a multi-directional JONSWAP wave spectrum with the following parameters respectively: ${H_S}$=1 m, ${T_{\rm{p}}}$=12 s, $\gamma $=1, s=15; ${H_{\rm{S}}}$=2 m, ${T_p}$=12 s, $\gamma $=1, s=15; ${H_{\rm{S}}}$=3 m, ${T_p}$=12 s, $\gamma $=1, s=15; ${H_{\rm{S}}}$=4 m, ${T_p}$=12 s, $\gamma $=1, s=15; ${H_S}$=5 m, ${T_p}$=12 s, $\gamma $=1, s=15; ${H_{\rm{S}}}$=6 m, ${T_p}$=12 s, $\gamma $=1, s=15; ${H_{\rm{S}}}$=7 m, ${T_p}$=12 s, $\gamma $=1, s=15; ${H_{\rm{S}}}$=8 m, ${T_p}$=12 s, $\gamma $=1, s=15; ${H_{\rm{S}}}$=9 m, ${T_p}$=12 s, $\gamma $=1, s=15; ${H_{\rm{S}}}$=10 m, ${T_p}$=12 s, $\gamma $=1, s=15; ${H_{\rm{S}}}$=11 m, ${T_p}$=12 s, $\gamma $=1, s=15. The calculation results were also summarized in Table 1. By carefully comparing and analyzing these calculation results, in all these 11 sea states, the mean value of the 1 200 s WEC absorbed power under the sea state of ideal linear irregular waves is always smaller than the corresponding mean absorbed power value when inputting nonlinear irregular waves.

In order to further investigate the influences of using different random seed numbers for generating irregular waves on the power performances of wave energy converters, the aforementioned WEC absorbed power values were sequently calculated when inputting 11 different irregular wave elevation time series simulated based on the same multi-directional JONSWAP wave spectrum (${H_{\rm{S}}}$=12 m, ${T_p}$=12 s, $\gamma $=1, s=15) using different random seed numbers. The calculation results were summarized in Table 2. By carefully comparing and analyzing these calculation results, in all these 12 sea states, the mean value of the 1 200 s WEC absorbed power under the sea state of ideal linear irregular waves is always smaller than the corresponding mean absorbed power value when inputting nonlinear irregular waves. In one specific case (No. 11) the predicted WEC absorbed power value with the nonlinear irregular waves as inputs is 15.8% larger than that predicted with the linear irregular waves as inputs. All the aforementioned calculation results highlight the vital importance of using the nonlinear irregular waves simulated based on a multi-directional wave spectrum when studying the power performances of wave energy converters.