1 Introduction

In recent years, much research attention has been devoted to wave-structure interactions in open seas and nonlinear waves. The numerical research on wave-structure interaction take interest in large-scale model and real-sea state simulations. This gives a big challenge on modeling the nonlinear wave propagation in large domain and the application of computational fluid dynamics (CFD). With the development of computer science, the CFD method becomes a popular approach due to its ability to resolve the complicated wave-structure interaction problems such as impact pressure on structures and wave patterns around high-speed ship. The ability of the CFD method to deal with large-scale model and refinement fluid simulations makes it become the main method in the future. However, the cost of propagating open-state waves or freak waves in the CFD method cannot be ignored. As is known, the viscous effect in incident wave is small, therefore employing the potential theory to do the fast wave generation is a good choice in numerical simulations. The traditional potential theories are mature in generating nonlinear waves, but a spectral method called high-order-spectral (HOS) method has salient ability on dealing with long-time wave evolution and large domain propagation in an efficient way. Yet its shortcomings are obvious: it has difficulties in simulating large amplitude motion of structures and complicated wave phenomena such as wave breaking. Therefore, combining the HOS method with the CFD method seems to be the best way to solve complex wave-structure interactions efficiently and accurately. This paper reviews the development of the HOS method and includes the coupling method in recent years. Besides, based on the strength as well as shortcomings of the combined method, this paper discusses the future research work of this combined method.

2 HOS Method

The HOS method is a robust numerical method for solving nonlinear wave equations. With the application of pseudo-spectral method and Fast Fourier Transform, the HOS method embodies the feature of efficiency and fast convergence. This method has another name called Dirichlet-Neumann Operator^[1]. Its accelerated version^[2] is equivalent to the HOS method^[3]. Ducrozet et al.^[4] gave a comparison and conclusion that HOS models are more efficient than advanced models with finite-difference discretization. Dommermuth and Yue^[5] compared spectral method with fully nonlinear mixed-Eulerian-Lagrangian (MEL)^[6], and this spectral method accomplishes a more rapid convergence in calculation as well as total energy conservation. With the computational effort increasing only linearly with modes and orders, the spectral method costs less computational time.

2.1 Theory of the HOS Method

Dommermuth and Yue^[5] gave the numerical description of the HOS method. They extended Zakharov equation^[7-8]/model coupling idea^[9-11] to steeper waves, and modified this equation and the approach to a higher order for gravity waves. It can simulate the wave nonlinearity into a specified order and use a large number of free Fourier modes in nonlinear simulations. The computational effort increases only proportional to the order and modes. The fluid is considered to be irrotational, homogeneous, incompressible, and inviscid. They considered a velocity potential which satisfies Laplace's equation. With the pseudo-spectral method, the equation only accounts for the free surface, thus they defined the surface potential as ϕ^s(x, t) = ϕ(x, η (x, t), t). The kinematic and dynamic boundary conditions on the free surface of ϕ^s are^[5]

$ {{\eta _t} + {\nabla _x}{\phi ^s} \cdot {\nabla _x}\eta - (1 + {\nabla _x}\eta \cdot {\nabla _x}\eta ){\phi _z}(\mathit{\boldsymbol{x}},\eta ,t) = 0} $

(1)

$ \begin{array}{l} \phi _t^s + \eta + \frac{1}{2}{\nabla _x}{\phi ^s} \cdot {\nabla _x}{\phi ^s} - \frac{1}{2}(1 + \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\nabla _x}\eta \cdot {\nabla _x}\eta )\phi _z^2(\mathit{\boldsymbol{x}},\eta ,t) = - {P_a} \end{array} $

(2)

They assumed a measure of the wave steepness ε as a small parameter. ϕ and η are O(ε) quantities. For a given order M in ε, ϕ can be written in a perturbation series based on ε^[5]. Thus ϕ^(m) is the same quantity with O(ε^m). They expanded each ϕ^(m) on z = η in a Taylor series, then^[5]

$ \begin{array}{*{20}{l}} {{\phi ^s}(x,t) = \phi (\mathit{\boldsymbol{x}},\eta ,t) = }\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \sum\limits_{m - 1}^M {\sum\limits_{k = 0}^{M - m} {\frac{{{\eta ^k}}}{{k!}}} } \frac{{{\partial ^k}}}{{\partial {z^k}}}{\phi ^{(m)}}(\mathit{\boldsymbol{x}},0,t)} \end{array} $

(3)

When expanding Eq. (3)^[12-13], they got a sequence of equations of ϕ^(m):

$ \begin{array}{l} {\phi ^s}(\mathit{\boldsymbol{x}},t) = {\phi ^{(1)}}(\mathit{\boldsymbol{x}},0,t) + \eta \frac{{\partial {\phi ^{(1)}}}}{{\partial z}}(\mathit{\boldsymbol{x}},0,t) + \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \frac{{{\eta ^2}}}{2}\frac{{{\partial ^2}{\phi ^{(1)}}}}{{\partial {z^2}}}(\mathit{\boldsymbol{x}},0,t) + \ldots + {\phi ^{(2)}}(\mathit{\boldsymbol{x}},0,t) + \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \eta \frac{{\partial {\phi ^{(2)}}}}{{\partial z}}(\mathit{\boldsymbol{x}},0,t) + \frac{{{\eta ^2}}}{2}\frac{{{\partial ^2}{\phi ^{(2)}}}}{{\partial {z^2}}}(\mathit{\boldsymbol{x}},0,t) + \ldots + \ldots \end{array} $

Collecting η in the same order, Eq. (3) is a Dirichlet boundary condition. The surface potential ϕ^s(x, 0) and η(x, 0) are known in

$ \begin{array}{l} {\phi ^{(m)}}(\mathit{\boldsymbol{x}},0,t) = - {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \sum\limits_{k = 1}^{m - 1} {{\kern 1pt} {\kern 1pt} {\kern 1pt} \frac{{{\eta ^k}}}{{k!}}} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \frac{{{\partial ^k}}}{{\partial {z^k}}}{\phi ^{(m - k)}}(\mathit{\boldsymbol{x}},0,t),\\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} m = 2,3, \ldots ,M \end{array} $

(4)

where m = 1, ϕ^s = ϕ^(m)(x, 0, t). Eq. (4) is the final formulation of the HOS method. While in a mode-coupling approach, except Dirichlet free-surface condition (4), Eq. (3) can be solved by using eigenfunction^[5].

The eigenfunction φ_n has analytical solutions for simple geometries. For example, the 2π-periodic boundary conditions can be represented as^[5]

$ \begin{array}{*{20}{l}} {{\phi ^{(m)}}(\mathit{\boldsymbol{x}},z,t) = }\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \sum\limits_{n = 0}^\infty {\phi _n^{(m)}} (t){\rm{exp}}[|{\mathit{\boldsymbol{k}}_\mathit{\boldsymbol{n}}}|z + {\rm{i}}{\kern 1pt} {\mathit{\boldsymbol{k}}_\mathit{\boldsymbol{n}}} \cdot \mathit{\boldsymbol{x}}]} \end{array} $

For constant finite deep water, with depth h^[5],

$ \begin{array}{*{20}{l}} {{\phi ^{(m)}}(\mathit{\boldsymbol{x}},z,t) = }\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \sum\limits_{n = 0}^\infty {\phi _n^{(m)}} (t)\frac{{{\rm{cosh}}[|{\mathit{\boldsymbol{k}}_\mathit{\boldsymbol{n}}}|(z + h)]}}{{{\rm{cosh}}(|{\mathit{\boldsymbol{k}}_\mathit{\boldsymbol{n}}}|h)}}{\rm{exp}}({\rm{i}}{\kern 1pt} {\mathit{\boldsymbol{k}}_\mathit{\boldsymbol{n}}} \cdot \mathit{\boldsymbol{x}})} \end{array} $

The study of Dommermuth and Yue^[5] is essential for developing the the HOS method. With its high efficiency and accuracy, it opens up a path in simulating nonlinear gravity waves for succeeding researchers. With the existence of the exponential profile, the HOS method can provide rapid convergence and consume less time. The boundary-value problems are solved rapidly through Fast Fourier Transform, which connects spatial derivatives in spectral representation and nonlinear products in physical space.

West et al.^[14] presented using the HOS method to study the free and bound waves on nonlinear ocean surface. The descriptions of the HOS method is similar to Dommermuth and Yue^[5]. However, there is a little difference in dealing with the vertical velocity in West^[14].They first assumed a power series in ε:

$ {{\phi _z}(\mathit{\boldsymbol{x}},\eta ,t) = \sum\limits_{m = 1}^M {\varphi _z^{(m)}} (\mathit{\boldsymbol{x}},0,t)} $

(5)

$ {\phi _z^{(m)}(\mathit{\boldsymbol{x}},\eta ,t) = \sum\limits_{k = 0}^{m - 1} {\frac{{{\eta ^k}}}{{k!}}} \frac{{{\partial ^{k + 1}}}}{{\partial {z^{k + 1}}}}{\varphi ^{m - k}}(\mathit{\boldsymbol{x}},0,t)} $

(6)

In this way, all the terms under boundary condition which contain the vertical velocity are consistent with the ordering ε, making the comparison with analytical data easier.

Although the advantages of the HOS method are obvious, some errors exist due to its features. Because the velocity potential is single value in the HOS method, it cannot simulate overturn phenomenon in waves. Besides, the HOS method has a limit in wave steepness, and it will fail when wave steepness is very large.

The first error is the truncation in the number of modes N and order M. For the singularity in the analytic continuation of velocity potential in the HOS method, the convergence breaks when the wave steepness is large. The second is amplification of round-off errors. The error at any order in high wavenumber modes will be amplified. The unstable growth of this error in many nonlinear free-surface simulations will occur when high wave number is used^[15-16]. Dommermuth and Yue^[5] used Longuet-Higgins and Cokelet's methods to make a low-pass filter to eliminate such high-wavenumber instabilities.

The aliasing errors^[17-18] happen because the nonlinear products in free-surface boundary conditions are performed in physical space instead of in spectral space. West et al.^[14] removed this error by changing the numbers of mesh points N_x, N_y in the physical plane to satisfy the condition

$ {N_x} > (M + 1){k_{{\rm{max}}}},{N_y} > (M + 1){l_{{\rm{max}}}} $

De-aliased computation can be performed by extending spectra in zero padding^[19-20]. The zero-padding applies half-rule and are defined on a number of points:

$ {N_d} = \frac{{p + 1}}{2}M $

2.2 Development and Application of the HOS Method

Liu et al.^[21] extended the HOS method to study nonlinear wave-body interactions. They considered the submerged body as a dipole distributed on the free surface and a source affected on the body surface. They showed the ability and accuracy of the spectral method in simulating submerged body. However, they also pointed out the method is not available when the body submergence is small and in wave breaking case.

Unlike many numerical models^[8,10,22] which have difficulties in simulating deep-water dispersive waves due to the complexity of the nonlinear terms, the HOS method is more suitable for strong nonlinearity terms. However, it has difficulty in simulating the progress of nonlinear wave which may give spurious high-frequency waves. It can be eliminated by either giving sufficient time when initialized with linear wave or changing initialization procedures. Dommermuth^[23] developed an adjustment for simulating the turbulent free-surface flows. In case of turbulent free surface flows, Dommermuth^[24] also developed an adjustment that allows nonlinear free-surface simulations to be initialized with linear solutions. By introducing the vertical component velocity ϕ_{z|z = η} = W^s(x, t), and letting the first order vertical velocity denote the leading-order, the free-surface Eqs. (1)-(2) can be expressed as^[24]

$ \begin{array}{*{20}{c}} {{\eta _t} - {W^{(1)}} = - {W^{(1)}} - {\nabla _x}{\phi ^s} \cdot {\nabla _x}\eta + }\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} (1 + {\nabla _x}\eta \cdot {\nabla _x}\eta ){W^s}}\\ {\phi _t^s + \eta = - \frac{1}{2}{\nabla _x}{\phi ^s} \cdot {\nabla _x}{\phi ^s} + \frac{1}{2}(1 + }\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\nabla _x}\eta \cdot {\nabla _x}\eta ){{({W^s})}^2}} \end{array} $

Then the nonlinear terms were all collected on the right-hand side. They applied a term T_a which is the duration of the transition period. They replaced right-hand terms into terms with T_a(named F and G), which smoothed the linear initial conditions to nonlinear computations. As the partial derivatives of the nonlinear terms F and G are equal to zero at any order, the procedure reduces the imbalance under the initial conditions. Dommermuth^[23,25] also used atmosphere pressure terms to balance the dynamic free-surface boundary conditions. It needs to be mentioned that the work of Refs. [23] and ^[25] requires the adjustment of dynamic boundary condition, while that of Ref. [24] requires adjustment of both dynamic and kinematic boundary condition.

The solution of the initial condition and less computational time of the HOS method make extreme waves with long time generation in open seas more available. Brandini and Grilli^[26] used the HOS method to do the first stage of growth of wave modulations and then applied the results as the initial conditions for extreme waves generation. It is worth noting that Brandini and Grilli presented the 2D HOS method as initial condition in transversal and longitudinal direction and it showed typical 3D wave modulations. Tanaka^[27] employed the HOS method as the basis numerical scheme, and in order to study the energy spectrum or time evolution of skewness and kurtosis of wave elevation, they introduced a complex amplitude function^[8]. The complex amplitude function b(k, t) is defined as^[8]

$ \begin{array}{*{20}{c}} {b(\mathit{\boldsymbol{k}},t) = \sqrt {\frac{{\omega (\mathit{\boldsymbol{k}})}}{{2\mathit{\boldsymbol{k}}}}} \hat \eta (\mathit{\boldsymbol{k}},t) + {\rm{i}}\sqrt {\frac{\mathit{\boldsymbol{k}}}{{2\omega (\mathit{\boldsymbol{k}})}}} \hat \phi (\mathit{\boldsymbol{k}},t)}\\ {\omega (\mathit{\boldsymbol{k}}) = \sqrt {gk} ,k = |\mathit{\boldsymbol{k}}|} \end{array} $

where $\widehat f\left(k \right)$ represents the Fourier transform:

$ f(\mathit{\boldsymbol{x}}) = \frac{1}{{2\pi }}\int {\hat f} (\mathit{\boldsymbol{k}}){{\rm{e}}^{{\rm{i}}kx}}{\rm{d}}\mathit{\boldsymbol{k}} $

ϕ_s and η are described with b(k, t) :

$ \begin{array}{*{20}{l}} {\eta (\mathit{\boldsymbol{x}}) = \frac{1}{{2\pi }}\int {\sqrt {\frac{k}{{2\omega }}} } \{ b(\mathit{\boldsymbol{k}}) + {b^*}( - \mathit{\boldsymbol{k}})\} {{\rm{e}}^{{\rm{i}}kx}}{\rm{d}}\mathit{\boldsymbol{k}}}\\ {\phi (\mathit{\boldsymbol{x}}) = \frac{{ - {\rm{i}}}}{{2\pi }}\int {\sqrt {\frac{\omega }{{2k}}} } \{ b(\mathit{\boldsymbol{k}}) - {b^*}( - \mathit{\boldsymbol{k}})\} {{\rm{e}}^{{\rm{i}}kx}}{\rm{d}}\mathit{\boldsymbol{k}}} \end{array} $

The complex amplitude b (k) is given by the value of k in discretized space, and for each k,

$ {b_k} = \sqrt {\frac{{{\omega _k}}}{{2k}}} {\hat \eta _k} + i\sqrt {\frac{k}{{2{\omega _k}}}} {\hat \phi _k},\quad {\omega _k} = \omega (\mathit{\boldsymbol{k}}) $

The energy density E can be given according to b_k on the wave field:

$ E \approx \sum\limits_k {{\omega _k}} |{b_k}{|^2} $

where E can also be given in terms of directional spectrum Φ(ω, θ) as

$ \begin{array}{l} E = \int\limits_0^{2\pi } {\int\limits_0^\infty \varPhi } (\omega ,\theta ){\rm{d}}\omega {\rm{d}}\theta = \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \int {\frac{{{g^2}}}{{2{\omega ^3}}}} \varPhi (\omega ,\theta ){\rm{d}}k \approx \sum\limits_k {\frac{{{g^2}}}{{2\omega _k^3}}} \varPhi (\omega ,\theta )\varDelta {S_k} \end{array} $

(7)

$ \varDelta {S_k} = \varDelta {k_x} \times \varDelta {k_y} $

(8)

Eqs. (7)-(8) gives the relation between |b_k| and Φ(ω, θ):

$ |{b_k}{|^2} \approx \sum\limits_k {\frac{{{g^2}}}{{2\omega _k^4}}} \varPhi (\omega ,\theta )\varDelta {S_k} $

(9)

Through Eq. (9), the directional spectrum Φ(ω, θ) can be determined by calculating the norm |b_k(t = 0)|, and the phase of b_k(t = 0) is determined by a random number in [0, 2π]. Φ(ω, θ) can be expressed with JONSWAP spectrum.

Ducrozet et al.^[28] followed the initializing condition of Dommermuth^[23] and simulated a forced focusing event wave based on the method of Tanaka^[27], and then a natural freak waves during a long time simulation were considered as well, as shown in Fig. 1. Their study gives a definite description on generating freak wave in an open-sea state, discusses the selection of parameters in directional spreading, and gives a proper solution to it. This encourages the exploration on the freak wave mechanics and the developments of the HOS method. Li et al.^[29] computed the sloshing wave in 3D tank which shows the efficiency of the HOS method in simulating focusing wave.

Zhao et al.^[30-32] applied an improved Adams-Bashforth-Moulton (ABM) numerical integration instead of Runge-Kutta forth (RK4) scheme to reduce simulation time and studied focusing wave simulation in experimental scale and real scale. The improved ABM scheme can reach a truncation error of O(Δt⁶), while RK4 scheme can only reach O(Δt⁵). However, as mentioned in Dommermuth and Yue's work^[5], they employed RK4 instead of ABM because RK4 can provide a lower global truncation error and a larger stability region. The improved ABM only increases the precise order, and they did not mention the global truncation error and numerical stability in their work. Therefore, it needs to be discussed in the future which numerical integration method has wider application scope.

Sergeeva and Slunyaev^[33] studied the characteristics of rogue waves by the HOS method. They produced space-time wave data sheets using the HOS method, and observed the detailed picture of individual rogue waves. They simulated 10 km field with 20 min sheets of data with surface elevation, fluid velocity, and so on. Fig. 2 illustrates rogue waves and rogue events, where markers indicate rogue waves and rectangles assemble them into rogue events. Rogue waves are difficult to observe for its long-time duration, short time appearance, and high crests with shallow troughs. Their study applies the HOS method to give the probability of appearance and lifetime of rogue waves, showing the possibility to study the rogue waves and wave-structure interactions in rogue waves. The application of the HOS method makes the study of rogue waves possible.

Xiao et al.^[34] studied the occurrence and dynamics of rogue waves employing the HOS method in three dimensions. They compared the HOS results to nonlinear Schrodinger equations, giving a quasi-stationary evolution in long-time prediction, which cannot be predicted by the nonlinear Schrodinger equations. They also introduced the definition of rogue wave occurrence which related to the occurrence probability, and discussed different situations of narrow spread seas. Along with the study mentioned above, the research on probability of rogue wave occurrence provides valuable reference for studying the structures which were emerged and changed in the rogue waves.

Those studies focused on unbounded domains of open-sea states, with periodic boundary conditions applied. However, unbounded domain means that the wave propagates without wave generation, absorption, and beach reflections, which makes experimental validations more difficult. Bonnefoy et al.^[35] developed a non-periodic HOS technique and presented a linear wave maker motion in the numerical wave tank. They considered an additional potential^[36] added on the original potential to make sure the potential satisfied the wave maker condition. By numerically solving the additional potential also in time domain, they extended the idea of Agnon and Bingham^[36], and adjusted a linear piston wavemaker to a generic wavemaker, as shown in Fig. 3. They presented 2D cases of focusing wave packet with extending second order initial decomposition and target phase signal iterative corrections.

Li^[37] employed the high-order-spectral tank (HOST) derived from Bonnefoy et al.^[35] to simulate 2-D focusing wave and made comparison with the experimental results. Li and Liu^[38] also applied the HOST method to generate directional focusing waves. Fig. 4 shows the comparison between experimental results and 2-D focusing wave generated by HOST. The surface elevation of 3-D focusing wave generation is displayed in Fig. 5. They studied the properties of the focusing wave and gave conclusions that shift of focusing point and maximum crest decreases when spreading parameter decreases in 3-D case. In addition, Li et al.^[39] used the HOST method to study the influence of current on the focusing wave. The potential velocity are expressed as the uniform current velocity added with wave potential velocity in their study. The numerical model revealed strong nonlinear interaction in the wave-current interaction.

Based on the previous works on HOST^[35] and its successful results compared with experiments, Ducrozet et al.^[40] developed this method into second-order wave maker motion. They combined the model developed by Bonnefoy et al.^[13,15] with the HOS method wave tank, and then extended HOST to second order. Fig. 6 shows the calculation process of this new HOST. A 2D focusing wave packet embedded in irregular wave fields was considered and showed apparent improvement in results compared with those in linear HOST.

Bonnefoy et al.^[41] applied both second-order wave maker motion and fully nonlinear numerical wave tanks and gave the validation results to experimental complex sea-states in 3-D.

Those improvements and developments of the HOS method lays the foundation for more application of the HOS method. In order to make the method more conveniently applied in many fields, Ecole Centrale Nantes, LHEEA Lab developed an open-source HOS model^[42] named HOS-Ocean. The software is available on the Github (https://github.com/LHEEA/HOS-ocean/wiki). The HOS method is realized in a numerical simulation way, and provides possibility to couple the CFD software.

Assuming the simulation based on a rectangular fluid domain with periodic boundary condition, the velocity potential of surface boundary condition can be written as Eqs. (5)-(6), where k_n = n(2π)/L_x is the wave numbers. Fig. 7 presents the brief description of code working in HOS-Ocean according to Ref.[42]. For the wave spectra used in HOS-Ocean is only JONSWAP spectrum, Song et al.^[43] made an attempt to add two new wave spectrum P-M spectrum and ITTC 2-parameter spectrum in it.

Meanwhile, Ducrozet et al.^[44] developed the open-source numerical wave tank based on HOS (HOS-NWT). The HOS-NWT is also available on Github (https://github.com/LHEEA/HOS-NWT/wiki). They extended wave maker modeling into three-order generation.

Seiffert and Ducrozet^[45] included the wave breaking mechanism into HOS-NWT. As the HOS method assumes free-surface as single-valued, they approximated the broken surface as a single value. They suggested a threshold of water particle velocity to crest velocity between a broken and an unbroken wave in HOS solvers. Seiffert et al.^[46] implemented a wave breaking onset criteria in HOS-NWT, which predicts good wave-breaking onset time and location compared with the experiments. Seiffert and Ducrozet^[47] validated this wave-breaking mechanism through energy dissipation by adding an eddy viscosity parameter under boundary conditions. This kind of wave breaking onset can be incorporated with other nonlinear potential solver to reduce the risk of numerical instabilities, thus increasing the scope of application.

The development of the HOS method shows a perspective in ocean and ship engineering. The ability of generating large-domain wave field with fast convergence and accuracy gives researchers the chance to study rogue waves. The difficulties in studying rogue waves are due to its transient time appearance and long-time with large space simulation, which requires huge investment in resources. With the help of the HOS method, the parameters of rogue waves and the position of its occurrence can be found out. This lays the foundation for simulating severe conditions in ocean and sail engineering.

The applications and extensions of the HOS method are not limited to wave simulations. Many physical mechanisms are included by employing and validating the the HOS method, which include nonlinear energy transfers^[48], bi-modal seas^[49], modulation instabilities^[50-51], and so on (cf.^[42]). There are also some developments based on the HOS method, including wind forcing effects^[52], jet current effects^[53], variable bathymetry^[54-58], and so on.

3 Combined Method of the HOS and CFD Method

In order to simulate the wave-structure interaction in a more realistic and flexible way, the CFD method is carried out. This approach functions extremely well in resolving complex surface problem, large amplitude motion of structures, and so on. However, even though the development of computational science is rapid, the time duration and cost of the CFD method is still nonnegligible. In fact, when wave-structure interaction is considered, wave loads extending on the structures can be focused even more. Besides, the viscous effects are not obvious in wave simulation, and long-time wave simulation in the CFD method may cost a long time and cause numerical dissipation. Therefore, combining the potential and viscous flow together can improve the simulation in an efficient and realistic way.

Li^[59] gave a view of two groups which divides the combined potential and viscous flow: Domain Decomposition and Function Decomposition, the schematic diagram is shown in Fig. 8. He mentioned Domain Decomposition sets viscous domain and potential domain apart, while Function Decomposition sets viscous part and potential part in viscous domain.

3.1 Function Decomposition

Many studies have been conducted^[60-64] to solve a wave-structure interaction employing function decomposition approaches. These studies applied BEM as potential solver and RANS or Lattice Boltzmann Method (LBM) as viscous solver. Later on, a new approach was carried out focused on solving diffracted flow only, which is named Spectral Wave Explicit Navier-Stokes Equations (SWENSE). It is derived from a frame of potential theory^[65-66] and splits all unknown variables into incident terms and diffracted terms^[67]. SWENSE modifies the variables in RANS equations as follows^[68]:

$ \left\{ {\begin{array}{*{20}{l}} {{u^\alpha } = u_{\rm{I}}^\alpha + u_{\rm{D}}^\alpha }\\ {p = {p_{\rm{I}}} + {p_{\rm{D}}}}\\ {h = {h_{\rm{I}}} + {h_{\rm{D}}}} \end{array}} \right. $

where α∈{1, 2}, the unknown variables u, p, and h stands for Cartesian components of velocity, pressure, and free-surface elevation, respectively. The subscripts I and D represent incident and diffracted variables, respectively.

According to Navier-Stokes equations^[59],

$ \begin{array}{*{20}{c}} {\nabla {u^\alpha } = 0}\\ {\frac{{\partial {u^\alpha }}}{{\partial t}} + {u^\alpha }\nabla {u^\alpha } = - \frac{{\nabla p}}{\rho } + \mathit{\boldsymbol{g}} + \mathit{\boldsymbol{v}}{\kern 1pt} {\kern 1pt} {\nabla ^2}{u^\alpha }} \end{array} $

and Euler equations^[59],

$ \begin{array}{*{20}{l}} {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \nabla u_{\rm{I}}^\alpha = 0}\\ {\frac{{\partial u_{\rm{I}}^\alpha }}{{\partial t}} + u_{\rm{I}}^\alpha \nabla u_{\rm{I}}^\alpha = - \frac{{\nabla p}}{\rho } + \mathit{\boldsymbol{g}}} \end{array} $

The transport equations are obtained by substituting the Euler equations from the NS equations^[59],

$ \begin{array}{l} \frac{{\partial u_{\rm{D}}^\alpha }}{{\partial t}} + (u_{\rm{I}}^j + u_{\rm{D}}^j)\frac{{\partial u_{\rm{D}}^\alpha }}{{\partial {x^j}}} - v\frac{{{\partial ^2}u_{\rm{D}}^\alpha }}{{\partial {x^{{j^2}}}}} + \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \begin{array}{*{20}{l}} {\frac{1}{\rho }\frac{{\partial {p_{\rm{D}}}}}{{\partial {x^\alpha }}} = u_{\rm{D}}^j\frac{{\partial u_{\rm{I}}^\alpha }}{{\partial {x^j}}} + \mathit{\boldsymbol{v}}\frac{{{\partial ^2}u_{\rm{I}}^\alpha }}{{\partial {x^{{j^2}}}}}}\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \frac{{\partial u_{\rm{D}}^j}}{{\partial {x^j}}} = 0} \end{array} \end{array} $

With those equations, the incident variables were explicitly computed, thus named the solver SWENSE Equations. The incident nonlinear regular wave was obtained by a stream function model^[69] while the nonlinear wave trains were derived from a spectral formulation combined with BEM^[70] first, and then from the HOS method. The SWENSE method using the spectral formulation combined with BEM was validated both in 2D cases^[71-72] and 3D cases^[73].

Luquet et al.^[74] employed the HOS method to play the role as incident wave model. Compared with the previous potential theories, the HOS method makes the simulation of arbitrary complex sea-state or freak waves more convenient and efficient. They used this method to simulate a Tension Leg Platform (TLP) in regular and irregular waves.

Other studies such as ship maneuvering with 6DOF, CALM buoy, and ship motion in 2DOF have been done in SWENSE model^[75-78]. A validation of ship motion in SWENSE in irregular head wave was carried out (Fig. 9).

However, those original SWENSE methods are based on one-phase CFD solver, which solve the free surface through the boundary condition. Li et al.^[79-80] extended the SWENSE model into two-phase SWENSE model. They modified Euler momentum equation by introducing an incident pressure to balance the governing equations both in water and air phases.

The SWENSE model divides the total field into two parts: incident field and diffraction field. In this way, the incident field was solved by potential theory, thus the computational grids can be a coarse mesh while only the grids near diffraction field need to be refined. It not only improves efficiency when using the HOS method as incident model, but also decreases the cost when simulating wave structure interactions due to less grids in far field.

Gatin et al.^[81-82] combined HOS with CFD RANS equation to solve 6DOF ships in extreme waves, as shown in Fig. 10 and Fig. 11. The combination was performed based on open source CFD software OpenFOAM^[83]. They obtained free surface elevation and velocity from the HOS method, thus the velocity field is presented as^[81-82]

$ \begin{array}{l} {v_x}(x,y,z,t) = \sum\limits_k {\sum\limits_l {{c_{k.l}}} } (t){\rm{i}}{K_k}\frac{{{\rm{cosh}}({K_{k.l}}({z^\prime } + d))}}{{{\rm{cosh}}({K_{k.l}}d)}}{{\rm{e}}^{{\rm{i}}{K_k}x}}{{\rm{e}}^{{\rm{i}}{K_l}y}}\\ {v_y}(x,y,z,t) = \sum\limits_k {\sum\limits_l {{c_{k.l}}} } (t){\rm{i}}{K_l}\frac{{{\rm{cosh}}({K_{k.l}}({z^\prime } + d))}}{{{\rm{cosh}}({K_{k.l}}d)}}{{\rm{e}}^{{\rm{i}}{K_k}x}}{{\rm{e}}^{{\rm{i}}{K_l}y}}\\ {v_z}(x,y,z,t) = \sum\limits_k {\sum\limits_l {{c_{k.l}}} } (t){\rm{i}}{K_{k.l}}\frac{{{\rm{cosh}}({K_{k.l}}({z^\prime } + d))}}{{{\rm{cosh}}({K_{k.l}}d)}}{{\rm{e}}^{{\rm{i}}{K_k}x}}{{\rm{e}}^{{\rm{i}}{K_l}y}}\\ {\rm{where}}{\kern 1pt} {\kern 1pt} {\kern 1pt} {z^\prime } = zd/\left( {d + \eta } \right) + d\left( {d/\left( {d + \eta } \right) - 1} \right). \end{array} $

where z' = zd/(d+η)+d(d/(d+η)－1).

They coupled the HOS method with the the CFD method using a two-phase SWENSE method^[84-86] in level set method^[87]. This two-phase SWENSE method is different from that in Li^[79-80]. It deals with incident wave as a source term in CFD meshes. Therefore, the incident domain in CFD mesh still needs to be refined to keep accurate.

3.2 Domain Decomposition

The domain decomposition sets computational domains apart into viscous flow sub-domain and potential flow sub-domain. In general, with advantages of solving complex free surface and viscous effect around structure, the CFD method takes the inner domain and potential flow method takes the outer domain. The two methods solve separately in each sub-domain and communicate through the common boundary. Therefore, the domain decomposition method can be two-way coupling or one-way coupling.

Before the emergence of the HOS method, either two-way coupling^[88-91] or one-way coupling^[92-94] employed BEM method in potential flow sub-domain. Kim et al.^[95] applied an overlay zone to do the communication between CFD and Euler wave, which is named Euler-Overlay Method (EOM). They combined EOM with mooring and riser model to do the simulation of Vortex-Induced Motion (VIM). Paulsen et al.^[96] combined Oceanwave3D^[97] with the CFD method and gave validations under different sea bed conditions.

Choi et al.^[98] combined the HOS method with OpenFOAM (foamStar) through a HOS wrapper program named Grid2Grid^[99] which can be obtained from GitHub (https://github.com/LHEEA/Grid2Grid). Grid2Grid reconstructs HOS wave results and interpolates at demanded time and position in viscous flow field. The frame of program of Grid2Grid is shown in Fig. 12. Choi et al.^[98] tested the new solver in empty wave tank and validated the stable and accuracy of this solver.

Zhuang et al.^[100] developed the HOS method into an in-house CFD solver naoe-FOAM-SJTU^[101-108]. They applied relaxation zone^[109] as the common boundary to communicate information, as shown in Fig. 13. They tested the ability and accuracy of the new combined solver in empty wave tank, and gave good results. Song et al.^[110] also applied this new solver to simulate freak waves in the empty tank. The CFD zone can obtain wave elevation and velocity from the HOS method in arbitrary time spot and positions, thus the new combined solver reduces the simulation cost in the CFD method.

The CFD solver naoe-FOAM-SJTU is based on OpenFOAM and develops many modules and functions in recent years, such as 6DOF module^[111], overset grids^[112-114], internal flow coupled with external flow^[115-117], and so on. Combined with the HOS method, the new solver has advantages on solving structures in open-sea states in an efficient way. A number of simulations have been carried out^[118] by the new combined solver, presenting the ability to simulate the cylindrical structure in multi-directional sea-state (Fig. 14) as well as focusing wave (Fig. 15).

In general, a ship sailing in the seas is assumed to be sailing in the encountered added current. But a ship sailing in a large domain constitutes realistic ship maneuvering. With the help of overset grids and the HOS method, ship can sail in open seas without large computational domain requirement, as shown in Fig. 16. KCS ship model in model scale is chosen as the ship model and the Lpp of the ship is 6.070 2 m. The ship speed is 2.017 m/s, and the wave is first order Stokes wave in λ/L = 1.15. Through the results of heave motion and pitch motion as shown in Ref. [118], the coupled method is closer to experimental results than pure CFD method.

Quinn^[119] combined HOS-NWT with OpenFOAM to study the breaking waves. The HOS-NWT replicates a 2D waves conducted in experiments with a 3D long-crest irregular wave and the results were analyzed with different breaking onset criteria. However, the breaking process was not visible in the coupled method. Quinn pointed out that this may due to the diminishing of crest value in HOS-NWT, which provided flows in OpenFOAM. Meanwhile, Quinn also simulated freak and rouge waves in HOS-NWT, and gave a conclusion that with the increase of the wave steepness, the discrepancy between HOS results and experimental results becomes large.

The review of application with the coupled method gives a huge possibility in wave-structure interactions. Nowadays, with the help of overset grids, the large amplitude of structures is achieved in the CFD method. The coupled method solved a big problem which the CFD method cannot solve at this moment. The HOS method gives the input information of large-scale or large domain of wave or sea states, the CFD method can simulate any situations in ocean and naval engineering. The given results of the coupled method are not enough though. Many research such as large-scale simulations, fine wave fields, and ship sailing in freak waves need to be done.

4 Conclusions

The HOS method has been developed in the 1980s. With the feature of efficient wave generation and fast numerical convergence, an increasing number of applications using the HOS method have been carried out. The HOS method is popular in recent years for its ability in simulating large-domain wave generation and long-time evolution. From the very beginning, the HOS method was applied for studying nonlinear gravity wave and wave-wave interaction. After that, the attempt in studying wave-body interaction in the HOS method was carried out. The extreme wave and long-time simulation was completed from 2D to 3D, along with the observation of occurrence of freak wave in open-sea state. Then the application of the HOS method extended to bounded boundary condition which is wave tank simulation. With the development of computer science, the open source software of using the HOS method is well developed. Although many studies have been conducted to eliminate numerical error or employed assumptions to solve the shortfalls, development of the HOS method should not be stopped.

The combined method of the HOS and CFD method is a new theory in numerical simulations. The idea of combining the HOS with CFD method was put forward in recent decades. The development and application of the integrated method is all about fixed structure to test the accuracy of the combination. Some studies have been done to simulate ships in freak waves, showing the ability of the combined method in dealing with large motion of structure in freak waves.

Although the combined method of the HOS and CFD method only rises in recent years, its achievement is a progress in numerical simulation. The ability and advantages of the combined method are listed as follows:

1) Because the HOS method can give incident wave information in arbitrary space and time spot, the CFD method can choose required position and time spot in wave propagation. For example, the HOS method can be employed to do a real-sea state generation. After observing the appearance of freak wave, the CFD method can be used to receive the freak wave signal to simulate the freak wave-structure interactions. It is well known that the capture of freak wave needs thousands of seconds of simulation, which is impossible for the CFD method, while the combined method provides the possibility to do the freak wave structure interactions.

2) The large-scale model simulation in the HOS method coupled with the CFD method can also be achieved. The HOS method has the advantage in generating wave in large field efficiently, thus it can be an economic way to simulate large-scale model, especially in two-way coupling. The domain of the CFD method can be minimized around the large-scale model, and the wave field outside can be solved by the HOS method. Therefore, it will save lots of time and resources compared with pure CFD method.

3) The combined method has the ability and advantages in simulating fine wave field. When the simulations are based on focusing wave or fully developed nonlinear wave, the CFD method needs to spend lots of time to make sure the nonlinear wave or focusing wave are fully developed. The viscous effects in the CFD method may diminish the wave amplitude during the wave propagation. In this way, it costs splendid time and resources when the fine wave field is considered. With the help of HOS wave field which is a fully developed wave field, the time during wave propagation can be ignored in the CFD method. Therefore, the meshes in CFD can be as fine as possible to achieve the accuracy as required.

As stated above, the combined method has great ability in simulating some kinds of wave-structure interactions, which include ships or platforms in freak waves that may encounter impact pressure, breaking wave, and air gap phenomenon; large-scale model simulations; realistic wave-structure interactions such as real-sea state simulation, and so on.

5 Future Works

Some research perspectives and corresponding solutions are given below.

1) There are two kinds of coupled method of the HOS and CFD method: one is function decomposition of two-phase SWENSE^[59,79-80], the other is SWENSE in Naval hydropack^{[81-82,84-86]} and domain decomposition^[98,100]. The first two-phase SWENSE model is limited to fixed structure up till now, therefore a 6DOF module and dynamic mesh needs to be applied. These modules can be implemented through OpenFOAM. The second SWENSE and the domain decomposition have the 6DOF module and dynamic/overset grids, but they have to refine the mesh in incident wave domain.

2) For coupled method of the HOS and CFD method in domain decomposition, it is one-way coupling at present. It does not need extra iterations between potential and viscous flow. However, it requires larger domain to realize complete radiation flow of wave-interaction influence than that in two-way coupling. Considering the ship maneuvering in the coupled method, the one-way coupling will end with misplaced wave patterns in the view of whole domain. Nevertheless, the two-way coupling of the HOS and CFD method needs to change the formulation of the HOS method and improvement of the common boundary condition, thus it is still under consideration.

3) The simulations based on the HOS method^[28,42,44] are conducted on single-processor, while simulations in the CFD method often need multiple processors. This leads to separate calculation in the HOS method and the CFD method. In domain decomposition, the CFD calculation in multiple processors will duplicate the results of the HOS method. When the internal memory of HOS results is large, the duplication in parallel running may cause a huge work in CFD simulation. The HOS method has no apparent solutions on parallel running at present, therefore the CFD method can solve this problem by sparing an extra processor for the HOS method. This will prevent the duplication of HOS results files.

4) The studies of breaking waves are often included in wave-structure interactions. With the assumption of single value of breaking waves^[45-47], the HOS method has the ability to predict breaking wave onset criteria. However, the assumption is not realistic and cannot be obviously visible in CFD domain. Therefore, the study of breaking phenomena needs to be done in viscous flow in the future. The CFD method can capture arbitrary time spot of HOS results, thus the thought is to map the breaking time spot of HOS to CFD domain, thus the breaking in waves happens under viscous flow solution. It requires further tests, because the incident information of the HOS method is absent, and the breaking wave in CFD domain may collapse or propagate in a wrong way.

5) There are some numerical instabilities with the combined HOS and CFD method. As all the variables are solved on free surface in the HOS method, the interaction with the CFD method will cause instability in sharp free surface. For function decomposition of two-phase SWENSE^[59,79-80], it applies Ghost Flow Method (GFM)^[120-122]to keep discontinuity in the CFD method. For domain decomposition of the combined method, the accumulated velocity on sharp free surface can be solved by an adding inversed velocity.