The resistance of a trans-critically accelerating ship in shallow water

ABSTRACT

The acceleration resistance of a vessel advancing in shallow water is investigated. Four acceleration intensities and two water depths are modelled using the CFD and potential flow methods. The results show a pronounced peak in resistance exists near the critical depth Froude number, but its location and magnitude are sensitive to the acceleration intensity and water depth. Excellent agreement between the results obtained from the CFD and potential flow methods is found in the low and high depth Froude number ranges regardless of water depth or acceleration, indicating that linear and unsteady methods can provide robust results at a low cost in those ranges. The magnitude of the resistance peak and its position are sensitive to nonlinear effects, evidenced by slight disagreements between the two adopted methodologies. The variation in the results produced by the two solvers is found to be sensitive to the parameters investigated.

KEYWORDS:

• Ship resistance
• URANS
• potential flow
• shallow water
• acceleration

1. Introduction

When sailing in shallow water, vessels are affected by a hydrodynamic interaction between the hull and seabed. This interaction changes the wave pattern and increases ship resistance. The former is of particular interest, because the so-called Kelvin wake depends on the combination of water depth and ship speed, expressed by the depth Froude number, $F_h=V/\sqrt{gh}$, where $V$ is the ship speed, $h$ is the water depth and $g$ is the gravitational acceleration. In very shallow water, the denominator of $F_h$ is equal to the wave speed. The depth Froude number is therefore the ratio of the ship and wave speeds.

If the critical value of the depth Froude number ($F_h=$1) is surpassed for example when a vessel advances over a change in the water depth shockwaves may be generated (Jiang et al. 2002). Under the aforementioned conditions, the vessel radiates large amplitude waves (Grue 2017). Upon passing over a reduction in water depth, the depth Froude number changes from a subcritical value ($F_h<1$) to a supercritical value ($F_h>1$). In shallow water hydrodynamics the water depth and/or the vessel speed may be manipulated to the same effect, therefore, one may expect that similar phenomena can be reproduced by either changing the speed or changing the depth.

The main aim of the present paper is to investigate the phenomena occurring as a ship accelerates past the critical depth Froude number for a constant water depth. The present study combines Unsteady Reynolds Averaged Navier-Stokes results with 3D panel-based potential flow solutions using MHydro (Li et al. 2019) of trans-critical ship performance in shallow water by varying the speed of the parabolic hull to address the above aim.

The remainder of this paper proceeds with a literature review where the main contributions in the field are investigated. Following this, case studies and details of the numerical solver and potential flow method are given. Next, results focusing on the resistance characteristics are discussed, followed by conclusions.

2. Background

The problem of predicting ship resistance has historically been treated as a steady. That is, a single resistance value is obtained for each condition, for example, for each combination of speed and water depth. There are however several contributions to the study of unsteady ship resistance estimation. To the best of the authors’ knowledge, Havelock (1949) was first to investigate the transient wave resistance on a submerged body. Havelock (1949) focused on a cylinder with uniform acceleration, expressing the wave resistance as the sum of a component at constant velocity and a component due to the effect of the acceleration. Following Havelock (1949), Lunde (1957, 1951) developed a mathematical theory to predict transient wave resistance. Later, Wehausen (1964, 1961) proposed an asymptotic formulation expressing the unsteady wave resistance whose oscillations are predicted to decay inversely with time. Calisal (1977) extended Wehausen's (1964) method to include the development of the wave pattern in unrestricted water. Pinkster (2004) developed a potential flow solver to investigate the ship-to-ship problem, where the unsteady hydrodynamic force can be expressed with the aid of an added mass term. More recently, Day et al. (2009) built on Wehausen’s (1964) work to include the effect of finite depth as well as width in a towing tank, which Wehausen (1964) considered infinite.

Experimental and theoretical investigations on unsteady resistance are primarily motivated by the performance of unconventional vessels such as hovercraft (Barratt 1965; Doctors and Sharma 1972; Yeung 1975; Haussling and van Eseltine 1978). These craft have historically been simplified as two-dimensional pressure distributions on the water surface (Doctors 1993, 1975). The aforementioned studies were motivated by the presence of oscillations in the time-history signal of craft towed in a tank. The key problem in such cases is that the oscillation amplitude in the time-history of the resistance can be significant and may complicate the accurate measurement of the steady value (Li et al. 2019). Although a multitude of unsteady wave drag theories have been proposed, as shown above, the subject of a ship accelerating past a critical depth Froude number is rarely investigated.

In shallow waters, the Kelvin wake and corresponding to the angle between each arm of the V-shaped wave pattern and centreline ($\theta$), and depends on the depth Froude number (Lee and Lee 2019). According to theory, when $F_h=1$ the wave pattern is transformed into a wave front that is perpendicular to the ship centreline (Havelock 1908). That is, the Kelvin half-angle is equal to 90° (Tunaley 2014) as shown in Figure 1.

Figure 1. Kelvin half angle as a function of the depth Froude number. The relationships used to construct this figure are given in Havelock (Havelock 1949). Experimental measurements of the Kelvin half angle can be found in Johnson (1957).

Many approaches to modelling ship drag theoretically are unable to re-create the phenomena occurring at $F_h=1$. For example, the well-known slender body theory (Tuck 1967, 1966) predicts a singularity at that speed. Lea and Feldman (1972) and later, Gourlay and Tuck (2001) addressed this gap, but the methodology is valid for a ship travelling at a uniform speed. There are also thought to be non-linear and dispersive effects one must handle at and around the critical depth Froude number. In addition, one must consider unsteady effects when a ship accelerates. The present study combines a linear unsteady potential flow solver with an unsteady Reynolds-Averaged Navier-Stokes (URANS) solver to establish the relative importance of the nonlinear and unsteady effects. Specifically, if unsteady effects dominate resistance with little contribution of nonlinear effects, then the solutions from the two aforementioned solvers must agree. If nonlinear effects are important, one expects to observe some disagreement between the two adopted methodologies.

To the best of the authors’ knowledge, Kevorkian and Yu (1989) who were first to focus exclusively on trans-critical transitions caused by changes in the water depth. Redekopp and You (1995) modelled a disturbance passing through the critical depth Froude number using the forced Korteweg-deVries (fKdV) equation. The above approaches model the problem in a single dimension. Jiang et al. (2002) combined a slender body approach with Boussinesq equations to model the two-dimensional free surface as a vessel advances over a restriction in the water depth. Later, Torsvik et al. (2006) used a set of forced Boussinesq equations to study the trans-critical transition of a pressure distribution on the water surface. The acceleration intensity was shown to govern whether a solitary wave is emitted from the vessel, with relatively low acceleration intensities allowing the solitary wave to detach and propagate independently of the hull.

Shockwaves produced when a craft breaches the critical depth Froude number are shown by Jiang et al. (2002) and Grue (2017) using Boussinesq-type modelling, however, both of these pieces of research employed a transition in the water depth to create a rapid and significant increase in the depth Froude number. The authors’ previous work on fully confined waterways (Terziev et al. 2020) investigated the water depth transition problem using a URANS approach.

As shown above, a number of studies have examined problems similar to that explored herein. However, no previous research has combined fully non-linear Navier-Stokes-based methodologies with potential flow solutions to investigate the problem at hand. By consequence, the viscous influence on the ship’s resistance under the conditions examined subsequently remains unknown. This gap in the literature alongside the fact that few studies have examined the trans-critical resistance of an accelerating ship will be addressed within the study presented here.

The hydrodynamic forces that depend on acceleration are commonly referred to as added mass. However, it has been established in many previous studies that the magnitude of added mass, or the coefficient of added mass, primarily depends on the shape of the object and its velocity, rather than the acceleration. For example, Flagg and Newman (1971) and Newman's (2018) research investigated the influence of the width of a ship and the water depth on the added mass while Ghassemi and Yari (2011) studied the added mass coefficient of sphere, ellipsoid and marine propellers. Wakaba and Balachandar (2007) studied the added mass force at finite Reynolds and acceleration numbers. They found the added mass coefficient is acceleration independent. A similar conclusion is also mentioned in Javanmard et al.'s (2020) research, which focused on using a CFD method to calculate the translational added mass coefficients of underwater vehicles. However, it should be noted that in this paper, accurately reflecting the contribution of added mass to wave-making resistance is challenging. This is because the hydrodynamic forces in Equation (12) are not decomposed to extract distinct added mass forces. Instead, the unsteady hydrodynamic forces are determined by directly solving for the velocity potential φ, where the influence of added mass is inherently incorporated within the fully unsteady velocity potential φ. In fact, in other studies (Wehausen 1961; Shebalov 1970; Doctors and Sharma 1972; Pinkster 2004) related to acceleration, the added mass as a component in the Bernoulli equation has indeed been mentioned. However, it is worth noting that it has not been specifically isolated and studied separately in those cases.

3. Case study selection

The Wigley hull, with particulars shown in Table 1, is chosen to perform case studies. To the best of the authors’ knowledge, no experimental campaign has studied a trans-critically accelerating ship in shallow water. Selecting the Wigley hull maximizes the potential for other researchers to re-use the results presented here in their future studies because the hull is easily modelled mathematically and is frequently used for potential flow-based studies (Yuan et al. 2014; Andrun et al. 2018; Li et al. 2019; Bašić et al. 2020).

Table 1. Wigley hull principal dimensions.

Parameter	Symbol	Value	Units
Length	$L$	3	m
Draft	$T$	0.1875	m
Beam	$B$	0.3	m
Displacement	$\mathrm{\nabla }$	0.076	m^3

The depth Froude number gives the parameters varied to produce a test matrix. Namely, the ship speed, $V$, and the water depth, $d$. Following Day et al. (2009) and Li et al. (2019) multiples of the gravitational acceleration are used to vary the ship acceleration. The aforementioned authors used values between 0.08$g$ and 0.02$g$, but their studies focused on a rapidly accelerating hull to a target speed and maintaining that speed. By contrast, the objective of the current study is observing transient phenomena during acceleration. For this reason, the maximum acceleration chosen in this study is 0.02$g$, that is, the lowest acceleration intensity used by Day et al. (2009) and Li et al. (2019). Additionally, acceleration intensities equal to 0.01$g$, 0.005$g$, and 0.002$g$ are used. These values are combined with depth-to-draft ratios of $d/T=$1.2, 1.5 to gauge the effect of water depth on the results. The test matrix is given in Table 2. Each case is allowed to develop up to and including a depth Froude number of 2. In addition to the accelerating case studies, steady, constant velocity cases were modelled using CFD to estimate the deviation between the fixed speed resistance and accelerating resistance. The depth Froude numbers modelled at constant velocities are also given in Table 2.

Table 2. Adopted case studies.

Case number	Acceleration	Depth Froude number	$d/T$
1	0.02$g$	Variable	1.2
2	0.01 $g$
3	0.005 $g$
4	0.002 $g$
5	0.02 $g$		1.5
6	0.01 $g$
7	0.005 $g$
8	0.002 $g$
9	0	0.6	1.2
10		0.8
11		0.89
12		0.95
13		1
14		1.05
15		1.2
16		1.4
17		1.7
18		0.6	1.5
19		0.8
20		0.89
21		0.95
22		1
23		1.05
24		1.2
25		1.4
26		1.7

4. Methodology

The methodology is split into two approaches used in this study. Specifically, the potential flow approach and CFD.

4.1. Potential flow method

The potential flow solver, MHydro (Yuan 2019), based on 3D boundary element method (BEM) by using a Rankine type Green function is introduced to study unsteady free-surface effects. This solver has previously been used in the study of deep and shallow water ship hydrodynamics, for example, by Yuan et al. (2019). In order to obtain the solution, it is necessary to correlate all results at each time step, so the unsteady terms in free-surface boundary conditions are preserved. These enable the simulation of unsteady effects on the free surface. The specific methodology is introduced in detail in the following sections.

4.1.1. Description of the problem

To describe the problem, two right-handed coordinate systems shown in Figure 2 are defined. As can be seen from that figure, the global coordinate system O-XYZ is fixed on the earth while the local coordinate system $o_0-x_0{y}_0{z}_0$ is fixed to the ship hull and moves forward with the ship. In both frameworks, the positive direction of the $x$-axis points to the ship bow, the positive direction of the y-axis points to the port side of the ship and the positive $z$-axis points to upward. The undisturbed calm-free surface is located at $z$= 0.

Figure 2. Sketch of the problem. In (a), L is the length of the ship and U(t) is the ship velocity which is changing with time. T represents the draft of the ship while the water depth is denoted as d. (b) shows a top view of the problem, where it can be seen that the distance between the front and rear boundaries of the computational domain and the midship point are both 3L. On the other hand, both the port and starboard boundaries are L away from the centreline of the ship. Finally, B in this figure represents the ship breadth.

4.1.2. Boundary value problem

There are only two boundary conditions in this problem: the body-surface condition and the free-surface condition. The body-surface condition is that the surface is impenetrable and the normal velocity of all cells on the body is equal to the speed $U(t)$ of the ship: (1) ${\displaystyle \frac{\mathrm{\partial }\phi }{\mathrm{\partial }\mathit{n}}=U{n}_1}$(1) where $n=(n_1,n_2,n_3)$ is the unit normal vector on body-surface.

The free-surface condition can be divided into two parts: the dynamic free-surface condition and the kinematic free-surface condition. Because the method described herein is based on potential flow, the fluid is assumed to be ideal which means it is inviscid, irrotational, and incompressible. Then, through the velocity potential $\phi (x,y,z,t)$ and wave elevation $\zeta (x,y,t)$, the free-surface conditions can be described as: (2) ${\phi }_t-U{\phi }_x+g\zeta +{\displaystyle \frac{1}{2}({\phi }_x^2+{\phi }_y^2+{\phi }_z^2)+{\displaystyle \frac{P}{\rho }=0,onz=\zeta }}$(2) (3) ${\zeta }_t-U{\zeta }_x+{\phi }_x{\zeta }_x+{\phi }_y{\zeta }_y-{\phi }_z=0,onz=\zeta$(3) where $g$ is the gravitational acceleration, $\rho$ is the fluid density and $P$ is the pressure on free surface, and ${\phi }_t$ and ${\zeta }_t$ are the respective time derivatives. Similarly, ${\phi }_x$, ${\phi }_y$ and ${\phi }_z$ are directional derivatives.

To simplify the solution process, the nonlinear terms are ignored, and the unsteady terms are preserved as mentioned previously. The simplified linear unsteady free-surface condition can therefore be expressed as: (4) ${\phi }_t-U{\phi }_x+g\zeta +{\displaystyle \frac{1}{2}({\phi }_x^2+{\phi }_y^2+{\phi }_z^2)=0,onz=\zeta }$(4) (5) ${\zeta }_t-U{\zeta }_x-{\phi }_z=0,\mathrm{on}z=\zeta$(5)

Since the boundary conditions have been obtained, the next step is to solve the resulting equations. To that end, a three-level scheme is used to discretise the free-surface conditions: (6) $({\phi }_t{)}_{i,j}^{m+1}={\displaystyle \frac{1}{\mathrm{\Delta }t}\left({\displaystyle \frac{3}{2}{\phi }_{i,j}^{m+1}-2{\phi }_{i,j}^m+{\displaystyle \frac{1}{2}{\phi }_{i,j}^{m-1}}}\right)}$(6) (7) $({\zeta }_t{)}_{i,j}^{m+1}={\displaystyle \frac{1}{\mathrm{\Delta }t}\left({\displaystyle \frac{3}{2}{\zeta }_{i,j}^{m+1}-2{\zeta }_{i,j}^m+{\displaystyle \frac{1}{2}{\zeta }_{i,j}^{m-1}}}\right)}$(7) where $m$ represents the $m^{\mathrm{th}}$ time step, $i$ and $j$ indicate location of the cell on the free surface. Substituting Equations (6) and (7) into the kinetic free-surface condition, allows one to obtain the value of $(\zeta {)}_{i,j}^{m+1}$: (8) ${\displaystyle \frac{1}{\mathrm{\Delta }t}\left({\displaystyle \frac{3}{2}{\zeta }_{i,j}^{m+1}-2{\zeta }_{i,j}^m+{\displaystyle \frac{1}{2}{\zeta }_{i,j}^{m-1}}}\right)-U{\displaystyle \frac{1}{\mathrm{\Delta }x}\left({\displaystyle \frac{3}{2}{\zeta }_{i,j}^{m+1}-2{\zeta }_{i,j+1}^{m+1}+{\displaystyle \frac{1}{2}{\zeta }_{i,j+2}^{m+1}}}\right)-{\displaystyle \frac{\mathrm{\partial }{\phi }_{i,j}^{m+1}}{\mathrm{\partial }z}=0}}}$(8) (9) $\left({\displaystyle \frac{1}{\mathrm{\Delta }t}-{\displaystyle \frac{3}{2}{\displaystyle \frac{U}{\mathrm{\Delta }x}}}}\right){\zeta }_{i,j}^{m+1}+2{\displaystyle \frac{U}{\mathrm{\Delta }x}{\zeta }_{i,j+1}^{m+1}-{\displaystyle \frac{1}{2}{\displaystyle \frac{U}{\mathrm{\Delta }x}{\zeta }_{i,j+2}^{m+1}={\displaystyle \frac{2}{\mathrm{\Delta }t}{\zeta }_{i,j}^m-{\displaystyle \frac{1}{2}{\displaystyle \frac{1}{\mathrm{\Delta }t}{\zeta }_{i,j}^{m-1}+{\displaystyle \frac{\mathrm{\partial }{\phi }_{i,j}^{m+1}}{\mathrm{\partial }z}}}}}}}}$(9)

Then, all variables in the dynamic free-surface condition can be updated with the results obtained: (10) ${\displaystyle \frac{1}{\mathrm{\Delta }t}\left({\displaystyle \frac{3}{2}{\phi }_{i,j}^{m+1}-2{\phi }_{i,j}^m+{\displaystyle \frac{1}{2}{\phi }_{i,j}^{m-1}}}\right)-U{\displaystyle \frac{\mathrm{\partial }{\phi }_{i,j}^{m+1}}{\mathrm{\partial }x}+g{\zeta }_{i,j}^{m+1}=0}}$(10) (11) ${\displaystyle \frac{3}{2}{\displaystyle \frac{1}{\mathrm{\Delta }t}{\phi }_{i,j}^{m+1}-U{\displaystyle \frac{\mathrm{\partial }{\phi }_{i,j}^{m+1}}{\mathrm{\partial }x}={\displaystyle \frac{2}{\mathrm{\Delta }t}{\phi }_{i,j}^m-{\displaystyle \frac{1}{2}{\displaystyle \frac{1}{\mathrm{\Delta }t}{\phi }_{i,j}^{m-1}-g{\zeta }_{i,j}^{m+1}}}}}}}$(11)

In the above process, as the number of iterations $k$ increases, all variables will be continuously updated until both $|{\phi }_{i,j}^{m+1,k}-{\phi }_{i,j}^{m+1,k-1}|<\epsilon$ and $|{\zeta }_{i,j}^{m+1,k}-{\zeta }_{i,j}^{m+1,k-1}|<\epsilon$ are satisfied.

In order to solve the above-mentioned boundary value problem, we developed a programme, namely MHydro, based on the Rankine source panel method. Once the unknown potential is obtained, we can obtain the pressure and forces/moments acting on the hull using Bernoulli’s equation, which reads: (12) $P=-\rho \left[{\phi }_t-U{\phi }_x+{\displaystyle \frac{1}{2}({\phi }_x^2+{\phi }_y^2+{\phi }_z^2)}\right]$(12) (13) $F_i=\underset{S}{\int \int }P{n}_{i}ds,i=1,2,\dots ,6$(13) where $i$ represents the degree of freedom, and (14) $n_i=\{\begin{array}{cc}n,& i=1,2,3\\ x\times n,& i=4,5,6\end{array}$(14) Different from the method used by Pinkster (Pinkster 2004), this study does not decompose the hydrodynamic forces in Equation (12) to extract separate added mass forces. Instead, the velocity potential $\phi$ is directly solved to obtain the unsteady hydrodynamic forces.

4.2. Numerical set up

In addition to the potential flow methodology described previously, we employ the commercially available RANS solver Star-CCM+, version 16.06.010-r8 is to perform numerical simulations of the same problem and compare results. The chosen solver makes use of the finite volume method which relies of the integral form of the governing equations and subdivides the computational domain into a finite number of adjoining cells.

The turbulence model employed for the present study is the standard k-ω Wilcox model (Wilcox 2008), as implemented in the RANS solver. The k-ω is chosen based on results obtained by Terziev et al. (2019a), where the selected turbulence closure was shown to be approximately 16% more computationally efficient than other two-equation eddy-viscosity turbulence models. Second-order discretisation is used to ensure accurate modelling of the flow features following recommendations by Andrun et al. (2018).

Time is advanced in intervals of $\mathrm{\Delta }t=$0.01 following Song et al. (2019). The temporal term in the governing equations is discretised using a first-order scheme to alleviate the requirements on the Courant number $C<1$. The discretisation uncertainty stemming from the choice of time step is examined in subsequent sections.

The water surface and the ship-generated disturbance are modelled using the volume of fluid (VoF) approach. To improve the sharpness of the air–water interface, the high-resolution interface capturing (HRIC) scheme is employed.

4.2.1. Computational domain and boundary conditions

A rectangular box is used in all numerical simulations to represent the computational domain. The computational domain contains a symmetry plane coincident with the ship centreline to reduce cell numbers. Inlet, outlet boundaries are placed four ship lengths upstream and downstream from the aft perpendicular, respectively. The side boundary is also placed at a distance of four ship lengths from the ship centreline. These distances are considerably greater than the recommendations of the ITTC (International Towing Tank Conference) (ITTC 2014) to ensure transient flow features can be captured without any influence exerted by the domain dimensions. The domain top is placed approximately one ship length from the undisturbed water level. Finally, the domain bottom is assigned a no slip wall and placed in accordance with the test matrix given in Table 2. The ship hull is also set as a no slip wall.

The inlet is assigned a velocity inlet boundary condition which maintains the water level, while the outlet boundary enforces the hydrostatic pressure through a pressure outlet boundary condition. The side boundary and symmetry boundary are symmetry planes. Wave damping is disabled from all boundaries. The resulting computational domain and boundary conditions are shown in Figure 3.

Figure 3. Computational domain dimensions and boundary conditions ($L$= 3 m). Top: side view; bottom: front view.

4.2.2. Ship motion modelling in the URANS framework

The approach taken to modelling forward ship motion in the URANS framework is to use moving frames of reference. Alternative approaches, such as invoking frame invariance to assign the fluid’s speed equal but opposite to that of the ship, cause problems in the modelling of acceleration. Within the moving frames of reference approach, the coordinate system is assigned the ship velocity. Adopting such an approach allows the surrounding flow and domain bottom boundary to remain static with respect to the earth-fixed coordinate system creating a relative velocity between the hull and the environment. This permits the velocity to change at each time step in accordance with the case studies given in Table 2 while maintaining no relative motion between the undisturbed fluid and seabed. In the frame of reference of the vessel, the domain bottom and fluid translate in the negative x direction at the time-varying velocity specified in Table 2.

4.2.3. Mesh generation

Mesh generation is performed within the automatic facilities of Star-CCM+. Hexahedral cells are used throughout all numerical simulations.

As shown in Figure 1, the Kelvin wake angle is expected to vary from the deep-water value of approximately ${19.47}^{\circ }$ to ${90}^{\circ }$ during the numerical simulation. Therefore, typical mesh refinements targeting the Kelvin wake are not sufficient to capture disturbance generated by the hull. The approach to meshing taken in this study is to enforce a cell size at the hull and assign a doubling in cell dimension every 20 cells in all directions in addition to the usual Kelvin wake refinement. In taking this approach, the flow features near the hull are modelled accurately, while disturbances extending far beyond the hull are gradually coarsened. Additionally, the slow growth in cell size ensures that flow features that do not fall within the Kelvin wedge are resolved near the hull.

All numerical simulations use a high $y^{+}$ approach, which determines the cell size near the hull. Since the ship speed will be 0 m/s at the start of the simulation, it is not possible to guarantee a $y^{+}$> 30 at all times. However, it is possible to confine the influence of the buffer layer (5 < y⁺< 30) to low depth Froude numbers which are of little interest to the present study. The approach described in Terziev et al. (2022) is used to confine the buffer layer below $F_h=$0.15. The aforementioned approach specifies the near-wall grid layers as a geometric series with a common ratio, $S$, defining the stretch factor between adjacent layers. The number of layers, $n$, are given in Equation (15): (15) $n=\log \left(1-{\displaystyle \frac{\delta (1-S)}{2\mathrm{\Delta }y}}\right)/\log (S)$(15) where $\delta$ is the total thickness over which layers are distributed equal to 0.01 m in this study, and $\mathrm{\Delta }y$ is the location of the first cell centre. To estimate $\mathrm{\Delta }y$, the ITTC correlation line is used: (16) $C_f=0.075/({\log }_{10}Re-2{)}^2$(16) where $Re=VL\rho /\mu$ is the Reynolds number with $\mu =8.8871\times {10}^{-4}$ Pa-s being the dynamic viscosity, $C_f$ gives the local shear stress ${\tau }_w=C_f\rho V^2/2$, where $\rho =997.561$ kg/m³ is the freshwater density. Once these parameters are known, the first cell centre is $\mathrm{\Delta }y=y^{+}\upsilon /u_t$ with $u_t=\sqrt{{\tau }_w/\rho }$ and $\upsilon =\mu /\rho$ (Peric 2019). Using these relations, the $y^{+}$ value is set to remain below 300 at the highest speed corresponding to $F_h=$2. The resulting computational mesh consists of approximately 5.3 million cells and is depicted in Figure 4.

5. Results and discussion

Three sets of resistance results are extracted from the CFD solution, consisting of the total resistance ($R_T$), frictional resistance ($R_F$), and pressure resistance ($R_P$). The latter components represent the normal and tangential components of the total, respectively. As such, the frictional resistance includes the effects of surface curvature, while the pressure resistance represents a combination of the viscous pressure resistance and wave resistance (Molland et al. 2017). The viscous pressure resistance typically attains a small portion of the total at all speeds (Terziev et al. 2021a). On the other hand, the potential flow theory used to compare results provides wave resistance. Since viscous pressure resistance is small, the results from the potential flow method and URANS method are compared directly. Splitting the viscous pressure resistance from the pressure resistance would require so-called double-body numerical simulations, where the water surface is replaced by a rigid symmetry plane. However, since the free surface deforms considerably near the critical depth Froude number, double-body simulations will not provide an adequate estimate of the viscous pressure resistance. In light of these facts, the numerically obtained pressure resistance is compared directly with the wave resistance calculated using the potential flow approach.

Figure 4. Computational mesh used to perform numerical simulations. Top: View of the entire domain; middle: close-up around the hull; bottom: surface mesh on the hull. Total cell number $\approx$ 5.3 million.

5.1. Pressure and wave resistance

This section compares the wave and pressure resistance values obtained as the ship accelerates for the two depth-to-draft ratios modelled. Figure 5 shows the comparison between the CFD method and MHydro for all cases listed in the test matrix (Table 2) in addition to the constant velocity cases modelled using CFD (Table 2). All forces are made dimensionless using the ship mass force $F_m$, following Day et al. (2009). A pronounced peak can be observed in all cases in the vicinity of $F_h=1$. In the steady problem (i.e. constant speed with no acceleration), such a peak is produced by the mutual interference of the fore and aft wave systems generated by the ship. The resistance curves in Figure 5 show that this kind of interference reaches its maximum between $F_h=1.05$ and $F_h=1.2$ depending on the value of the dimensionless acceleration ($a/g$).

Figure 5. Changes of resistance during ship acceleration. The tiles on the left show d/T = 1.2 while the tiles on the right show d/T = 1.5. F_x/F_m (with F_m equals to mg) is the dimensionless wave resistance obtained by MHydro and pressure resistance obtained through the URANS method described previously. It also should be noted that the results of CFD here are pressure resistance while MHydro results are only wave making resistance.

It should be noted that because a limited number of points are selected for testing the constant speed case, the corresponding peak values shown in Figure 5 will not necessarily indicate the exact peak, since that curve can only provide an approximate range of where the peak value will appear. That increase in resistance is clearly influenced by the water depth, as evidenced by the higher F_x values when d/T = 1.2 relative to d/T = 1.5. This phenomenon changes when acceleration is added. In other words, the unsteady wave generated by the existence of acceleration has a significant impact on the resistance.

The depth Froude number at which CFD-derived resistance peaks occur for d/T$=1.2$ decreases by 10.25% when the acceleration halves from $a=$0.02$g$ to $a=$0.01$g$. A further halving of the acceleration to 0.005$g$ reduces the $F_h$ at which the resistance peak occurs by a further 5.13% relative to the 0.01$g$ case. Reducing the acceleration further has a marginal effect of 0.94% in the depth Froude number. A similar monotonic decrease in the position of the peak is found when d/T$=1.5$, but the range reduces to 2.21% between $a=$0.005$g$ and $a=$0.002$g$. In other words, at larger water depths, the position of the peak in resistance shows greater variation and sensitivity to acceleration than in very shallow waters.

The difference between the magnitude of the peaks predicted by CFD and MHydro is between 5.54% and 12.44% for $d/T=1.2$, and 6.23% and 11.75% for $d/T=1.5$. When the water depth corresponds to $d/T=1.2$ the difference between the two predictions grows by approximately 3% at each step between $a=$0.02$g$ to $a=$0.005$g$. At $a=$0.002$g$, the disagreement between the two solvers peaks at 12.44%. In the case where d/T$=1.5$, that effect is only partially replicated. That is to say, the disagreement grows with decreasing acceleration, but this is only up to $a=$0.005$g$, where the discrepancy peaks at 11.75%. A further reduction in the acceleration decreases the difference by approximately 3%.

In terms of the disagreement in the position of the two peaks, the solvers predict remarkably close results for high accelerations. Interestingly, the difference between the predictions grows monotonically from 0.75% and 0.68% to 10.13% and 5.54% for d/T$=1.2$ and d/T$=1.5$, respectively. Further investigations are necessary to determine why the two sets of predictions performed better for the lowest acceleration when the water depth is slightly larger. Such an effect could be explained through the presence of non-linear effects which are known to grow in importance with decreasing water depth, but such conclusions require more detailed investigations, preferably using experimental data to eliminate any bias in numerical or potential solutions.

When comparing the results obtained from MHydro and CFD, it can be concluded that no matter whether d/T = 1.2 or d/T = 1.5, the peak of results from MHydro will appear before the peak predicted by CFD, and the resistance value obtained from MHydro is always larger. The phenomenon is especially obvious when the acceleration is small. It is postulated that this difference arises due to viscous and nonlinear terms. Nonlinear terms are ignored in MHydro, which means that MHydro provides a linear and unsteady result while the CFD results are nonlinear and unsteady. As the acceleration increases, the unsteady effect increases in relative importance eclipsing nonlinear effects. Therefore, the unsteady effect occupies a dominant position when the acceleration is sufficiently large, which is why the results obtained from MHydro and CFD show less agreement when the acceleration is small. Conversely, the difference between the results obtained from the two prediction methods is smallest when the acceleration is large. The discussion above is also supported by observations made elsewhere regarding the linear potential flow solutions near the critical depth Froude number (Tuck and Taylor 1970; Lea and Feldman 1972; Tuck 1978).

5.1.1. Acceleration effects

As stated previously, there are two main influencing factors in the problem examined herein: the magnitude of the acceleration and the water depth. These are discussed in turn. It can be seen from Figure 6 that after the introduction of acceleration, both the peak value of the resistance and the position where this peak appears change. Having established good agreement between CFD and MHydro, additional case studies are explored in Figure 6 with the potential flow solver only due to its faster turnaround time and lesser computational requirements.

Figure 6. (a) Changes of fully unsteady resistance with varying acceleration; (b) Resistance peak value for different accelerations; (c) Peak location at different acceleration; Cases depicted correspond to d/T = 1.2.

Figure 6 shows that the greater the acceleration, the smaller the resistance peak and the later it appears along the examined depth Froude number range. It can therefore be concluded that the unsteady effect induced by the acceleration has a significant influence on the results. The shifting in position of the peak in resistance can be explained by the shock wave at the bow caused by the acceleration. In the initial stages when the ship speed is low a shock wave will be emitted forward causing additional wave/pressure resistance. As the speed of the ship increases to the same speed as the wave, the wave making resistance reaches a maximum. This phenomenon appears as the peak in Figure 5.

In other words, as the acceleration increases, the velocity of the resulting shock wave also increases, so the ship speed must be higher before the shockwave can be overtaken. Thus, the position of the peak is delayed at higher accelerations. Following this, as the ship continues to accelerate, the shock wave will be overtaken by the ship, so the resistance drops sharply for a narrow band of $F_h$ values, followed by a linear increase.

Another observation on the results presented herein is that regardless of acceleration, the time-history of the resistance passes through a single point, located close to F_h= 1. In all cases, the growth in resistance beyond that point is approximately linear until the peak of resistance is reached. These results are depicted in Figure 7(a and b) for d/T = 1.2 and d/T = 1.5, respectively.

Figure 7. Resistance time-history for all acceleration intensities while solid lines show the results obtained using CFD and dashed lines show the results obtained using MHydro. (a) Resistance time-history at d/T = 1.2 for CFD. (b) Resistance time-history at d/T = 1.5 for CFD. (c) Resistance time-history at d/T = 1.2 for MHydro. (d) Resistance time-history at d/T = 1.5 for MHydro.

When the water depth corresponds to d/T = 1.5, the CFD results show that resistance time-histories coincide at approximately $F_h=$0.906, while the potential flow results predict this point at approximately $F_h$ = 0.935. Similarly, when d/T = 1.2, CFD predicts a shift in this point to approximately $F_h$ = 0.92, while MHydro places this at approximately $F_h$ = 0.947. Prior to and following the point where resistance time-histories intersect, resistance follows unique paths. To the best of the authors’ knowledge, such phenomena have not been shown in the existing literature.

In order to study the cause of this phenomenon, wave fields under different accelerations near the intersection point are shown in Figure 8. When F_h or the acceleration is different, wave elevations at the bow and stern also change accordingly. The wave at the bow is mainly determined by the velocity of the ship, while the condition is different for waves at the stern. As shown in Figure 9(a), the wave elevation at the ship stern is inversely proportional to the acceleration. A greater acceleration will bring more pronounced unsteady effects, and when superimposed on the original wave field at the ship stern, results in smaller wave peaks and troughs. It should also be mentioned that the waves at the ship stern show peaks before the critical speed and troughs after it. This means that when $F_h<1$, the smaller the acceleration, the closer the wave elevation at the bow and stern, and the less the wave resistance. However, after this point, as the acceleration decreases, the wave trough at the stern will increase significantly. According to Bernoulli's equation, the pressure field at the stern of the ship will also decrease, resulting in greater wave resistance. This is why we observe the intersection of the wave resistance in Figure 7, and also why the peak value of the resistance decreases with an increase in acceleration.

Figure 8. Wave elevations on the port side of ship with different accelerations (y/L = 0.0167). These results were computed using MHydro for d/T = 1.2.

Figure 9. Wave elevations around the ship stern and bow with different accelerations. (a) Stern wave (WS); (b) Bow wave (WB); (c) Wave difference between ship bow and stern. These results were computed using MHydro for d/T = 1.2.

5.1.2. Finite depth effects

As evidenced by the agreement between the constant speed and accelerating results shown in Figure 5(d), the effect of acceleration at $a=0.002g$ and $d/T=1.2$ is weak. However, there is a significant difference between the corresponding values at d/T$=1.5$. The aforementioned observation, in conjunction with the fact that the peak shifts along the $F_h$ range with varying acceleration indicates that the depth Froude number is not the best parameter to represent the time-varying velocity. Figure 10 shows a set of example cases with a/g$=0.02$ for d/T$=1.2-4$ computed with MHydro. Under this framework, larger depth-to-draft ratios push the peak towards higher $F_h$ values and reduce its magnitude. That effect is most perceptible between d/T$=2$ and d/T$=4$. The resistance penalty induced as a shallow water effect at higher $F_h$ decreases rapidly with increasing water depth. If the water depth were increased beyond $d/T=4$, the characteristic peak of resistance observed in shallow waters would vanish.

Figure 10. Wave resistance and associated peaks at different water depths, computed using MHydro; (a) shows the time-history of the resistance, (b) shows the magnitude of the peak for varying $d/T$ ratios, (c) show the position of the peak as a function of $d/T$ while (d) shows the resistance value at $F_h=1$.

A second observation from Figure 10 is that shallow water may be advantageous in terms of reduced resistance when accelerating a vessel. Specifically, the values depicted for $F_h=1$ in Figure 10(d) demonstrate that the higher the water depth, the higher the resistance. For the given acceleration, the increase in resistance between very shallow water ($d/T=1.2$) and moderately shallow water ($d/T=4$) is more than 100%. That is a consequence of the delay in the appearance of the resistance peak to higher depth Froude numbers as shown in Figure 10(c).

5.2. Unsteady free-surface effects

The results by MHydro shown in the above sections are obtained by solving the fully unsteady Boundary Value Problem (BVP) defined in Equations (1), (4) and (5). It would be interesting to investigate the effect of unsteadiness on free-surface boundary condition. A quasi-steady solution is therefore proposed, in which the time-dependent terms in Equations (4) and (5) are ignored at each time step, while the body-surface boundary condition in Equation (1) remains the same. For both unsteady and quasi-steady methods, the unsteady Bernoulli's equation in Equation (12) is used to calculate the hydrodynamic force, more specifically, the wave-making resistance in the present study (Figure 11).

Figure 11. Comparison of resistance values for different accelerations and water depths; (a) and (b) show $d/T=1.2$ and $d/T=1.5$, respectively; (c) and (d) show the magnitude of the resistance peak as a function of the acceleration and its position along the $F_h$ range, respectively.

Figure 12(a) shows a direct comparison of the results obtained from quasi-steady and unsteady free-surface solutions mentioned above. The discrepancies indicate the importance of the unsteady free-surface condition on the present acceleration problem. The free-surface effects are predominant at near critical speed region. Due to the unsteady nature of the acceleration problem, the unsteady terms on free-surface condition must be considered.

Figure 12. Wave-making resistance obtained by MHydro and double-body flow method. (a) MHydro results with and without considering the time-dependent terms on free-surface boundary condition; (b) the unsteady component in MHydro results; (c) double-body flow results without considering acceleration; (d) double-body flow results with considering acceleration. It is worth mentioning that the result in (c) represents a fully steady state, while the results in (d) are quasi-steady.

The impact of unsteady effects on wave-making resistance is more evident in Figure 12(b). It is observed that the influence of unsteadiness is minimal when the velocity is below the critical speed (F_h <1). However, near the critical speed, there is a significant increase in unsteady effects. Once the velocity exceeds the critical speed (F_h >1), the intensified unsteady effects, accompanied by increasing acceleration, lead to a balanced wave field around the ship hull, resulting in a reduction in wave-making resistance. This conclusion aligns with the earlier analysis presented in Figure 9.

Figure 12(c and d) illustrate the results of the double-body flow method with and without considering acceleration, respectively. It is evident that regardless of considering or not considering acceleration, the results obtained from the double-body flow method are significantly smaller than those obtained from MHydro. This indicates in the case of a ship continuously accelerating in extremely shallow water, it becomes evident that the contribution of the wave field, especially the unsteady waves, dominates the wave-making resistance. This is also why the double-body flow method fails to capture the variations in wave-making resistance closely associated with the intricate wave effects when the free surface is absent.

Figure 12(c and d) present the results of the double-body flow method with and without considering acceleration, respectively. It is clearly evident that regardless of whether acceleration is considered or not, the wave-making resistance values obtained from the double-body flow method are significantly smaller compared to those obtained from MHydro. This observation highlights the pronounced influence of the wave field, particularly the unsteady waves, on wave-making resistance when a ship undergoes continuous acceleration in extremely shallow water. Furthermore, it explains why the double-body flow method fails to accurately capture the intricate variations in wave-making resistance associated with the presence of the free surface.

5.3. Frictional resistance

Previous sections examined the pressure and wave resistance results obtained through CFD and MHydro, respectively. However, an important component of the resistance is due to friction. That component can only be predicted through the CFD results because friction stems from viscous effects. This sub-section examines possible acceleration effects in the frictional component. The results obtained in the accelerating cases and steady cases are depicted in Figure 13.

Figure 13. Frictional resistance predicted by the CFD solver for all cases, and comparison with the shallow water friction line proposed by Zeng et al. (2019).

It should be noted that at high accelerations, the frictional resistance grows smoothly with no perceptible variation near the critical depth Froude number. However, when the acceleration is low ($a\le$0.005$g$) or when the speed is constant, a pronounced hump spanning approximately $0.8<F_h<1.2$ is visible. This may be explained by the change in the wave field and therefore wetted area of the hull as one switches from subcritical to supercritical flow regimes.

When the depth Froude number is less than unity, increasing the acceleration causes a monotonic increase in the frictional resistance. For example, a cubically interpolated value of the frictional resistance for $F_h=0.6$ at $d/T=1.5$ increases by 0.64% as a result of an increase in acceleration from $a=$0.002$g$ to $a=$0.005$g$. A further increase in acceleration magnifies friction by 3% relative to $a=$0.005$g$. That effect is replicated for both water depths examined in the present study. However, the effect is more pronounced when the water depth is smaller, ranging from 0.83% to 4.01% at $F_h=$0.6.

For a critical depth Froude number, the change in friction due to acceleration is approximately 1% or less depending on the water depth. Further increases in speed reverse the relationship discussed above at low depth Froude numbers. Namely, the frictional resistance decreases monotonically, but always by about 1% or less. To the best of the authors’ knowledge, such effects have not been demonstrated elsewhere.

To provide a reference for comparison, Figure 13 includes the results obtained using the friction line proposed by Zeng et al. (2019). The difference in predictions grows with increasing depth Froude numbers. Such an observation is expected because the friction line of Zeng et al. (2019) was devised through a regression using double-body CFD results for relatively low $F_h$. For that reason, our results agree well in the subcritical range. Free-surface effects on the frictional resistance have been documented previously with similar magnitude (Terziev et al. 2021b, 2019b), meaning that that the gap between our results and the Zeng et al. (2019) friction line is expected. It should also be noted that a better agreement is found for the higher water depth examined.

In summary, the acceleration effects on frictional resistance can be ignored. However, the free-surface effects may influence the magnitude of the frictional resistance by way of changes produced by ship-generated waves. Application of the friction lines at very high speeds should be done with care since such cases are typically outside of a friction line’s range of applicability, particularly when such relationships are established through CFD-based regressions.

5.4. Verification study

All solutions derived from numerical solvers rely on the mapping of the governing continuous partial differential equations onto discrete intervals in space, and when unsteady solvers are used as is the case herein, time (Roy and Oberkampf 2002). Doing so introduces a discretisation error which must be quantified (Freitas 2020). In the present paper, an extension of Richardson Extrapolation (RE), known as the Grid Convergence Index (GCI) is used to estimate the discretisation uncertainty. The uncertainty is a symmetrical band around the solution showing the 95% confidence interval achieved with the mesh and time step choices. The GCI methodology is described by Celik et al. (2008), is openly available in the literature (ITTC 2008) and is described qualitatively only in this study.

The procedure begins by coarsening the target grid and time step by a constant factor, known as the refinement ratio, $r=\sqrt{2}$, following the recommendations by the ITTC (2017) and American Society of Mechanical Engineers (ASME 2009). Coarsening the solution mesh and time step twice sequentially gives rise to the fine, medium, and coarse solutions. The approach taken to coarsening follows Burmester et al. (2020) who recommended the Courant number be kept the same across the so-called triplet of solutions. Therefore, when the grid dimension is multiplied by $r$ and $r^2$, the time step is also multiplied by the same factor. The resulting grid numbers were approximately 5.3, 2.8, and 1.6 million for the fine, medium and coarse grids, respectively. Case 8 from Table 2 ($d/T=$1.5 accelerated by 0.002$g$) is chosen for the uncertainty assessment, whose results are given in Figure 14. It is assumed that other case studies will show a similar level of discretisation uncertainty because identical numerical settings are used throughout. In addition, as shown previously, the low acceleration cases exhibit the highest changes in the time-history of the resistance. It is therefore expected that these cases will be most challenging for the numerical solver, and consequently represent a good choice for the discretisation uncertainty estimation exercise.

Figure 14. Fine (solid line), medium (broken line), and coarse (dotted line) numerical solutions for d/T = 1.5 and an acceleration intensity of 0.002$g$. Black, blue, and red represent the total, pressure, and fricitonal resistance, respectively. Each parameter’s discretisation uncertainty is given as a shaded area around the fine solution estimated using the Grid Convergence Index method.

The uncertainty, shown in Figure 14, is dominated by the pressure resistance, whose influence is translated into the uncertainty in the total resistance. The highest level of uncertainty is equal to approximately 3.5 N, 0.0048 in coefficient form, or 8.49% for the total resistance at the peak (located at $F_h=$1.256). These levels of maximum uncertainty are tolerable considering that the total uncertainty band throughout the grid triplet spans the distance between approximately 40 and 45 N, or 0.0545–0.0613 in coefficient form. Locations where the uncertainty is elevated correspond to inflection points in the resistance curve. Specifically, near the peak and subsequent dip in the pressure and total resistance. The frictional resistance exhibits little uncertainty throughout the depth Froude number range; a result replicated in other studies (Korkmaz et al. 2021; Terziev et al. 2021a).

The 95% confidence interval of less 8.49% at the peak points to the possibility that further mesh refinements may bring the CFD results closer to those of MHydro, which was 12.44% larger for the examined case. However, it should be noted that the GCI procedure employed to predict the uncertainty is incapable of predicting the bias in the numerical model, meaning that further refinement is as likely to reduce the gap between the two methods as it is to increase it. The present study did not attempt further refinement due to resource restrictions: a mesh, refined systematically by a further factor of $r=\sqrt{2}$ would likely consist of between 10 and 11 million cells. The further reduction in the time step places too great a strain on the resources available for this study. In addition, the generated mesh produced results in excellent agreement across the accelerations and depth-to-draft ratios examined.

6. Conclusions and recommendations for future work

A craft accelerating past the speed of the fundamental wave in a given medium experiences a fluctuation in resistance. To the best of the authors’ knowledge, such effects had not previously been documented for the case of a steadily accelerating ship in shallow water, where the wave speed is a constant value which is a function of the gravitational acceleration and water depth only. The present study aimed to fill this gap in the literature by combining URANS and potential flow solutions.

Four acceleration intensities and two water depths were modelled using the Wigley parabolic hull. In all cases, a pronounced peak in the resistance was observed near the critical depth Froude number. The results show a maximum disagreement in the prediction of the resistance peak’s location along the examined depth Froude number between 5.54% and 12.44% for d/T = 1.2, and 6.23% and 11.75% for d/T = 1.5, respectively.

Results showed excellent agreement in the low and high-speed range, past the critical transition, indicating that the resistance in that range is linear and dominated by unsteady rather than nonlinear terms. Some deviation was found between the two solvers near the critical speed, which was interpreted as an indication that nonlinear effects influence the solution in that range. A point near but less than $F_h=1$ was identified for both water depths modelled where all acceleration intensities cross, creating an intersection point. The existence of this point was explained through the proximity of the trans-critical boundary and the fact that no steady flow is possible when the depth Froude number is unity.

Although we have confidence in the results presented herein in view of the agreement between the two solvers, experimental data of trans-critically accelerating ship resistance would be invaluable in demonstrating the exact level of accuracy. In addition, models of trans-critical acceleration in fully confined water are of particular interest since the blockage ratio would also be an important parameter. Such cases are also interesting because the trans-critical range can occupy a significant portion of the depth Froude number depending on the blockage (Lataire et al. 2012).

Acknowledgements

Results were obtained using the ARCHIE-WeSt High-Performance Computer (www.archie-west.ac.uk) house at the University of Strathclyde.

Disclosure statement

No potential conflict of interest was reported by the author(s).