Extreme Wave Simulation with Iterative Adaptive Approach in Numerical Wave Flume

Extreme wave is highly nonlinear and may occur due to diverse reasons unexpectedly. The simulated results of extreme wave based on wave focusing, which were generated using high order spectrum method, are presented. The influences of the steepness, frequency bandwidth as well as frequency spectrum on focusing position shift were examined, showing that they can affect the wave focusing significantly. Hence, controlled accurate generation of extreme wave at a predefined position in wave flume is a difficult but important task. In this paper, an iterative adaptive approach is applied using linear dispersion theory to optimize the control signal of the wavemaker. The performance of the proposed approach is numerically investigated for a wide variety of scenarios. The results demonstrate that this approach can reproduce accurate wave focusing effectively.


Introduction
Extreme wave, also called freak wave, which is unpredictable in the open sea (Osborne et al., 2000), is always known as an event with highly asymmetric and steep wave crest and leading to significant damages to ships, offshore and coastal structures. Many catastrophic accidents caused by extreme waves have been reported in the last few decades (Didenkulova, 2020;Lavrenov, 1998;Nikolkina and Didenkulova, 2012). As a consequence, studies on extreme waves and their impact are still active in the academia.
To understand the complex phenomenon, associated research has been conducted on the physical mechanisms of extreme wave generation. Trulsen (2001) conducted the work on formation process of the New Year wave using nonlinear Schrödinger (NLS) equation and found that the freak wave did not appear suddenly but gradually evolved from an intense wave group, which suggests that the dispersive focusing mechanism played an important role in the generation of extreme waves. Divinsky et al. (2004) researched the extreme waves which occurred in the Black Sea using the NLS equation and found that the extreme wave also appeared from the evolution of wave packet with a long time, while the extremely high wave just emerged in a short limited distance. Slunyaev (2006) simulated several extreme waves measured in the North Sea during a storm by the method of the linear model, the NLS equation and the Dysthe equation, respectively. The extreme waves were also found to occur gradually from intensive wave packets with a mechanism of modulational instability. Both the dispersive focusing and modulational instability can contribute to large waves. However, as noted by Slunyaev (2006), weakly nonlinear models, such as the NLS equation and the modified NLS equations, struggle to predict the evolution of steep wave trains. Therefore, numerical models with higher-order nonlinearity should be employed. Goullet and Choi (2011) found out that the high-order spectral method (HOSM) with third-order nonlinearity is a reliable method in simulating the evolution of nonlinear wave packets. Recently, Wang et al. (2020) investigated the evolution of extreme waves with HOSM demonstrating that the third-order non-linearity changes the dispersion relationship of the wave components. As a result, the higher-order models should be used for strong nonlinear wave packets. To date, the existing methods for the simulating of extreme waves usually apply linear wave theory (Baldock et al., 1996;Janssen, 2003;Swan, 2001, 2003;Liu et al., 2016;Rapp and Melville, 1990). The problem of harmonics does not satisfy the linear dispersion relation, resulting in extreme wave occurring prematurely compared with theoretically predicted result (Liu et al., 2016). This means that linear dispersion relation is employed to back scatter the target wave elevation to the wave maker, when a wave packet is generated, nonlinear wave-wave interaction leads to obvious increase of wave crest, therefore extreme wave occurs before defined focusing position. Additionally, the presence of spurious components generated by the wavemaker can also contribute the wave focusing prematurely (Sriram et al., 2013). Hence, the study of an accurate extreme wave occurring at the desired position in the flume is very important.
In this study, an iterative adaptive approach is proposed to achieve the specified wave profile allowing better control over the generation of extreme wave with different steepness and frequency spectra. The objectives of this paper are: (a) to analyze the parameters of the wave group affecting the focusing position, (b) to demonstrate the capability of the iterative adaptive approach to reproduce accurate wave focusing based on high order spectral method in the present study, (c) and to assess the performance of iterative adaptive approach by means of using different parameters.

Governing equations
It is assumed that the fluid is inviscid and incompressible, the flow is irrotational during the whole study. The domain and coordinate system are shown in Fig. 1. A Cartesian coordinate system is employed where the z = 0 line corresponds to the still water level, z direction is positive upwards and x direction follows the propagation of wave. The origin is defined at the mean position of the wave-maker. For the two-dimensional surface wave propagation problem, the velocity potential Ф (x, z, t) in the fluid domain satisfies the Laplace equation, i.e., (1) The fully nonlinear free surface boundary conditions expressed by the velocity potential on the water surface can be written as: where P a is the atmospheric pressure and can be set zero here, . As usual, lateral boundaries of the domain are either fixed or moving ones. For this wave flume, waves are generated by a wave-maker with motion X(t) at the left boundary. As a result, the boundary condition for wave making with linear method can be written as: On other fixed boundaries (bottom, side and end walls of the tank), the condition may be simply written: n where is a vector normal to the considered boundary.

Numerical procedure
In this study, the equations are solved by High Order Spectral (HOS) method proposed by Dommermuth and Yue (1987). For the numerical wave flume, the wave making system is essential, which introduces the non-homogeneous boundary condition leading to the HOS method being difficult to use. In this way, the velocity potential Φ will be split into the sum of an unknown period component Φ f and a prescribed non-periodic components Φ w according to Agnon and Bingham (1999) where both the prescribed non-periodic components Φ w and period component Φ f satisfy the Laplace equation and Eq.
(2) and (3) and following Zakharov (1968), the free surface boundary conditions can be written as: A spectral expansion method following Bonnefoy et al. (2004) is employed to simulate the non-periodic components, which has the advantage of less calculation and application for three-dimensional wave basin. Therefore, the nonperiodic components can be written as: where k q = (2q−1)π/(4h), h is the water depth and L x is the length of the two-dimensional wave flume. The unknown amplitude a q can be obtained by the wave maker boundary Eq. (4) using the method of FFT.
As to the period component Φ f , it can be calculated by the HOS method and the procedure is briefly shown as follows.
With Taylor expansion, the velocity potential can be expanded into a series of potentials of order m on z = 0, The Dirichlet boundary condition is satisfied when Φ f and η are known. In this way, can be expanded as follows to solve the Dirichlet problem, where the spectral modes Ψ n satisfy periodic conditions and can be expressed as: a (m) n where k n = (nπ)/L x . The amplitudes can also be obtained by FFT like a q . At this time, Eqs. (7) and (8) can be solved by Rung−Kutta scheme step by step.

Wave maker and absorbing boundary
For a wave maker boundary, the velocity potential Φ w should be pre-solved before the time integration is conducted. According to linear wave maker theory, the incident velocity can be calculated by the following equation for a specified wave.
where η is the expected water waves. ω = 2πf is the circular wave frequency satisfying dispersion equation. T(k) is the transfer function and can be calculated by the following equation for piston type wave-maker.
where k is the wave number. It is worth mentioning that a ramping function is applied to the water surface elevations so that the flow will change slowly when the computation starts.
At the end of the flume the sponge layer proposed by Larsen and Dancy (1983) is used to mitigate reflections. In the sponge layer the surface and velocity potential are divided by μ(x, y) after each time step. The expression of μ(x, y) is given as follows: where d z is the width of sponge layer; d is the distance between the points in sponge layer and the boundary; Δd is the dimension of the grids; α and n are the constant parameters which are 4 and 5, respectively.

Generation of extreme wave
In this study, the dispersive focusing method was employed to generate the extreme waves. According to the linear wave theory (LWT), the surface elevation of different frequencies at an arbitrary point (x, y) can be expressed as: where N is the total number of wave components being 64 in this study; and a i , k i , ε i and ω i are the corresponding wave amplitude, wave number, phase shift and radian frequency of the i-th wave component, respectively. The wave number of each component k i is given by the linear, gravity wave dispersion relation ω where g is the gravitational acceleration and h is the water depth. It was assumed that the waves were focused at a spatial position x f and at a specified time t f . Hence, at the focused point, the phase shift for each wave component is ε By substituting Eq. (18) into Eq. (16) and letting n = 0, the wave elevation at an arbitrary point can be written as: That is, the theoretical displacement of the paddle for generating focused waves can be described as: (20) 2.5 Verification of numerical model Before conducting the numerical simulations, a convergence check has been conducted (Li et al. , 2012) to ensure the numerical approach accurately simulates the physics. Thus, the numbers of points in the x direction used for simulation work during this study is 2000, where an acceptable result can be obtained. Additionally, in order to assess the accuracy of this numerical model, classical test cases are first used which involve the transformation of 5-th order Stokes waves in deep water. The details of parameters are shown in Table 1. Stokes waves of period T = 2 s with different steepness were calculated in the horizontal part of a wave flume with water depth h = 4 m. For each case, data were collected with a record length of 81.92 s and a sample interval of 0.02 s. As seen in Fig. 2, synchronized time series of free-surface elevations are available at location x = 10 m compared with theoretical results. It shows that they have a good overall agreement. However, the discrepancies began to be obvious when wave steepness ε is larger than 0.15. It is worth noting that the time series are not in phase well indicating that the linear dispersion employed here may contribute the error. In addition, the linear wave maker boundary used here may also contribute some discrepancies, for which the second-order wavemaking theory may be studied further.
Another challenging test cases are then considered, and the ability of the numerical model to describe wave focusing process involving the effect of the distribution of frequency spectra and increasing of wave steepness can be seen in Table 2. The first kind of frequency spectrum adopted is called the constant-wave-steepness (CWS) distribution and is a spectrum with a sharp peak which can obtain a steady focusing wave group and inhibit breaking prematurely according to Perlin et al. (1996). In this case, the steepness of each component is assumed to be constant, i.e. k i a i = const. The second is called the constant-wave-amplitude (CWA) distribution, that is a spectrum with a "top hat" shape, which means that the amplitude of each component a i is constant. In addition, a focusing wave experiment whose focusing point x f and time t f were 7 m and 15 s respectively was also conducted in a two-dimensional wave flume of 50 m long with water depth of 0.45 m for comparison in the State Key Laboratory of Coastal and Offshore Engineering, Dalian University of Technology. Figs. 3−6 show the comparison of wave surface elevations along the flume between experimental and numerical results. It is obvious from the figures that the numerical model employed here can generally simulate the focusing process well, particularly the main waves were predicted much better. However, it seems that the width of frequency band and dis-   tribution of frequency spectrum have more influence on the simulating result, while the trough elevations near focusing position are much larger for numerical results. It seems possible that this result is due to the strong nonlinearity near focusing position with a significant effect on the evolution of the surface elevations.
In summary, these two different examples can demonstrate that the numerical model has the capability to simulate wave generation and focusing and can be used to study further.

Numerical study of generation of extreme wave
In this section, the extreme wave was generated by the method of wave focusing. As is discussed above, the char-acteristics, i.e. frequency width, distribution of frequency spectrum as well as focusing amplitude, having an obvious effect on wave focusing were investigated in this part. The detail information is shown in Table 3. In all cases, the wave steepness of each component k i a i is smaller than 0.15. It is assumed that the extreme wave occurred at the position x = 7 m and at the time t = 15 s according to linear theory. The wave surfaces were collected with a record length of 81.92 s and a sample interval of 0.02 s along the flume seen in Fig. 1.

Focusing position and relative phasing
Predicting occurrence of extreme wave is one of the major challenges in the study of ocean waves. Thus, the meth-  od of wave focusing considered as a good manner to study the feature of extreme wave is usually adopted by scholars, and waves with the different frequencies can focus at a spatial position as control. Actually, the real focusing position also changes during the experimental and numerical studies using the linear focusing theory (Baldock et al., 1996). Fig. 7 shows the shift of focusing position with different parameters. It can be seen that the real focusing position downstream shifts and the shift increases with the increase of focusing amplitude. For the case with smaller focusing amplitude (A = 0.01 m), it is obvious that there is little shift occurring indicating that the effect of nonlinear interaction between different wave components is negligible. It is also consistent with the experimental results of Liu et al. (2016). What is more, the shift of the focusing position is also depend on the distribution of frequency spectrum, among   which the CWA spectrum contributes more shift. In addition, the width of frequency band is another main factor that can affect the wave focusing process. It can be found in Fig. 7 that the shift increases largely as the frequency bandwidth decreases. These results suggest that the shift of focusing position is closely related to the focusing amplitude, frequency bandwidth as well as frequency spectrum, which means that the change in the position focusing depends on the global nonlinear interactions of the wave group (Johannessen and Swan, 1997). It is clear from the shift of focusing position that there is a redistribution of wave components within the focused wave group. To investigate this behavior, the relative phasing of wave components can be determined from a Fourier transform of the wave surface elevation. In this paper, the method of all-phase FFT (Huang et al., 2008), which can achieve high-accuracy phase measurement without any correcting operations that the computation amount is decreased largely, is adopted to determine the phase of every component. This method can prevent spectrum leakage effectively so that the initial wave phase is invariant. The surface eleva-tion was firstly windowed by a (2M−1)-length window, which (2M-1) is equal to 4 095 in this paper. Then every two samples with the space between M intervals (except the middle element) were summed up thus a new M-length data can be formed. Finally, the FFT was implemented on the new data to generate the outcome of all-phase FFT. Fig. 8 shows the wave phase changing along with the flume concerning the CWA spectrum (Case 1) with a relatively smaller focusing amplitude. The phase of the free waves is considered at four locations fixed relative to the linear focus position. In this case, it is obvious that the simulated data are in good agreement with the linear result and in particular, the phases are matched at the focusing position. It is consistent with the result above, which suggests that the linear dispersion relation works well when the focusing amplitude is small. In contrast, Fig. 9 again concerns the CWS spectrum (Case 2) and express the similar result corresponding to the focusing amplitude of A = 0.05 m. Near the wavemaker, the wave phases obtained from the simulated data still match well with the linear results presented in Fig. 9a. However, as the waves propagating downstream the nonlinear wave-wave interactions dominate and thus, the relative phase of wave components change. As a result, the waves are no longer in phase at the focusing position. This result is consistent with previous conclusions of the focusing position shift shown in Fig. 7. Additionally, it is interest that the phase changes exhibit considerable variation in high frequency band (f > 1 Hz) for this wide band frequency, while the phase mostly changes in the low frequency. These local phase changes indicate that the phase velocity of wave components has already changed by the nonlinear interaction. Indeed, Fig. 9c provides that the phase velocity increasing identified via Fourier transform for the high frequency band compared with small amplitude case  LIU Dian-yong China Ocean Eng., 2021, Vol. 35, No. 1, P. 61-71 ( Fig. 8c). In the meantime, Fig. 10 shows another phase change for the case with narrow-band frequency. Similarly, the phase exhibits little variance in comparison with the linear result when the probe is located near wavemaker. But in this case, it can be seen that the phase of all wave components in the frequency band changes, which is consistent with the theoretical investigation revealed by Wu et al. (1979) concluding that the velocity of nonlinear wave group increases for steep narrow-banded waves. This may explain why the narrow band frequency has more significant effect on the shift of focusing position.
In the end, it is obvious that the wave phase change by nonlinear wave-wave interaction downstream of the focusing position is perhaps permanent, which has been demonstrated by Peregrine (1983) that the nonlinear interaction of waves travelling at different velocities produces a permanent phase change.
3.2 Generation of extreme wave at fixed position From the results presented above, it can be known that it is the phase of wave component changes resulted from nonlinear interactions that contribute to the shift of focusing position. By this way, if the phase matches during the nonlinear wave-wave interactions the extreme wave would occur at the fixed position as we need. Therefore, it is obvious that the linear dispersion relation used for the generation of strong nonlinear wave group is no longer suitable. Madsen and Fuhrman (2012) made use of the third-order theory to simulate multi-directional irregular wave, which incorporates the effect of nonlinear dispersion relation. Through this  manner, the shift of focusing position may disappear and it will be studied in the future. In this paper, an iterative adaptive approach was proposed based on all-phase FFT algorithm to modify the phase at the theoretical focusing position, where the detailed steps can be applied as follows: Firstly, the target wave profile at the predefined location is generated by means of linear focusing theory without considering any nonlinear interactions. Secondly, the all-phase FFT spectral analysis method is used to estimate the phase difference of each wave components between simulation and linear results using above wave elevations. Finally, the phase difference from Step 2 is back transformed to Eq. (20) to modify the initial phase of the paddle for generating focusing waves. After several iterative procedures, the control signal will be improved and a required target wave is obtained at the predefined location. Thus, the missing nonlinear interaction that was not taken into consideration will be automatically taken care. Fig. 11 shows the comparison of the phase between linear and numerical data after iterative adaptive approach used for the case with focusing amplitude A = 0.05 m. Clearly, the phases of wave components are in good agreement after one iteration (red circles in Fig. 11) with linear result at the theoretical focusing position, where the standard deviation is 0.044, implying that the fluctuation of wave phases is very small. The result indicates that the iterative adaptive approach simulates well the wave phase in the focusing event. In addition, as shown in Fig. 11d, the simulated wave phase downstream of the focusing position matches well the linear case. However, the wave phases for relatively high frequencies near the wavemaker from numerical and linear data generate large discrepancies compared with Fig. 9a. This further demonstrates that high frequencies play a significant role in nonlinear wave-wave in-teractions contributing to phase mismatch. Fig. 11 also provides a comparison after three iterations (red triangles in Fig. 11), showing that the wave phases obtained by simulation are almost the same as the linear result. The standard deviation at the focusing position is 0.005 at this time. Overall, the iterative adaptive approach does a great job in simulating the focusing process at fixed position. Fig. 12 shows the surface elevations at the focusing position and the linear data are also imbedded for comparison when the iterative adaptive approach is used. As shown in the figure, great focusing results are achieved for different frequency spectra and frequency bandwidth. But it can be seen that there are large differences between the numerical result and linear solution. The crest from simulation become higher and narrower, while the troughs become wider and shallower compared with linear solution, suggesting that the nonlinear interaction between the wave components has a significant effect on the wave evolution. While comparing the wave surface elevations simulated by iterative adaptive approach, they are almost the same, which means that only one iteration result is still acceptable here. But during the physical experiment two or three iterations may be needed due to a more complicated process. Furthermore, the application of the iterative adaptive approach also correct the focusing time at the same time, by which Baldock et al. (1996) found the focusing time is also dependent upon focusing amplitude and frequency bandwidth.

Conclusions
In this paper, an improved numerical model based on the high-order spectral method is employed to simulate the extreme wave in a two-dimensional flume. The focusing wave characteristics, in particular for the shift of focusing position and change of relative phasing are numerically in- Fig. 11. Comparison of phase between numerical and linear data after iteration for Case 2 with focusing amplitude A = 0.05 m.
LIU Dian-yong China Ocean Eng., 2021, Vol. 35, No. 1, P. 61-71 69 vestigated by this numerical model. It can be concluded that: (1) As we can see the focusing amplitude, frequency bandwidth as well as distribution of spectrum have an important effect on the wave groups evolution, contributing to the shift of focusing position downstream. It reveals that larger amplitude, narrower frequency bandwidth as well as CWA spectrum where the nonlinear wave-wave interactions dominate, and the relative phases of wave components change significantly, resulting in more shift.
(2) The iterative adaptive approach which can correct wave phases of each wave components is tested to be an acceptable approach for reproducing accurate wave focusing effectively. However, the proposed approach used in this study should still be verified by laboratory experiment in the future.