1. Introduction
Some studies have considered the frontal collision of two equal solitary waves both numerically and experimentally. Maxworthy [
1] conducted experiments on the frontal collision of two solitary waves and found that the maximum amplitude of these waves can reach twice the initial amplitude. Su and Mirie [
2] considered the frontal collision between two solitary waves on the surface of an inviscid uniform fluid and found that the waves generated by the collision maintained their original characteristics with third-order accuracy. Using a set of approximate equations derived by Su and Gardner [
3], Mirie and Su [
4] presented a numerical calculation of the collisions of two solitary waves and checked the phase shifts and maximum amplitude of the collision with a corresponding perturbation calculation. In addition, Akhmediev and Ankiewicz [
5] considered a linear coupler formed by two intersecting solitons. When the soliton amplitudes are equal, the transmission coefficient depends on the collision angle; however, this coefficient does not depend on the relative phase of the solitons. Cook et al. [
6] simulated the high-amplitude reflection of solitary waves on a vertical wall based on the boundary element method, and found that solitary waves can be converted into a series of dispersive waves after hitting the wall. Craig et al. [
7] proposed a numerical analysis based on the pseudo-spectral method to solve completely nonlinear equations and study the residual generated by the head-on collision of two solitary waves. Chamberel et al. [
8] numerically studied head-on collisions of two equal and two unequal steep solitary waves, and found that when the normalized amplitude
a0/
h of a solitary wave is greater than 0.60, a thin residual jet will appear. Recently, Slunyaev et al. [
9] studied strongly solitary wave groups of collinear surface waves using numerical simulations and laboratory experiments with Euler equations. The maximum displacement of a head-on collision is the sum of the two sets of peak amplitudes. Wave trains with modulational instability [
10,
11], which transform into envelope solitons, may lead to huge wave events, and the nonlinear interaction of solitons is a possible physical mechanism of freak wave (the wave height exceeds at least 2.2 times the significant wave height) generation [
12]. Interactions of solitary waves were investigated by simulating their collisions in a variety of wave contexts, such as plasma [
13], optical fibers [
14], and Bose–Einstein condensates [
15,
16]. Head-on collisions between two solitary waves in plasma have been studied through experiments [
17,
18] and numerical simulations using the particle-in-cell method [
19]. Zhang et al. [
19] found that the maximum amplitude during the collision process is less than the sum of the two amplitudes of the envelope solitary waves, and no phase shift was observed after the two envelope solitary waves collided head-on. The head-on collision of two envelope solitary waves was also studied analytically [
20].
Some new exact traveling wave solutions and solitary wave solutions for nonlinear partial differential equations have been obtained [
21,
22,
23]. These new solutions have many applications in other branches of physics and applied science. However, they are rarely used to systematically study the collision of solitary waves: the peak value of surface elevation,
ζmax, is affected by the relative water depth and wave steepness. For a long time, NLS equations have been used to simulate the propagation of waves in finite deep water. Benilov et al. [
24] studied the solitary wave solution of NLS with variable coefficients, which controls the weakly nonlinear envelope of surface gravity waves on topography. Grimshaw and Annenkov [
25] developed higher-order NLS equations with variable coefficients to describe how a water wave packet will deform and eventually be destroyed as it propagates shoreward from deep to shallow water. Rajan and Henderson [
26] considered a variable coefficient dissipative NLS equation, found a solution with uniform amplitude in an appropriate reference frame, and checked its linear stability. Grimshaw [
27] deduced a two-dimensional long wave equation with small but limited amplitude in variable depth; the differential equation describing the slow change in the parameters of the solitary wave was derived and solved when the solitary wave evolved from a uniform depth region. Benilov and Howlin [
28] studied the evolution of surface gravity wave packets in channels with terrain, where the terrain scale is far smaller than the limit of nonlinear/dispersion scale. Rajan et al. [
29] analyzed narrow-banded, periodic one-dimensional surface gravity carrier waves propagating on water of variable depth. Many contemporary researchers have carried out work on soliton collision but have not considered the role of water depth and wave steepness in soliton collision, which is the main purpose of this paper. In this study, the effects of surface tension were neglected; the analysis focused on gravity waves. It is well known that FSM is a powerful tool for solving partial differential equations (PDEs) both theoretically and numerically [
30,
31]. With the emergence of fast Fourier transform (FFT), Fourier spectral methods can provide a numerical discretization with so-called finite order convergence and high accuracy. In fact, for high-dimensional problems with periodic boundary conditions, the efficiency of the Fourier spectral method is equivalent to that of the finite difference method.
Based on the highly accurate FSM, we solved the NLS equation for gravity waves in finite water depth, as detailed in
Section 2. Results of the numerical studies on generation and interactions of the solitary groups are presented in
Section 3. The final discussion is reported in
Section 4.
2. Numerical Simulation for the NLS Equation
In dimensional variables, the following form is used:
where the complex wave envelope
A =
A(
x,
t);
x is the horizontal coordinate; and
t is time. The coefficients in the equation depending on the local water depth have the following definitions [
32]:
where
ω0 is the carrier wave angular frequency, which follows the following linear dispersion relationship:
where
Here,
k0,
h, and
c denote the wave number, water depth, and wave phase velocity, respectively.
cg denoting the group velocity can be written as follows:
The dimensionless variables are presented as follows [
33]:
where
ε0 = k0a0 is the wave steepness,
a0 is the free background amplitude parameter, and
γ0 is a scale factor that renders the computational domain in
ξ to 2π. The evolution of the complex envelope of surface elevation
A (
ξ,
η) can be described using the spatial NLS equation:
where
The dimensionless form of the focusing NLS equation is convenient to use:
which results from Equation (8) under the following transformation (
k0h > 1.363):
The localized solution of the NLS equation is the envelope soliton, which is presented in the following expression [
34]:
Therefore, Equation (8) yields the soliton solution in the following form:
The surface elevation
ζ (
x,
t) to the leading order is presented as follows:
A highly accurate Fourier spectral method was used to obtain numerical solutions for Equation (8). The spectral method is introduced by approximating the function as a sum of smooth basis functions [
30]:
where Φ
j(
ξ) are polynomials or trigonometric functions, and each aj must be determined.
In practice, the basis functions have many feasible options, such as Φ
j(
ξ) =
eikξ (the Fourier spectral method). The continuous Fourier transform of
A′(
ξ), which is transformed into Fourier space with regard to
ξ and its inverse Fourier transform, is defined as follows:
Similarly to the continuous Fourier transform, we define the discrete Fourier transform: let the solution interval be [0, 2π]; let
A′(
ξ) be the solution of Equation (8), then transform it into the discrete Fourier space as follows:
Using the inverse Fourier transform, we can obtain the following equations:
where
ξj = 2
πj/
N,
N is the discrete number, noting that
N is an even number. The spatial derivatives can be approximated using fast Fourier transform, whereas the time evolution is obtained using the fourth-order Runge–Kutta integrator [
31,
35]. Using the discrete Fourier transform and its inverse
F and
F−1, the Fourier transform operator
F can obtain a linear operator:
where
A′
(n)(
ξ) denotes the
n-th order derivative of
A′(
ξ). When
n = 2,
Using the Fourier transform, Equation (8) becomes
with
Under periodic boundary conditions on the spatial interval [0,
L],
L is the spatial domain. Then, the fourth-order Runge–Kutta scheme is used to solve the system of complex ordinary differential Equation (22). Letting
Equation (22) can be written in the following vector form:
where
F defines the right-hand side of (22). An initial condition is given to the abovementioned problem:
The classical fourth-order Runge–Kutta methods for the system (22) is presented as follows:
Figure 1 shows the profile of a soliton with
k0 = 4.8, ε
0 = 0.157,
γ0 = 0.1,
h = 1,
N = 2048, and step size ∆
η = 0.0006, 1 ≤
m ≤ 20,000, where
m denotes the
mth time step. We carried out sensitivity analyses on grid size and time step; the results are not dependent on the selection of grid size and time step. However, for simplicity, relevant images are not included here. In the transformed domain, the amplitude and shape of the soliton remain invariant during propagation. Firstly, the accuracy of the evolution envelope equation in the isolated subspace was investigated. Slunyaev et al. [
36] studied the stability of deep water waves and the existence of strong nonlinear groups through laboratory and numerical simulations.
Figure 2 shows a comparison between experimental and numerical simulated surface elevations at five different locations. The small wave steepness obtained in a physical basin may be due to the specificity of the boundary conditions of the wave generator (Equation (14)) or wave dissipation near the wave generator. Under the same initial conditions, compared with numerical simulation, the standing wave group in laboratory experiments ultimately appears to have a smaller steepness. Therefore, the experimental cases shown in
Figure 2 (
ε0 = 0.2) and
Figure 3 (
ε0 = 0.30) are compared with the numerical simulations of groups
ε0 = 0.157 and
ε0 = 0.25, respectively. Some deviations between the experimental and numerical results are observed at a relatively larger distance. NLS equation may be very useful in describing the general characteristics of initial narrow band water wave sequences that evolve in deep and intermediate-depth water. However, the focusing characteristics of the deep water NLS equation inevitably led to significant spectral broadening, which violates the basic assumptions used in the derivation of the NLS equation [
37], and can lead to deviations between numerical simulation and experimental results. The time frequency analysis shown in
Figure 4 and
Figure 5 clearly demonstrates the stationary waves.
In this study, the Morlet wavelet was selected as the mother wavelet. The wavelet spectrum was non-dimensionalized by the maximum wavelet energy at each location. Most of the wave energy was distributed near the center frequency of the carrier wave, and the instantaneous frequency of the intermediate wave packet in the wave train was a triangular signal. However, the measured wave heights covered higher-frequency components, and the energy of higher harmonics became significant. Therefore, the numerical simulation was consistent with the experimental measurement results. In addition, the NLS equation can be used to describe the general characteristics of initially narrow band water wave trains that evolve in deep and medium-deep water [
37].
3. Interaction of Solitons
Wave collision is a common phenomenon in water waves. The collision between two solitary waves manifests in two forms: overtaking collisions (solitary waves propagating in the same direction), and head-on collisions (solitary waves propagating in opposite directions). We selected the following initial conditions to study the interactions between two solitons [
5]:
where
b is the amplitude of each solitary wave, and
a1 and
a2 are the “velocities”. The values of
φ1 and
φ2 are the initial phases of the wave; when
φ1-
φ2 is 0 or
π, the solitary waves are initially in phase.
ξ1 and
ξ2 are the locations of solitons. We only considered the collision of the solitons with equal amplitude and equal velocity but with different directions of movement, and a wave steepness
ε0 =
k0A0 (
A0 = 0.01
b). A head-on collision between two solitary waves is depicted in
Figure 6,
Figure 7,
Figure 8,
Figure 9,
Figure 10 and
Figure 11. We observed that the amplitudes of the solitary waves were not altered after the interaction, and the solitary waves passed through each other without losing their identity but with a slight phase shift. During collision, the peak value of the amplitude was initially equal to the sum of the peak values of the two solitons. As shown in
Figure 7, the water depth, h, increased, whereas the other parameters were kept fixed, and the widths of the two solitary waves became smaller than those shown in
Figure 6. Therefore, the shallow-water region contains more waves than the deep water region, and no phase shift in solitary waves is observed in deep water (
Figure 7b). Thus, the propagation velocity of the solitary waves is constant.
Figure 8 and
Figure 9 show the surface elevation of the corresponding head-on collision between the two solitary waves shown in
Figure 6 and
Figure 7. With the increase in water depth, the distance that the solitons must travel before collision will decrease significantly (i.e., the evolution of the modulational instability is accelerated).
The effect of water depth,
k0h, on the surface elevation of a corresponding head-on collision between the solitary waves is presented in
Figure 10. During the collision, the value of
ε0 is initially equal to the sum of the peak values of the solitary waves, and the wave profile is based on the principle of linear superposition. However, with the increase in water depth,
k0h,
ζmax is less than the sum of the peak values of the two solitary waves, which is consistent with the results reported by Slunyaev et al. [
9] (the effect of finite depth reduces the amplitude of the eventual solitary group).
ζmax exhibits fluctuations with increasing water depths.
Figure 11 presents the peak value of surface elevation versus wave steepness,
ε0. In the case of smaller steepness,
ζmax is approximately equal to the sum of the peak values of the two solitary waves, and it exhibits significant fluctuations for
ε0 > 0.10. These results indicate that the water depth and wave steepness play important roles in the interaction of solitons. We also investigated the interaction of four solitary waves, and the initial condition is given as follows:
As shown in
Figure 12, head-on collisions occurred between the four solitary waves (
ε0 = 0.08,
k0 = 5,
h = 0.28,
a1 =
a2 = 21,
a3 =
a4 = 7,
b = 1.6,
ξ1 =
ξ2 = 6,
ξ3 =
ξ4 = 2,
φ1 =
φ2 =
φ3 =
φ4 = 0), where the four solitary waves retained invariant shapes, except for the phase shifts before and after the collision. During collision, the peak value of surface elevation,
ζmax, was initially equal to the sum of the peak values of the four solitary waves. In addition, the amplification factor was equal to six, which indicated a freak wave. Therefore, the four solitary waves interacted for a short time, and the superposition principle did not hold at the point of interaction in deep water.
Figure 13 shows the effect of water depth,
k0h, on the surface elevation of the corresponding head-on collision among the four solitary waves. During collision, the peak value of surface elevation,
ζmax, was initially equal to the sum of the peak values of the solitary waves; however,
ζmax exhibited significant fluctuations with increasing water depths.
4. Discussion
In this study, we modelled the solitary wave interactions by solving the NLS equation numerically. First, we solved the NLS equation to determine the spatial evolution of a solitary wave in finite water depth; the results showed good agreement between the numerical simulations and experimental measurements. However, the measured wave elevation covered evidently higher frequency components, and the energy of the higher harmonics became significant. Then, we considered collisions of the solitons with equal amplitudes and equal velocities, but with different directions of movement. During the collision, the solitary waves passed through each other and completely preserve their shapes and velocities, but a slight phase shift was observed. In contrast, no phase shift of the solitary waves was found in deep water. In addition, the peak value of surface elevation, ζmax, exhibited significant fluctuations with increasing water depths, k0h. In the case of relatively small steepness, ε0, ζmax was approximately equal to the sum of the peak values of the two interacting solitary waves, and it exhibited significant fluctuations for ε0 > 0.10 for the head-on collision. However, the mechanism of generation of this phenomenon remains an open problem. Similar results have been observed in studies of head-on collisions of four solitary waves. These results indicate that the water depth and wave steepness play important roles in soliton interactions, and the effects of the interactions of intense wave groups are important to study the mechanisms and manifestations of freak oceanic waves. These new results have many applications in the field of physics and other branches of applied science; such results show great potential for applications in the instability study of plasma and optical fibers. Meanwhile, FSM simulations provide a realistic description of the dynamics of nonlinear waves.