2.1. The HOS model

The HOS model is based on the potential flow formalism, where the fluid is assumed to be homogeneous, inviscid and incompressible, and the flow is assumed to be irrotational. The horizontal plane is defined by axis $\boldsymbol {x}=(x,y)$, and $z$ is the vertical axis pointing upwards. The flow is described by a velocity potential $\phi$, which satisfies the Laplace equation,

(2.1)\begin{equation} \Delta{\phi}=0. \end{equation}

The velocity potential on the free surface is defined by $\phi ^S(x,y,t)=\phi (x,y,\eta (x,y,t),t)$, such that the nonlinear free-surface boundary conditions are

(2.2)$$\begin{gather} {\eta_t}=(1+\boldsymbol{\nabla}\eta\boldsymbol{\cdot}\boldsymbol{\nabla}\eta)\phi_z-\boldsymbol{\nabla}\phi^S\boldsymbol{\cdot}\boldsymbol{\nabla}\eta, \quad {\rm on} \ z=\eta(x,y,t), \end{gather}$$

(2.3)$$\begin{gather}{\phi^S_t}=-g\eta+\frac{1}{2}(1+\boldsymbol{\nabla}\eta\boldsymbol{\cdot}\boldsymbol{\nabla}\eta)\phi_z^2-\frac{1}{2}\boldsymbol{\nabla}\phi^S\boldsymbol{\cdot}\boldsymbol{\nabla}\phi^S-\frac{p}{\rho}, \quad {\rm on} \ z=\eta(x,y,t), \end{gather}$$

where $\eta$ is the free-surface elevation, $g$ and $\rho$ are the gravitational acceleration and water density, respectively. Here $\boldsymbol {\nabla }=({\partial }/{\partial {x}},{\partial }/{\partial {y}})$. The moving pressure distribution $p$ is applied on the free surface, and it is introduced in § 2.2 in detail. Periodic boundary conditions are used in the horizontal plane.

The bottom boundary condition for a variable topography is as follows:

(2.4)\begin{equation} \boldsymbol{\nabla}{\phi}\boldsymbol{\cdot}\boldsymbol{\nabla}{\beta}-\phi_z=0, \quad {\rm on}\ z=-h_0+\beta(x,y), \end{equation}

where the total water depth is written as the sum of the average water depth $h_0$ and a depth variation $\beta$. A flat topography is adopted when $\beta =0$.

The detailed solution process for the velocity potential $\phi$ and the free surface elevation $\eta$ in HOS incorporating variable topography was introduced in Gouin, Ducrozet & Ferrant (Reference Gouin, Ducrozet and Ferrant2017). Here, a brief summary of the solution process is provided.

The velocity potential $\phi$ can be expressed as a power series expansion using components $\phi ^{(m)}$, and is truncated at order $M$, i.e. $\phi =\sum _{m=1}^M\phi ^{(m)}$ (see Dommermuth & Yue Reference Dommermuth and Yue1987; West et al. Reference West, Brueckner, Janda, Milder and Milton1987). Here $M$ is the order of free-surface nonlinearity. To account for the variable water depth, the velocity potentials at each order of $m$ are written as

(2.5)\begin{equation} \phi^{(m)}=\phi_{h_0}^{(m)}+\phi_{\beta}^{(m)}, \end{equation}

where $\phi _{h_0}^{(m)}$ is derived at the average depth $h_0$ and satisfies the boundary condition for a flat bottom; while $\phi _{\beta }^{(m)}$ generated by the variable topography can be further expanded with respect to $\xi =\beta /h_0$ under the assumption that $\xi \ll {1}$, and the expansion is truncated at order $M_b$ ($M_b$ is the order of bottom nonlinearity), i.e. $\phi _{\beta }^{(m)}=\sum _{l=1}^{M_b}\phi _{\beta }^{(m,l)}$. Here $\phi _{\beta }^{(m)}$ satisfies the boundary condition (2.4) for a variable topography and a Dirichlet condition on the still water surface,

(2.6)\begin{equation} \phi_{\beta}^{(m)}=0, \quad {\rm on} \ z=0. \end{equation}

The two potential terms $\phi _{h_0}^{(m)}$ and $\phi _{\beta }^{(m)}$ can be expressed by an eigenfunction expansion of free modes which satisfy the boundary conditions,

(2.7)$$\begin{gather} \phi_{h_0}^{(m)}=\sum_{p}\sum_{q}A_{pq}^{(m)}(t)\frac{\cosh\left[k_{pq}(z+h_0)\right]}{\cosh(k_{pq}h_0)}{\rm e}^{{\rm{i}}({k_{xp}x+k_{yq}y})}, \end{gather}$$

(2.8)$$\begin{gather}\phi_{\beta}^{(m)}=\sum_{p}\sum_{q}B_{pq}^{(m)}(t)\frac{\sinh(k_{pq}z)}{\cosh(k_{pq}h_0)}{\rm e}^{{\rm{i}}({k_{xp}x+k_{yq}y})}, \end{gather}$$

where $k_{xp}=2{\rm \pi} {p}/L_x$, $k_{yq}=2{\rm \pi} {q}/L_y$ and $k_{pq}=\sqrt {k_{xp}^2+k_{yq}^2}$. Here $L_x$ and $L_y$ are the lengths of computational domain in the $x$ and $y$ directions, respectively. The modal amplitude $A_{pq}^{(m)}(t)$ can be solved based on (2.2), (2.3) and (2.4) with $\beta =0$ by a fast-Fourier-transform method. By a Taylor expansion of the bottom boundary condition (2.4) around $z=-h_0$ truncated at order $M_b$, the modal amplitude $B_{pq}^{(m)}(t)$ can be computed as a function of $A_{pq}^{(m)}(t)$ (see Gouin et al. Reference Gouin, Ducrozet and Ferrant2017).

After the vertical velocity on the free surface $\phi _z(x,y,z=\eta )$ is computed based on (2.7) and (2.8), $\phi$ and $\eta$ can be time stepped based on (2.2) and (2.3). A fourth-order Runge–Kutta integrator is adopted for the time stepping. The sawtooth instabilities are mitigated using a Savitzky–Golay (SG) filter. The set-up of the SG filter is outlined in Appendix A.

The in-house HOS code, originally developed for a flat bottom based on the method by Dommermuth & Yue (Reference Dommermuth and Yue1987), has been optimized for computational efficiency by transferring input data to graphics processing units (GPU) and employing GPU-accelerated computing techniques. This code has been successfully applied to investigate phenomena such as freak waves in crossing sea states and four-wave resonant interactions in finite water depth (Liu et al. Reference Liu, Waseda, Yao and Zhang2022a; Liu, Waseda & Zhang Reference Liu, Waseda and Zhang2022b). Additionally, building upon the approach by Gouin et al. (Reference Gouin, Ducrozet and Ferrant2017), the in-house HOS code was adapted for a variable bottom. This updated version has been effectively utilized to study the impacts of free-surface and bottom nonlinearities on mini-tsunamis (Yao et al. Reference Yao, Bingham and Zhang2023).

2.2. Dynamic correction strategy

The travelling ship is modelled by a moving pressure distribution $p(x,y,t)$, which is shown in the last term of (2.3). A static pressure distribution can be computed directly from the desired hull shape $\eta _s(x,y,t)$,

(2.9)\begin{equation} p(x,y,t)=\rho{g}\eta_s(x,y,t)={\rho{g}df(x,t)q(y,t)}, \end{equation}

with

(2.10)\begin{gather} f(x,t)=\left\{ \begin{array}{@{}ll@{}} \cos^2\left[\dfrac{{\rm \pi}\left(\left|x-x^*(t)\right|-\frac{1}{2}\alpha_s{L}\right)}{(1-\alpha_s){L}}\right] & \frac{1}{2}\alpha_s{L}\leq{\left|x-x^*(t)\right|}\leq\frac{1}{2}{L},\\ 1 & \left|x-x^*(t)\right|\leq{\frac{1}{2}\alpha_s{L}}, \end{array} \right. \end{gather}

(2.11)\begin{gather} q(y,t)=\left\{ \begin{array}{@{}ll@{}} \cos^2\left[\dfrac{{\rm \pi}\left(\left|y-y^*(t)\right|-\frac{1}{2}\beta_s{B}\right)}{(1-\beta_s){B}}\right] & {{\frac{1}{2}\beta_s{B}\leq{\left|y-y^*(t)\right|}\leq\frac{1}{2}{B}}},\\ 1 & {{\left|y-y^*(t)\right|\leq{\frac{1}{2}\beta_s{B}}}}, \end{array} \right. \end{gather}

where $L$, $B$ and $d$ are the ship length, width and draft, respectively. The centre of the water plane of the ship is located at $(x^*,y^*)$, which varies as a function of time $t$. The two coefficients $\alpha _s$ and $\beta _s$ determine the shape of the wetted hull surface in the $x$ and $y$ directions, respectively (Ertekin et al. Reference Ertekin, Webster and Wehausen1986; Wu Reference Wu1987; Shi et al. Reference Shi, Malej, Smith and Kirby2018). If the ship length, width, draft and the static displacement $V_0$ are determined, the values of the two coefficients $\alpha _s$ and $\beta _s$ can be calculated based on $V_{0}=\iint {p}/(\rho {g})\,{\rm d}\kern0.06em x\,{\rm d} y$. The ship region is a rectangular region defined by the ship's length and width. The ship pressure is zero outside the ship region.

To include the dynamic effects of a moving ship, a dynamic correction for the pressure distribution is derived as follows. By introducing a local coordinate moving with the ship velocity $U$,

(2.12)\begin{equation} \left. \begin{gathered} x'=x-Ut,\\ y'=y,\\ z'=z, \end{gathered} \right\} \end{equation}

the dynamic boundary condition on the free surface (2.3) is written as

(2.13)\begin{equation} {\phi^S_t}=-g\eta+\frac{1}{2}(1+\boldsymbol{\nabla}\eta\boldsymbol{\cdot}\boldsymbol{\nabla}\eta)\phi_z^2-\frac{1}{2}\boldsymbol{\nabla}\phi^S\boldsymbol{\cdot}\boldsymbol{\nabla}\phi^S+U(t)\phi^S_x-\frac{p}{\rho}, \quad {\rm on} \ z=\eta(x,y,t).\end{equation}

The pressure distribution with dynamic correction, which is used to model a travelling ship, is set to

(2.14)\begin{align} p&=-\rho{g}\eta_s-\frac{\rho}{2}\left[\boldsymbol{\nabla}\phi^S\boldsymbol{\cdot}\boldsymbol{\nabla}\phi^S-\phi_z^2(1+\boldsymbol{\nabla}\eta\boldsymbol{\cdot}\boldsymbol{\nabla}\eta)\right]+{\rho}U(t){\phi^S_x +\rho{g}(\eta-\eta_s)} \end{align}

(2.15)\begin{align} &=-2\rho{g}\eta_s+\rho{g}\eta-\frac{\rho}{2}\left[\boldsymbol{\nabla}\phi^S\boldsymbol{\cdot}\boldsymbol{\nabla}\phi^S-\phi_z^2(1+\boldsymbol{\nabla}\eta\boldsymbol{\cdot}\boldsymbol{\nabla}\eta)\right]+{\rho}U(t){\phi^S_x}. \end{align}

The first part of the pressure in (2.14) simply cancels the dynamic pressure in (2.13), to ensure that $\phi ^S_t=0$ as long as $\eta =\eta _s$. The last term is a penalty term which is added to help enforce the condition that $\eta =\eta _s$, and it tends to improve the stability of the method. Further discussion of this penalty term is provided in Appendix B. We note that this is very similar to the scheme originally proposed by Lindberg et al. (Reference Lindberg, Glimberg, Bingham, Engsig-Karup and Schjeldahl2013).

For clarity, the static pressure distribution in (2.9) is referred to as ‘static pressure’, whereas the ship pressure with a dynamic correction in (2.15) is referred to as ‘dynamic pressure’ in subsequent discussions.

To compute the resistance of the ship, the pressure integration method, wave-cut analysis and far-field methods are adopted. These three methods are described in Appendix C.

3.1. Test cases for solitons above a flat bottom

To investigate the properties of ship-induced solitons over a flat topography, Ertekin et al. (Reference Ertekin, Webster and Wehausen1984) conducted a series of experiments in the Ship Model Towing Tank of the University of California, Berkeley. The dimensions of the tank are $L_x=61$ m by $W=2.44$ m. Some of these experimental cases are used in the present study to investigate the effects of a dynamic correction on ship-induced waves over a flat bottom. A Series 60 ship with a block coefficient $C_B=0.8$ travels along the centreline of the channel. The ship length, width and draft are $L=1.524$ m, $B=0.234$ m and $d=0.075$ m, respectively. The water depth over a flat topography is denoted by $h$. Two different water depths, $h=0.125$ m and 0.15 m are considered. The water depth to ship draft ratios are $h/d=1.67$ and 2, respectively. The Froude number $Fr=U/\sqrt {gh}$ ranges from 0.7 to 1.1. A narrow channel and a wide channel, whose widths satisfy $W/B=5.21$ and 20.85, respectively, are adopted. Because the ratio of the width of the experimental tank to the ship width is $W/B=10.43$, a false wall is placed along the centreline of the experimental tank to ensure the channel width satisfies $W/B=5.21$ in test cases with a narrow channel. In the test cases with a wide channel, a half-ship model is towed along the tank wall so that the channel width is equivalent to $W/B=20.85$. Sketches of the set-up for the experiments for both the narrow and wide channels are shown in figure 1. Wave gauges (WGs) are placed in the channel to record the upstream solitons. The specific locations of these gauges are marked in figure 1.

Figure 1.

Sketch of the experiments on solitons over a flat topography in the (a) narrow and (b) wide channels (Ertekin et al. Reference Ertekin, Webster and Wehausen1984). The ship moves along the centreline of the channel with a speed $U{(t)}$. The blue lines indicate the ship. Here $L$ and $B$ denote the length and width of the ship, respectively. Red circles with crosses at the centre represent the WGs.

In the numerical simulation, the ship, modelled by a pressure distribution, starts impulsively and then moves steadily along the $x$-axis with a constant speed. The set-up of the numerical simulation is illustrated in figure 2. In the numerical simulation, the channel width $W$ can be assigned any value, eliminating the need to simulate the false wall and the half-ship used in the experiment. The width of the computational domain is set to match the channel width, specifically $W/B=5.21$ with the number of nodes $N_y=128$ for the narrow channel and $W/B=20.85$ with $N_y=512$ for the wide channel. Because periodic horizontal boundary conditions are used in the HOS model and the ship travels along the centreline of the channel, the two boundaries along the $x$-axis can be regarded as the two channel walls. A detailed illustration about this can be found in Appendix D. The length of the computational domain should be sufficient to prevent waves propagating along the $x$-axis from exiting one boundary and re-entering through the other. The length of the computational domain in a non-dimensional form is $L_x/L=80.63$, with the number of nodes $N_x=8192$. The nonlinear order of the free surface is set to $M=3$, meaning that interactions up to third-order nonlinearities are included. To analyse the soliton characteristics and compare the numerical results with the experimental data, four numerical wave gauges (NWGs) are placed in the computational domain; NWGs 1, 2 and 3 are positioned identically to WGs 1, 2 and 3 in the experiments for both the narrow and wide channels; NWG4 is positioned at the same position as WG4 in the experiments with the narrow channel.

Figure 2.

Sketch of numerical simulation for the experiments conducted by Ertekin et al. (Reference Ertekin, Webster and Wehausen1984). The ship moves along the centreline of the channel at a speed $U{(t)}$. The blue lines indicate the ship. Here $L$ and $B$ denote the length and width of the ship, respectively. The four red circles with crosses at the centre indicate the positions of the four NWGs.

When the ship sails over a flat topography with a near-critical or supercritical speed ($Fr\geq 1$), both wake waves and upstream solitons are generated. The wake waves form a stationary wave pattern in the ship-fixed frame, and the resistance induced by the wake waves, denoted as $R_{ww}$, can be computed by the wave-cut analysis. As shown in figure 3, a TWC plane moving with the ship velocity $U$ is positioned aft of the stern. The distance between the wave-cut plane and the stern is set at three times the ship length. However, upstream solitons are periodically generated, and the resistance induced by upstream solitons $R_{uw}$ is computed by the far-field method. A far-field control plane, also moving at the velocity of the ship, is established ahead of the bow. The distance between this plane and the bow is configured to be five times the length of the ship.

Figure 3.

The positions of the wave-cut plane and the far-field plane used for calculation of wave resistance. Here $L$ is the ship length, and $U$ is the ship velocity.

3.2. Test cases for mini-tsunamis from a variable bottom

Mini-tsunamis are generated by a ship sailing from deep to shallow water. As shown in figure 4(a), the ship travels along the centreline of the channel which is restricted by the two sidewalls parallel to the sailing track of the ship. The coordinates of the centre of the water plane of the ship are $(x^*,y^*)$, where $x^*$ varies with time, but $y^*$ is constant. The ship initially accelerates from a state of rest and then proceeds to sail at a constant speed. The profile of the ship speed variation is defined by

(3.1)\begin{equation} U(t)= \left\{ \begin{array}{@{}ll} U_0\sin(t/T_a), & {{t/T_a\leq{\dfrac{\rm \pi}{2}}}},\\ U_0, & {{t/T_a>\dfrac{\rm \pi}{2}}}, \end{array} \right. \end{equation}

where $T_a$ is a coefficient related to the acceleration duration, and $U_0$ is the stable travelling velocity after the acceleration phase.

Figure 4.

Sketch of the test cases on mini-tsunamis over a variable topography. (a) Plan view: the ship moves along the centreline of the channel. The red square represents the centre of the water plane, whose coordinates are $(x^*,y^*)$. The ship accelerates in the relatively deep region, and moves at a constant velocity through a transition in water depth towards the shallower region. The regions between the two blue dashed lines are the transition zone for the water depth. Here $x_a$ and $x_b$ are the centre positions of the two depth changes. (b) Side view: $d$ is the ship draft. The initial position of the ship is denoted by $x_s$. Here $h_1$ and $h_2$ are the depths in the deep and shallow regions, respectively.

The bathymetry changes solely along the channel direction and remains uniform across the channel width. The bathymetric variation is shown in figure 4(b) and is described by $-h_0+\beta (x)$, where $h_0$ is the average water depth and

(3.2)\begin{equation} \beta(x)=-\tfrac{1}{2}\Delta{h}+\tfrac{1}{2}\Delta{h}\left[\tanh\left(\gamma\left(x-x_a\right)\right)-\tanh\left(\gamma(x-x_b)\right)\right], \end{equation}

where $\gamma$ determines the slope of the depth variation, $x_a$ and $x_b$ denote the centre positions of the two depth changes, respectively. The water depths in the deep and shallow regions are denoted by $h_1$ and $h_2$, respectively, with the depth difference given by $\Delta {h}=h_1-h_2$.

In the numerical simulations involving a variable bottom, the ratio of the ship length to ship width is $L/B=6$, and the ratio of ship width to ship draft is $B/d=6$. The two coefficients governing the shape of the wetted hull surface in (2.10) and (2.11) are $\alpha _s=0.5$ and $\beta _s=0.25$, respectively. In the relatively deep region, the water depth is $h_1/d=11.64$, while in the relatively shallow region, the water depth is $h_2/d=4.36$. Therefore, the average water depth is $h_0/d=8$. The size of the computational domain is $(L_x, L_y)=(300h_0,3h_0)$, and the number of nodes is $(N_x, N_y)=(8192,512)$. Because the width of the computational domain is equal to the channel width, the ratio of the channel width to ship width is $W/B=4$. The centre positions of the two depth changes in (3.2) are $x_a=90h_0$ and $x_b=210h_0$, respectively. The initial position of the ship is $x_s=30h_0$, and the coefficient controlling the acceleration duration is $T_a/\sqrt {h_0/g}=90$. The ship finishes the acceleration before reaching the centre position of the depth change $x_a$, and proceeds through the depth change with a constant speed. Wave characteristics and wave resistance are studied across multiple scenarios, wherein the Froude number based on the average water depth $Fr_0=U/\sqrt {gh_0}$, and the slope of the depth variation $\gamma$ are varied. All of the test cases featuring variable topography are outlined in table 1.

Table 1.

Main parameters used in the variable bathymetry test cases. Here $Fr_0$ denotes the Froude number based on the average depth and $\gamma$ controls the slope of the depth variation.

It should be noted that if the order of bottom nonlinearity $M_b=0$, a flat bottom corresponding to the average water depth $h_0$ is used. If $M_b\geq 1$, the bathymetric variation $\beta (x,y)$ is considered, and the corresponding water depth is expressed by $-h_0+\beta (x,y)$. The increase of $M_b$ does not affect the variable bathymetry but decreases the truncation error of the velocity potential, i.e. the numerical results have an exponential convergence rate with respect to $M_b$ (Gouin et al. Reference Gouin, Ducrozet and Ferrant2017; Yao et al. Reference Yao, Bingham and Zhang2023). In all the test cases with a variable bottom in the present study, the amplitude and wavelength of mini-tsunamis increase with $M_b$. As $M_b$ increases from 2 to 3, the increase in the wave amplitude is less than 6.3 %. The increase is less than 2 % as $M_b$ increases from 3 to 5 and is less than 0.2 % as $M_b$ increases from 5 to 6. The increase in the wavelength is less than 0.3 % as $M_b$ increases from 2 to 5. Therefore, within the HOS model, the order of bottom nonlinearity is established as $M_b=5$, which is enough to consider the impact of bottom variation on waves. The numerical results also have an exponential convergence rate with respect to the order of free-surface nonlinearity $M$. In all the test cases, the increase in the wave amplitude can reach 20 % as $M$ increases from 2 to 3, but is less than 3.9 % as $M$ increases from 3 to 4. The increase in the wavelength is less than 5.8 % as $M$ increases from 2 to 3, and is less than 0.75 % as $M$ increases from 3 to 4. Therefore, $M=3$ is chosen within the HOS model.

4.1. Validation of the wave resistance calculation

To estimate the wave resistance of an air-cushion vehicle (ACV), Doctors (Reference Doctors1993) introduced a solution based on potential-flow theory and compared their results with model test experiments. The ACV is represented by a moving static pressure distribution,

(4.1)\begin{align} p_{ACV}&= \frac{p_0}{4}\times\left\{ \tanh\left[\alpha_a\left( \frac{x-x^*}{L}+\frac{1}{2} \right)\right]-\tanh\left[\alpha_a\left( \frac{x-x^*}{L}-\frac{1}{2} \right)\right] \right\}\nonumber\\ &\quad \times\left\{ \tanh\left[\beta_a\left( \frac{y-y^*}{L}+\frac{1}{2} \right)\right]-\tanh\left[\beta_a\left( \frac{y-y^*}{L}-\frac{1}{2} \right)\right] \right\}, \end{align}

where $\alpha _a=24$, $\beta _a=24$ and $p_0=0.0127\rho {g}L$. The ship has $L/B=3/2$ and travels on a flat bottom. The water depth satisfies $h/L=3.58$.

The wave resistance $R_w$ of the ACV is computed at various Froude numbers based on the ship length $Fr_L=U/\sqrt {gL}$ within the range 0.4–1.14. When the ship moves at constant speed, only wake waves are generated, no upstream waves are present, so the wave-cut analysis method is used to compute the wave resistance $R_w$. A TWC plane moving with the ship velocity is positioned after the stern at a distance of $3L$. As shown in figure 5, the wave resistance coefficients $C_w=\rho {g}R_w/p_0^2B$ computed by the pressure integration and wave-cut analysis are compared with the experimental data. The results from the two numerical methods exhibit strong agreement with each other, and they both align with the experimental results.

Figure 5.

Wave resistance coefficient $C_w=\rho {g}R_w/p_0^2B$ for an ACV as a function of Froude number $Fr_L=U/\sqrt {gL}$. The ACV has $L/B=3/2$ and travels in open water with $h/L=3.58$. The red circles indicate pressure integration results while the blue crosses show the wave-cut analysis results. The lines are splines through the point values. The black triangles show the experimental data in Doctors (Reference Doctors1993).

4.2. Validation of the dynamic correction strategy

Several experimental cases in Ertekin et al. (Reference Ertekin, Webster and Wehausen1984) are chosen to validate the dynamic correction strategy in the HOS model. The experimental set-up and numerical model configuration have been detailed in § 3.1. Figure 6 shows the data from WG4 in different test cases where the water depth is $h/d=1.67$ and the Froude numbers are 0.8, 1.0 and 1.1, respectively. The outcomes from the NWG4, represented by the blue solid lines, are derived using the dynamic pressure distribution. In each plot, the peak of the first main wave is aligned to ${t'/\sqrt {h/g}}=0$. Behind the main wave, several short waves propagate with a mean free-surface level above zero. The amplitude of the upstream solitons and the mean free-surface level rise with the Froude number. The data from the NWG4 correspond reasonably well to the experimental data, although there are some phase differences observed in the short waves behind the main wave, particularly when the Froude number $Fr$ is large. The phase differences could be attributed to the idealization of the hull shape. The bow and stern shapes described by (2.10) and (2.11) are not identical to the actual bow and stern of the ship model used in the experiments.

Figure 6.

Comparison of the wave gauge records computed by the HOS ($M=3$) with the dynamic correction strategy and those from experiments. The water depth satisfies $h/d=1.67$. The Froude numbers $Fr=U/\sqrt {gh}$ are (a) $Fr=0.8$, (b) $Fr=1.0$ and (c) $Fr=1.1$. The blue solid lines are computed results and the black circles are the experimental data from Ertekin et al. (Reference Ertekin, Webster and Wehausen1984).

As shown in figure 7, the amplitude of the first upstream wave generated by the dynamic pressure distribution in the water depth $h/d=2$ is compared with the experimental records. The Froude number varies from 0.7 to 1.1. The numerical results agree well with the experimental records, providing confidence that the dynamic correction strategy in the HOS model accurately captures the underlying physics, despite the geometric simplification of the hull form.

Figure 7.

The variation of the non-dimensional amplitude $A/h$ of the main upstream wave with the Froude number. The water depth is $h/d=2$. The black circles represent the experimental records from Ertekin et al. (Reference Ertekin, Webster and Wehausen1984). The blue squares represent the numerical results computed by HOS ($M=3$) with the dynamic correction strategy.

5.1. Wave patterns

To examine the impact of dynamic correction on ship-induced waves over a flat topography, the upstream waves and wake waves generated by the static and dynamic pressure distributions are compared. The ship travels along the centreline of a wide channel ($W/B=20.85$). The water depth is $h/d=1.67$. Two test cases with different Froude numbers, $Fr=1$ and 1.1, are used. The detailed set-up can be found in § 3.1.

Once the ship starts to move with $Fr=1$, the surface elevation along the centreline of the ship region is shown in figure 8. Since the ship is stationary at the initial time $t=0$, the hull shapes modelled by the dynamic and static pressure distributions are identical, and they correspond to the desired hull shape. After the ship moves impulsively, the ship bottom modelled by the static pressure shows notable oscillations, with the mean level of the ship bottom exceeding the initial ship draft. The volume of the ship modelled by the static pressure is increased by approximately 28 % during the initial acceleration of the ship, as shown in figure 9. On the other hand, the ship bottom modelled by the dynamic pressure stays flat and the ship volume is very nearly constant with time. The dynamic correction strategy accurately preserves the desired hull shape both while accelerating and when moving at constant speed.

Figure 8.

The surface elevation $\eta$ along the centreline of the ship region at the initial stage of ship movement. The initial position of the ship is at $(x-x_s)/L=0$. Here $d$ denotes the ship draft. The ships are modelled by (a) the dynamic pressure distribution and (b) the static pressure distribution, respectively.

Figure 9.

Variation of the ship volume $V$ at the initial stage of ship movement. The initial position of the ship is at $(x-x_s)/L=0$. Here $V_0$ is the original ship volume.

As shown in figure 10, the data of the NWG3 computed by the two types of moving pressure distribution are compared. The WG3 is placed near the channel wall, and its location is indicated in figure 2. To access the accuracy of the waves induced by the two types of moving pressure distribution, the experimental data of WG3 are also shown in figure 10. Significant depressions are observed at $195\leq {t/\sqrt {h/g}}\leq 221$ in figure 10(a) and $186\leq {t/\sqrt {h/g}}\leq 204$ in figure 10(b), respectively. The depression emerges along the ship body, and propagates towards the channel wall. The lengths of the depressions generated by the dynamic pressure align closely with those observed in the experiments, albeit being shorter than those induced by the static pressure. The amplitude of depression in the numerical results is smaller than that observed in the experiments. This discrepancy can be attributed to the idealization of the hull shape. Due the difference in the length of the depression, the phase of the wake waves generated by the static pressure is very different from that in the experiments. However, the length of the depression and the phase of the wake waves induced by the dynamic pressure are much closer to the experimental observations. The amplitudes of the wake waves generated by the static pressure are larger than those induced by the dynamic pressure and those in the experiments, particularly for the waves between $266\leq {t/\sqrt {h/g}}\leq 283$ in figure 10(a). The upstream solitons induced by the two types of pressure distribution and those in the experiments are comparable.

Figure 10.

Time series of free surface elevation $\eta$ recorded by WG3. The whole gauge records comprise the upstream waves, the depression and the wake waves. The black dash–dotted line represents the experimental data from Ertekin et al. (Reference Ertekin, Webster and Wehausen1984). The blue and red dashed lines denote the results computed by the dynamic pressure and static pressure, respectively. Here (a) $Fr=1.0$, the depression is between $195\leq {t/\sqrt {h/g}}\leq 221$; (b) $Fr=1.1$, the depression is between $186\leq {t/\sqrt {h/g}}\leq 204$.

To quantify the disparities between the numerical results and the experimental records, a correlation coefficient $\epsilon$ is employed. The correlation coefficient is defined as

(5.1)\begin{equation} \epsilon[\eta_{num}(t),\eta_{exp}(t)]=\frac{\int_t \left[\left(\eta_{num}(t)-\overline{\eta_{num}}\right)\left(\eta_{exp}(t)-\overline{\eta_{exp}}\right)\right] {\rm d}t} {\left\{\left[ \int_t \left(\eta_{num}(t)-\overline{\eta_{num}}\right)^2 \,{\rm d}t\right] \left[\int_t \left(\eta_{exp}(t)-\overline{\eta_{exp}}\right)^2 {\rm d}t \right]\right\}^{1/2}}, \end{equation}

where $\eta _{num}(t)$ is the time series of free surface elevation computed by the numerical methods, $\eta _{exp}(t)$ is that recorded in the experiments. The time-averaged values of $\eta _{num}(t)$ and $\eta _{exp}(t)$ are denoted by $\overline {\eta _{num}}$ and $\overline {\eta _{exp}}$, respectively. The complete gauge records shown in figure 10 encompass the upstream waves, the depression and the wake waves. In the two test cases, the depressions are between $195\leq {t/\sqrt {h/g}}\leq 221$ and $186\leq {t/\sqrt {h/g}}\leq 204$, respectively; the upstream waves are between $169\leq {t/\sqrt {h/g}}<195$ and $160\leq t/ \sqrt {h/g}<186$, respectively; the wake waves are between $221<{t/\sqrt {h/g}}\leq 292$ and $204<{t/\sqrt {h/g}}\leq 292$, respectively. The values of correlation coefficient $\epsilon$ for upstream waves, wake waves and the whole gauge records are shown in table 2. The correlation coefficients of upstream solitary waves exceed 0.94 for both methods, indicating that the upstream solitary waves generated by both methods closely correspond to the experimental observations. However, the wake waves are more sensitive to the dynamic correction. Due to the phase differences and larger wave amplitudes, the correlation coefficient of the wakes waves for the static pressure is lower than 0.4, whereas for the dynamic pressure it exceeds 0.6. The dynamic correction strategy enhances the correlation coefficient for wake waves by approximately 0.25. Moreover, the correlation coefficient of the entire gauge records for the dynamic pressure is also higher than that for the static pressure. Therefore, the dynamic correction strategy improves the accuracy in simulating the ship-induced waves.

Table 2.

Correlation coefficients $\epsilon$ between the numerical and experimental results shown in figure 10. The correlation coefficients of upstream waves, wake waves and the whole gauge records for the dynamic and static pressure distributions are calculated, respectively.

5.2. Wave resistance

Periodic upstream solitons and wake waves are generated at near-critical and supercritical conditions ($Fr\geq 1$). The free surface along the centreline when solitons are generated at the bow is shown in figure 11, where the ship travels along the centreline of a narrow channel with a width of $W/B=5.21$, in relative water depth of $h/d=2$, at $Fr=1.0$. The blue and red lines show the dynamic and the static pressure results, respectively. The region between $11.55\leq (x-x_s)/L\leq 12.55$ is the ship region where the ship modelled by the dynamic pressure distribution has a flat bottom, while oscillations are shown on the bottom of the ship modelled by the static pressure distribution. Therefore, with the dynamic correction strategy, the described hull shape is preserved while the ship is moving. At around $(x-x_s)/L=13.56$, the amplitudes of upstream solitons induced by the dynamic and static pressure distributions are comparable. However, differences in both wave amplitude and phase are shown in the wake waves exhibited at $(x-x_s)/L\leq 11.55$. The wake waves induced by the dynamic pressure have a smaller amplitudes than those generated by the static pressure.

Figure 11.

The free surface elevation $\eta$ along the centreline of the channel with a width of $W/B=5.21$ at $t=133\sqrt {h/g}$. The water depth satisfies $h/d=2$.

The resistance induced by the wake waves, denoted as $R_{ww}$, and the resistance induced by the periodic upstream solitons, denoted as $R_{uw}$, are computed by the wave-cut analysis and the far-field method, respectively. The positions of the TWC plane and the far-field control plane are given in § 3.1. The total wave resistance can be calculated from $R_w=R_{ww}+R_{uw}$. Correspondingly, the resistance coefficients are expressed by $C_w=\rho {g}R_w/p_0^2B$, $C_{ww}=\rho {g}R_{ww}/p_0^2B$ and $C_{uw}=\rho {g}R_{uw}/p_0^2B$, respectively. The symbol $p_0=\rho {g}d$ where $d$ is the ship draft shown in (2.9).

Because of the periodic generation of upstream solitons, the variation of resistance computed on the far-field control plane over time is also periodic, as shown in figure 12. Here $C'_{uw} (t)$ is the instantaneous resistance coefficient generated by the upstream solitons. The period $T_s$ of the solitons is $44\sqrt {h/g}$, and the discrepancy in resistance over different periods is relatively small. The time-averaged value of $C'_{uw} (t)$ over one period is used as the resistance coefficient induced by the upstream waves, denoted as $C_{uw}$.

Figure 12.

The variation of the instantaneous resistance coefficient generated by the upstream solitons $C'_{uw} (t)$ with time. Here $T_s$ is the period of the solitons. The channel width and the water depth satisfy $W/B=5.21$ and $h/d=2$, respectively.

The fluctuation of wave resistances induced by the dynamic and static pressure distributions over one soliton period is compared in figure 13, and the resistance coefficients obtained by the different methods are summarized in table 3. The wake waves form a stable wave pattern in the ship-fixed reference frame, such that the resistance corresponding to wake waves remains relatively constant over time. Due to the larger amplitudes induced by the static pressure distribution, the resistance coefficient $C_{ww}$ is larger than that for the dynamic pressure distribution. Additionally, because the upstream waves generated by the two types of pressure are comparable, the values of $C_{uw}$ are also nearly identical. Therefore, the total wave resistance induced by the dynamic pressure distribution is smaller than that for the static pressure distribution. To validate our results, the total resistance coefficient $C_w$ computed by the pressure integration (PI) method and from experimental (EXP) records are also presented in figure 13 and in table 3. The periodic variation of the instantaneous resistance coefficient $C'_{uw} (t)$ implies that the instantaneous coefficient of the total resistance ${C'_{uw} (t) + C_{ww}}$ also exhibits periodic behaviour. The total resistance coefficient calculated by $C_w=\overline {C'_{uw} (t)+C_{ww}}$ corresponds well with that obtained by the pressure integration method. However, the total wave resistance induced by the dynamic pressure distribution aligns well with the experimental results. This comparison illustrates that a more accurate wave resistance can be obtained through the dynamic correction strategy.

Figure 13.

Variation of the wave resistance coefficients over one soliton period $T_s$. The ships are modelled by (a) the dynamic pressure distribution and (b) the static pressure distribution, respectively. The blue solid line represents the resistance coefficient from the wake waves $C_{ww}$. The red solid line represents the coefficient of the instantaneous resistance induced by the upstream waves $C'_{uw} (t)$, and its time-averaged value in one period $C_{uw}=\overline {C'_{uw} (t)}$ is denoted by the red dashed line. The coefficient of the instantaneous total wave resistance computed by $C'_{uw} (t)+C_{ww}$ is denoted by the black solid line, its time-averaged value in one period $C_{w}=\overline {C'_{uw} (t)+C_{ww}}$ is denoted by the black dashed line. The magenta dashed line and the green dash–dotted line represent the total wave resistance coefficients $C_w$ obtained by the pressure integration (PI) method and the experiments (EXP), respectively.

Table 3.

Comparison of the resistance coefficients computed by the numerical models and those in the experiments. Here $C_{ww}$ and $C_{uw}$ are the coefficients of resistance induced by the wake waves and upstream waves, respectively. Here $C_w$ is the total wave resistance coefficient. Calculation methods are indicated in parentheses. Here ‘WA’, ‘FM’, ‘PI’ and ‘EXP’ denote ‘wave-cut analysis’, ‘far-field method’, ‘pressure integration’ and ‘experiment’, respectively.

6.1. Generation of mini-tsunamis

Figure 14 shows an example of the mini-tsunami generation process using the conditions outlined in Test 3 from table 1. The ship, modelled by a dynamic pressure distribution, travels at $Fr_0=0.5$. Due to the narrowness of the channel, i.e. the ratio of channel width to ship width $W/B=4$, reflection of wake waves is significant. Figure 14(a) shows a time before the ship reaches the depth change centred at $x_a$, so the upstream waves are negligible. As the ship passes the depth variation, mini-tsunamis are induced at the bow (shown in figure 14b), and then propagate ahead of the ship (shown in figure 14c).

Figure 14.

The surface elevation $\eta$ at different times in Test 3: (a) $t/\sqrt {h_0/g}=139$; (b) $t/\sqrt {h_0/g}=186$; (c) $t/\sqrt {h_0/g}=224$. Here $h_0$ is the average water depth, and $x_a$ is the centre of the depth change.

Figure 15 shows the free surface along the centreline of the ship region while the ship is sailing through the depth variation. The centre of the depth variation is at $(x-x_a)/L=0$. The ships modelled by the dynamic and static pressure distributions are shown in figures 15(a) and 15(b), respectively. As the ship travels from the deep water to the shallow water, the bottom of the ship modelled by the dynamic pressure stays flat, and the ship volume at different places is nearly the same as the original volume $V_0$, which is shown in figure 16. However, substantial oscillation is shown at the bottom of the ship modelled by the static pressure. Because of the deformation of the hull shape, the mean level of the ship bottom is much deeper than the original ship draft which is at $\eta /d=-1$, inducing a change in the ship volume, which is denoted by the blue line and squares in figure 16. Due to the deformation of the hull shape during the initial acceleration phase, the ship volume is around $V=1.09V_0$ in the relatively deep region. As the ship travels into the shallow-water region, additional deformation occurs, and the ship volume is increased by approximately 27 % compared with that in the relatively deep water. After the ship has passed the depth change and travels in shallow water, the three-dimensional free surface of the ship region modelled by the two types of pressure are compared in figure 17. From these comparisons, the dynamic correction strategy is able to preserve the described hull shape and ensures that the ship has a flat bottom and a stable ship volume.

Figure 15.

Surface elevation $\eta$ along the centreline of the ship region as the ship travels through the depth variation. The centre of the depth change is at $(x-x_a)/L=0$. Here $d$ denotes the ship draft. The ships are modelled by (a) the dynamic pressure distribution and (b) the static pressure distribution, respectively.

Figure 16.

Variation of the ship volume $V$ as the ship travels over the depth variation. The centre of the depth change is at $(x-x_a)/L=0$. Here $V_0$ is the original ship volume.

Figure 17.

Three-dimensional surface $\eta$ in the ship region after the ship has passed the depth change. Here $x_a$ represent the centre of the depth change. (a) Dynamic pressure distribution and (b) static pressure distribution.

During the generation process of the mini-tsunamis, the free surface elevation along the centreline of the channel is shown in figure 18. The wave amplitude of the mini-tsunamis induced by the static pressure distribution is approximately $0.016h_0$. Because the static pressure distribution has an uneven bottom and a larger ship volume than that of the dynamic pressure distribution, the wave amplitude of the mini-tsunamis induced by the dynamic pressure distribution is smaller than that induced by the static pressure distribution. To be specific, the wave amplitude computed by the static pressure distribution is 45 % larger than that obtained using the dynamic pressure distribution. This illustrates that the dynamic correction has a significant impact on the wave amplitude of mini-tsunamis.

Figure 18.

Surface elevation $\eta$ along the centreline of the channel as the mini-tsunamis are generated from the bow: (a) $t/\sqrt {h_0/g}=182$; (b) $t/\sqrt {h_0/g}=190$; (c) $t/\sqrt {h_0/g}=207$.

6.2. Mini-tsunami wave characteristics

To systematically study the impact of the dynamic correction on the wave amplitude $A$ and wavelength at half-crest height $\lambda _{1/2}$ of mini-tsunamis, the test cases from table 1 are considered.

Figure 19 compares the generated wave amplitudes from the two methods as a function of increasing Froude number and bottom slope parameter. At $Fr_0=0.40$, the difference in wave amplitudes induced by the two methods is small, and both results are around $A/h_0=0.005$. As the Froude number increases, the wave amplitude increases, as does the induced fluid acceleration around the ship, and the dynamic effects become increasingly significant. For example, the dynamic correction strategy decreases the wave amplitude by approximately 42 % at $Fr_0=0.55$. This indicates the importance of a dynamic correction when the ship speed is large. Moreover, based on the results for the dynamic pressure and static pressure, fitted curves corresponding to $Fr_0^\alpha$ are obtained, and are denoted by the blue and red lines in figure 19(a), respectively. The value of $\alpha$ is 3.9 for the dynamic pressure distribution and 5.0 for the static pressure distribution. As shown in figure 19(b), the wave amplitude also increases with the slope of the depth variation $\gamma$. However, the decrease in wave amplitude due to the dynamic correction is nearly constant with a value of approximately $0.005h_0$ as the bottom slope changes from 0.7 to 10.15.

Figure 19.

Influence of dynamic correction on the amplitudes $A$ of mini-tsunamis in cases with different (a) Froude numbers based on the average water depth $Fr_0$ and (b) slopes of the depth variation $\gamma$.

Figure 20 shows the wavelength at half-crest height $\lambda _{1/2}$ of mini-tsunamis generated by the dynamic and static pressure distributions. The wavelength diminishes with increasing Froude number and remains nearly constant with varying bottom slope. The dynamic correction results in a slight reduction in wavelength, but the effect is not sensitive to either Froude number or bottom slope. The decrease is less than 6.9 % in all the tested cases.

Figure 20.

Influence of dynamic correction on the wavelength at half-crest height $\lambda _{1/2}$ of mini-tsunamis in cases with different (a) Froude numbers based on the average water depth $Fr_0$ and (b) slopes of the depth variation $\gamma$.