3.1.4. Time-History of the Solution

An example time-history of the resistance of the ship is given in Figure 6. The figure is characterised by two distinct regions. Firstly, the ship velocity is linearly increased up to the target velocity. This is done by linking the time-step and velocity. The specific approach adopted requires the ship speed to reach its steady value at the end of the first 1000 time-steps. In other words, the ship's velocity is increased by *Utarget*/1000 each time-step. Thus, after an initial oscillatory behaviour, the resistance time-history exhibits oscillatory convergence towards its steady-state value. The oscillations are liked with the reflections of waves from the side walls, which also impact the observed resistance [39]. All final values reported in this work are obtained by averaging over one period of oscillation of the resistance curve. The specific point where averaging is performed is chosen to be sufficiently far from the acceleration phase to eliminate its effect.

**Figure 6.** Example time-history of the solution (depicted: *Fh* = 0.57).

#### 3.1.5. Verification

This section contains the numerical verification of the case study in the rectangular canal, *Fh* = 0.57. The numerical uncertainties are estimated via the grid convergence index (GCI) method. This is the standard way to report numerical uncertainties in ship CFD [40]. It is assumed that the remaining cases are characterised by similar levels of numerical uncertainty.

The GCI method requires three systematically coarsened solutions for the same case. In this study, the recommendations of the ITTC [40] are followed. Specifically, the refinement ratio, *<sup>r</sup>* <sup>=</sup> <sup>√</sup> 2 is adopted. The refinement ratio is used to coarsen and lessen the grid size and time-step, respectively. The GCI method assumes that all three solutions are close to the asymptotic range and are sufficiently different, which may be difficult to achieve in practice. The proximity to the asymptotic range is typically characterised by the convergence ratio, *p*, which is shown in Equation (1).

$$p = \ln(\varepsilon\_{32}/\varepsilon\_{21}) / \ln r \tag{1}$$

where ε<sup>32</sup> = *f3* − *f2* and ε<sup>21</sup> = *f2* − *f1.* Here, *fn* represents the *n*-th solution, generated by a systematic coarsening/lessening of the input parameter (grid or time-step). If *p* = 2, then the grid or time-step can be deemed asymptotic [41]. The convergence properties, however, can be determined in a different manner, which also carries information on what type of convergence/divergence is achieved with refinement. This is known as the convergence ratio, *R* = ε21/ε<sup>32</sup> [42]. Based on the value of *R*, the following may be interpreted:

1. Monotonic convergence is observed if 0 < *R* < 1;


The GCI method is only applicable in case (1). Next, an error estimate (*e*21) is defined as [43]:

$$c\_{21} = (f\_1 - f\_2) / f\_1 \tag{2}$$

Once the error is known, the numerical uncertainty can be calculated as shown in Equation (3) [44]:

$$\text{GCI} = 1.25e\_{21} / \left(r\_k^p - 1\right) \tag{3}$$

where *k* represents the *k*-th input variable (grid or time-step). The factor 1.25 in the numerator of the expression, defining the numerical uncertainty represents a factor of safety. This has been devised to ensure that the true solution lies within the bracket provided by the GCI with 95% confidence. The results from the convergence study can be seen in Table 4.


**Table 4.** Grid and time independence (rectangular canal, *Fh* = 0.57). Experimental result: 4.5047 N.

The successive grid coarsening resulted in 10,955,825 and 4,155,326 cells for the medium and coarse solutions, respectively. In terms of spatial dependence, the solution exhibited rapid ranges with reduction in cell numbers. This 'superconvergence' can be deduced by examining the order of accuracy, *p*. While in the case of grid coarsening it is approximately 9, when the time-step is lessened, the solution changes according to *ptime* = 0.463. According to Eca and Hoekstra [45], orders of accuracy between in the range 0.5 and 2 can still be treated as asymptotic. Eca and Hoekstra [45] devised a procedure based on a least-squares fit to estimate the numerical uncertainty. Their method is not employed in this study, because it requires a minimum of four solutions. Further coarsening of the computational mesh resulted in divergent behaviour in the simulation. The consequence of the time-step exhibiting the above order of accuracy means that the GCI method predicts large numerical uncertainties, even though the overall change between the coarse and fine solution is less than 3% of the fine solution's value.

It should be noted that the time-step is kept at the smallest value while coarsening the grid. Conversely, the finest grid was maintained throughout the temporal convergence analysis. To ensure that the ratio of overset cell to background cell dimension is kept constant, both domains were coarsened simultaneously. In this study, it is assumed that all examined cases will exhibit similar levels of spatial and temporal dependence. For this reason, the above procedure is not repeated.

#### *3.2. Spectral Representation of the Wave Field*

In this section, the method devised by Caplier et al. [20] and Gomit et al. [22] is briefly examined. This method has previously been used to determine a ship's speed via satellite imagery [46,47]. The essence of the approach is to process an available water surface in Fourier space. To achieve this, a 2D Fourier transform is used to represent the disturbance in the spectral space (*kx*, *ky*). These arise by defining the angular wavenumber, *k*, as a vector containing x and y components, *k* = *k*2 *<sup>x</sup>* + *k*<sup>2</sup> *<sup>y</sup>*. The spatial equivalent of these components are used to calculate the extents of the spectrum. The x-direction length of the entire water surface (*Lx*) becomes *kx*,*max*; similarly, using the extent in the y-direction (*Ly*) one obtains *ky*,*max*, as shown in Equation (4). Likewise, the resolution of the water surface in real space dictates the steps in Fourier space Δ*kx* and Δ*ky* in the x and y directions, respectively (Equation (5)).

$$k\_{\rm x,max} = 2\pi / L\_{\rm x} \; k\_{y,max} = 2\pi / L\_{y} \tag{4}$$

$$
\Delta k\_{\rm x} = \frac{1}{2} \times \frac{2\pi}{\Delta X}, \Delta k\_{\rm y} = \frac{1}{2} \times \frac{2\pi}{\Delta Y} \tag{5}
$$

where Δ*X* and Δ*Y* are the resolutions in the x and y directions, respectively.

If surface tension is ignored, the dispersion relation in shallow water may be expressed via the angular frequency (ω) of the waves as shown in Equation (6). This is justified, because surface tension becomes important only in waves characterised by wavelengths smaller than 7 cm. Alternatively, the travelling disturbance should propagate with a speed higher than 0.23 m/s [48] for surface tension to be negligible.

$$
\omega(k) = \pm \sqrt{\text{gktanth} h} \tag{6}
$$

A moving ship will cause the waves to be Doppler-shifted, and setting ω = ω − *Ukx* and *k* = *k*, the resulting dispersion relation becomes [20]:

$$
\omega'(k) = \pm \sqrt{\text{gktambkh}} - llk\_{\text{x}} \tag{7}
$$

To obtain the locus of the dispersion relation, Equation (7) is solved for ω (*k*) = 0 [49], which yields:

$$\ln l^2 k\_x^2 - g\sqrt{k\_x^2 + k\_y^2} \tanh\left(h\sqrt{k\_x^2 + k\_y^2}\right) = 0\tag{8}$$

Equation (8) is symmetrical with respect to both axes [50]. This is demonstrated in Figure 7, which includes computed loci (this is used to represent the solutions of Equation (8)) for *Fh* = 0.57, 0.77, according to the adopted cases. Alongside these, the critical depth Froude number is depicted to demonstrate the effect of ship speed on the dispersion relation in shallow water. It should be noted that the dispersion relation is speed-independent in deep water [20]. The arms of the loci always begin at *ky*,*x*/ *g*/*U*<sup>2</sup> = 1 in deep water. This is also the cut-off wavenumber. In shallow waters on the other hand, the cut-off wavenumber varies with speed. This can be seen by consulting Figure 7, specifically, where the curves cross the abscissa. Here, the case for *Fh* = 0.47 is not shown, as it is practically impossible to distinguish it from the *Fh* = 0.57 case. Deep and shallow water cases are essentially identical when *kh* 1.

**Figure 7.** Solutions to Equation (8). Figure depicts the examined loci alongside the critical depth Froude number (*Fh* = 1) to demonstrate the effect of speed.

The cut-off wavenumber separates the near-field disturbance from the far-field waves generated by the ship. Thus, useful analysis with applications to loads on coastal structures may be performed by removing the near-field disturbance, which does not propagate away from the ship. To determine the cut-off wavenumber (*k<sup>c</sup> <sup>x</sup>*) in shallow water, Caplier et al. [20] solved Equation (9).

$$
\hbar L^2 k\_x^\varepsilon - \mathcal{g} k\_x^\varepsilon \tanh(\hbar k\_x^\varepsilon) = 0 \tag{9}
$$

Finally, the Kelvin half-angle may be determined by computing tan θ = *dky*/*dkx* −<sup>1</sup> at the inflection point. According to Nakos and Sclavounos [51], numerical errors will manifest near cut-off wavenumbers. This can be deduced by examining Figure 7. Even a small deviation in the intersection between the locus and the abscissa will lead to large errors. The specific example given demonstrates the relatively short distance between the intersection points of *Fh* =0.57 and 0.77. On the other hand, as the speed is increased further, the locus approaches the origin. At the critical depth Froude number, the locus will transition into crossing the origin and progressing into a quadrant characterised by opposite signs of the abscissa and ordinate.

The manner in which a numerical free surface generates this pattern is in terms of maxima of the spectrum in Fourier space. This can be extracted and compared to the theoretical prediction provided by Equation (8). Caplier et al. [20] demonstrated that the relationship holds well despite its neglect of nonlinear and three-dimensional terms. Thus, if one can prove that a numerically generated free surface (once processed to the spectral domain) provides maxima near the locus, then the free surface can be considered validated. This is explored in the following section.

#### **4. Results and Discussion**

The presentation and analysis of the results generated in this study are split into two major parts. The first relates to the computed ship resistance coefficients, while the second subsection relates to the spectral analysis of ship waves.

#### *4.1. Ship Resistance*

In this subsection, the numerically obtained resistance is briefly discussed and compared to the experimental data. To begin with, the total resistance coefficients are presented in Figure 8 for all cases. The experimental data of Elsherbiny et al. [23] is included alongside each numerical result to enable comparison. Figure 8 clearly indicates the numerical prediction has a well-defined tendency to underpredict the experimental data. Here, the subscripts refer to depth Froude number. As one might expect, the ship's resistance at higher speeds becomes more challenging to predict by CFD. This is evident, especially in the resistance characteristics at the highest depth Froude number for each case.

The overall agreement between the experimental and numerical data is encouraging. This is the first sign that the constructed towing tank is capable of providing good predictions for the resistance of a ship. Possible sources of discrepancy are suspected to stem from the fact that the numerical approach did not model sinkage and trim. In shallow waters, their combined effect, termed ship squat, is attributed greater relative importance than in deep waters. Thus, the imperfect modelling adopted may have been partly the cause of the observed levels of discrepancy between the experimental and numerical results. Moreover, as the ship speed is increased, the difference also grows. This matches the pattern observed in the sinkage and trim curves for a ship, both theoretically [52,53] as well as numerically [54]. Turbulence modelling is also a source of error in the results presented herein. Based on previous research, it is expected and indeed observed that the *k*–ω model has a tendency to underpredict resistance [31]. More sophisticated modelling approaches, such as those accounting for transition or those resolving large parts of the turbulent kinetic energy spectrum, should be used to assess the sensitivity of the adopted approach to turbulence modelling. Recent studies have demonstrated that the latter approach can be highly beneficial in terms of accuracy [55].

**Figure 8.** Comparison of total resistance coefficients for all cases (R indicates the rectangular canal case, whereas S indicates the Suez Canal). Subscripts refer to depth Froude number. Relative error computed as (EFD − CFD)/EFD × 100.

In practice, the cell numbers used in this study tend towards being prohibitively high. As stated earlier, a virtual towing tank where the principle of Galilean invariance has been utilised will consist of no more than 2–3 million cells. Such simulations can be performed within a few days on a standard computer. Thus, the adopted approach of virtually towing the ship may not become widespread soon. Nevertheless, the additional information that may be extracted from a case such as this can be useful. An example of this is given in the following subsection.

### *4.2. Spectral Analysis of the Numerical Free Surface*

We begin by examining the highest speed, simulated in the virtual towing tank, *Fh* = 0.77 in the rectangular canal. The numerical free surface is depicted in Figure 9. Here, the far-field waves generated by the ship are clearly visible. Due to the lateral restriction, waves have reflected approximately 2.5 ship lengths aft of the ship. It should be noted that Figure 9 and other figures are reflected around the central symmetry plane to enable a better visualisation.

**Figure 9.** Generated wave field in the rectangular canal at *Fh* = 0.77.

At this point, it is useful to attempt to determine the wave angle. According to Havelock's [17] method, the Kelvin half-angle is θ = 21.58◦. Figure 10 depicts an attempt at solving this problem by projecting a (dashed) straight line from the forward perpendicular (FP) to the sides (and its reflection) at an angle of 21.58◦. The line 'lands' at a wave trough on the canal wall—clearly, this is not the correct approach. The solid line in Figure 10 represents the same process but in this case beginning from the nearest peak downstream at the wall and projecting in both directions. The line intersects the ship approximately 1/4 *L* from the forward perpendicular. Then, the broken line is initiated at the highest peak at the wall, where wave reflection occurs. This intersects the ship approximately 3/4 *L* from the FP. Finally, the dotted line shows the same process. It originates at the point of the ship defined by [min(*x*), max(*y*)], representing the point where the aft wave system is generated.

**Figure 10.** Generated wave field in the rectangular canal at *Fh* = 0.77 and corresponding half-angle according to Havelock [17]. Dashed line originates at FP; solid line originates at the nearest downstream peak, where the dashed line is reflected. Broken line originates at the highest wave elevation on the wall; dotted line originates at the ship coordinates representing the point where the aft wave system is generated.

Clearly, neither line in Figure 10 accounts for the wave angle well. This is not surprising since the method used to estimate the half-angle is linear and devised for a point source. In this case, the spectral representation may be used to approach the problem. As explained earlier, the first step is to calculate the Fourier transform of the wave field. This is performed in MATLAB, which uses grayscale images. For this reason, the images used henceforth to represent the free surface will be shown in the grayscale format used to perform the analysis. By doing so, other researchers may cross-reference results obtained herein by analysing the provided images.

Figure 11 depicts the process used in this study. Specifically, the image is first represented in Fourier space. Then, for each column of the matrix defining the Fourier transform, a maximum is identified. The *kx, ky* components of these maxima are then compared to the theoretical relationship provided by Equation (8). A polynomial fit is constructed from these points to demonstrate the accuracy of large and small *kx* on the fit. In the present case, it is apparent that as one progresses in the *kx* range to higher values, agreement deteriorates quickly. According to Nakos and Sclavounos [51], insufficient grid resolution will be manifest as numerical dispersion in the *kx, ky* plane being curved towards high values of *kx*, eventually forming closed curves. Since this is not what is observed in the present study, it may be concluded that the numerical wave field is represented with sufficient grid resolution. Nonlinear phenomena will be revealed in the appearance of additional branches in the spectrum. In each quadrant, a branch, emanating from the origin and propagating linearly to the edge of the plot, where it reflects, is observed. Moreover, smaller branches of the dispersion relation are observed, with origins at higher *kx* values, indicating nonlinearity [21]. These lead to deformation of the wake in the real space.

**Figure 11.** Processing of the wave field. Depicted: *Fh* = 0.77 in the rectangular canal. Top: the raw image, where real space extents are 32 m in the stream-wise and 4.6 m in the span-wise directions. Bottom left: the Fourier representation of the wave field. Bottom right: detected maxima and fit, superimposed onto the theoretical relationship, Equation (8).

Now, it is important to deduce the origin of the apparent disagreement in the high *kx* region, as well as its effect on the predicted Kelvin half-angle. One source of disagreement inevitably stems from the fact that the dispersion relation used for comparison is linear [56]. Ship waves, particularly in shallow water, are nonlinear. Even in deep waters, Ma et al. [57] found significant nonlinear influence on the dynamic pressure of the KCS. On the other hand, the experimental data presented in Caplier et al. [20] suggest that the theory approximates real ship waves very well for the Wigley hull. This may not be a fair comparison, because nonlinear effects for the Wigley hull are known to be small [57–59]. In other words, using the parabolic hull plays to the method's strengths.

The neglect of nonlinear terms is chiefly manifest in the near-field disturbance, close to the ship. Although a modification to the Kelvin half-angle may be produced as a consequence of modified pressure in the near-field, the magnitude of such an effect is not known. Interference between transverse and divergent waves, generated at the ship's bow, may be one cause of the observed disagreement [60]. This interference effect, along with the nonlinearity exhibited by the KCS, and the influence of viscosity are all thought to be the contributing to the observed discrepancy.

There is one more aspect of the solution that one should consider carefully. This relates to the curve fit used to approximate the numerical dispersion for higher *kx* values than maxima were detected for. In shallow waters, the arms of the locus are typically not well developed [20]. Thus, it is difficult to extract a sufficient number of points to perform the analysis. For this reason, the only fair assessment recommended is within the range where maxima have been detected from the Fourier transform. The range *kx*/(*g*/*U*2)<sup>∈</sup> [1, 2.5] is used to perform all subsequent analysis. This includes part of the fit over which no maxima have been detected to illustrate the effect of limited data samples.

Figure 12 contains the derivatives d*ky*/d*kx* for deep and shallow water, alongside the numerically generated fit from CFD. This is shown because upon evaluating d*ky*/d*kx* at the inflection point, the Kelvin half-angle can be obtained. Moreover, Nakos and Sclavounos [51] recommend the examination of these derivatives to highlight differences between numerical and theoretical dispersion relations. Figure 12 also includes the area under each curve for reference. Clearly, assessing solely the area under each curve is not a good approach to determine an apparent disagreement, or error, which is −3.901% in this case. This is the case because different parts of the *kx* range over which the derivative is shown may attenuate or reinforce the total favourably. On the other hand, the RMS (root mean square) of the difference between the theoretical and numerical curve is predicted as 0.456. The effect of this on the predicted half-angle is illustrated in Figure 13, where the consequences of the previously examined differences are highlighted.

**Figure 12.** Derivatives d*ky*/d*kx* for deep water, shallow water and numerical shallow water cases.

**Figure 13.** Predicted and theoretical half-angles.

The net effect of the difference between the fit and theoretical curves is translated into a difference of approximately 2.6◦ in the predicted half-angle. This is not considered as a substantial discrepancy. However, the locations where inflection occurs are significantly different between the two sets of data. This occurs at *kx,theory*/(*g*/*U*2) = 1.06 according to the theoretical relationship, whereas CFD predicts this at *kx,CFD*/(*g*/*U*2)<sup>≈</sup> 1.46. The identification of this half-angle does not help in visualising the wake better. Plotting a line with origins at the bow with an angle of 24.1◦ causes an intersection with the wall earlier than what is shown in Figure 10. The relative error between the two predictions is 5.93% (arrived at by computing |(θ*theory* − θ*CFD)*/ θ*theory* × 100)|).

There are several aspects of this technique that should be improved. Firstly, the range over which it is acceptable to find maxima of the Fourier transform should be defined. The only way to accomplish this is via an extensive experimental campaign. If such an interval is known, then it may be possible to define a metric expressing the degree to which waves are correctly modelled. It may also be possible to link specific parts of the computational free surface with increased error in the representation of ship waves. The only manner in which this can be achieved is via experimental work, which should demonstrate the validity of the assumptions as they apply to waves generated by a ship, rather than point sources. Specifically, this concerns ships causing highly nonlinear flows, thereby generating waves that may not resemble the method's predictions well.

This work proceeds with the next aspect of the solution, which one may obtain via the spectral representation. Specifically, this concerns splitting the near-field from the far-field components. This is illustrated for the rectangular case study at *Fh* = 0.77 in Figure 14, where the cut-off wave number is *kc <sup>x</sup>* = 4.7885. Here, the shape of the ship leaves a small effect onto the corresponding near- and far-field systems because the outline of the vessel forms part of the free surface itself. The intensity of the spectrum depends on the input and can be changed based on the brightness of the supplied input. The range of the spectrum is therefore not shown.

**Figure 14.** Splitting of near- and far-field components via manipulations of the spectrum (cut-off wave number *k<sup>c</sup> <sup>x</sup>* = 4.7885). Top: original free surface; (**a**) indicates the far-field component, whereas (**b**) indicates the near-field disturbance and their corresponding Fourier representations. Longitudinal extent: 32 m.

The near-field disturbance is not confined in the immediate vicinity of the ship in Figure 14, contrary to expectations. Instead, it is shed from the ship downstream, with its influence being clearly visible near the domain walls. This representation also allows the detection of the far-field waves, as well as their reflection from the side walls, with ease. It is apparent that the wave system is convex with respect to the ship centreline. Once reflected, this is not as clearly visible. The full spectrum for this case can be consulted in Figure 11.

Figure 15 shows the spectral decomposition process as applied to the rectangular canal for *Fh* = 0.57. Here, the arms of the spectrum, previously used to extract maxima and compare with the theoretical relationship, are not formed. This is consistent with findings of other researchers [20]. Specifically, the higher the depth Froude number is, the more clearly the arms of the spectrum are formed. The physical origins of this relate to the relatively small far-field disturbance generated by the ship in the examined speed range. Simultaneously, speeds corresponding to *Fh* ≥ 1 are impractical, which is why they have not been investigated. Experimental data in terms of resistance are also not

available for the adopted case studies at the aforementioned speeds. In any case, several features of the spectrum can be observed. Firstly, the low-intensity arms, observed in Figures 11 and 14, extending into the far field are reproduced. However, in Figure 15, the arms are not reflected from the boundaries of the plot; instead, they exhibit periodic structures, which vanish near the limits.

**Figure 15.** Free surface disturbance in the rectangular canal, *Fh* = 0.57; (**a**) Computed free surface. (**b**) Far-field and (**c**) near-field representations in the real and spectral space (*k<sup>c</sup> <sup>x</sup>* = 9.5796). Longitudinal extent: 16.5 m.

It is now appropriate to shift the focus onto the Suez Canal and the spectral representation of the wave field obtained. As before, the higher speed (*Fh* = 0.57) is examined first, as shown in Figure 16. An immediately apparent difference relates to the structure of the wave field. Specifically, the slopped canal banks have caused a rundown of the water surface. Since the theoretical relationship used to plot the solid line in the plot can only account for a single depth, it is not seen to represent the Fourier representation of the numerical wave field well. Interestingly, the spectrum contains maxima arranged in semi-circular arcs. An interpretation of this is not attempted at present; instead, this is left for a more theoretical piece of work. Such research would need to determine the form of the dispersion relation in nonconstant water depths. The cut-off wave number used to produce Figures 15 and 16 is the same (*k<sup>c</sup> <sup>x</sup>* = 9.5796 m<sup>−</sup>1) for this reason. As was the case for *Fh* = 0.77, the near-field disturbance is trapped near the canal walls and shed over a great distance downstream.

**Figure 16.** Free surface disturbance in the Suez Canal, *Fh* = 0.57; (**a**) Computed free surface. (**b**) Far-field and (**c**) near-field representations in the real and spectral space (*k<sup>c</sup> <sup>x</sup>* = 9.5796). Longitudinal extent: 13 m.

The lower speed investigated in the Suez Canal and the Fourier representation of its wave field are depicted in Figure 17. As expected, the disturbance generated by the ship at *Fh* = 0.47 is significantly smaller than that produced at *Fh* = 0.57 (Figure 16). The spectrum exhibits a similar structure to what was previously observed for *Fh* = 0.57. For both results shown in Figures 16 and 17, the near-field hydrodynamic response has caused bright parts of the spectrum, which are periodically broken. These correspond to what is identified as a near-field wave by the method, trapped at the lateral extents of the tank. This is primarily the case due to the relative size of these disturbances. Namely, their wavelength is of the order of magnitude of the ship itself. Whether this classification itself is correct probably requires further research. However, their effect on the Fourier representation is clearly visible, especially in Figure 16, where a high-intensity patch can be seen undulating along the ordinate.

**Figure 17.** (**a**) Computed free surface in the Suez Canal, *Fh* = 0.47. Far-field (**b**) and near-field (**c**) representations in the real and spectral space (*k<sup>c</sup> <sup>x</sup>* = 14.1833). Longitudinal extent: 20 m.

#### **5. Summary and Future Work**

This paper presented a towing tank which does not make use of on the principle of Galilean relativity. This was accomplished via the overset method, which was used to actually 'tow' the KCS model in a virtual environment. The main benefits of adopting such an approach were identified in terms of the reduced number of assumptions needed to perform the analysis. Specifically, all boundaries (except those prescribed as symmetries or overset boundaries) are no-slip walls, which is more physically consistent than 'traditional' virtual towing tanks. Since the ship advances over a static fluid, the approach presented in this paper does not require the definition of inlet turbulent properties, as is usually the case. Thereby, one major source of modelling error and uncertainty is removed.

The adopted case studies numerically replicated recently published results in a rectangular canal and in the New Suez Canal [23]. The computed resistance of the ship was compared to the experimentally obtained values. Satisfactory agreement was found, although some discrepancy persisted in the highest speeds examined. The source of the difference between the experimental and numerical results is primarily attributed to the fixed sinkage and trim used in the virtual towing tank. The study was supplemented by a recently developed method to decompose the wave field and determine the Kelvin wake angle. In terms of the former, it was discovered that near-field disturbances propagate outwards towards the canal sides and are shed by the ship downstream. Their effect persisted over a significant distance. This is of practical interest, because near-field disturbances are typically linked with strong pressure variations. Thus, information extracted via the spectral decomposition method may be used to assess the optimum slope and positioning of canal sides to avoid excessive forces linked with bank erosion.

On the other hand, it was shown that it is difficult to identify the boundaries of the Kelvin wake in a narrow canal. Values for the half-angle computed via linear point-source methods were compared with those obtained by CFD. The effects of nonlinearity and interference of wave systems shed by the bow and stern were identified potential sources of discrepancy. The numerical (θ*CFD* = 24.1◦) and theoretical (θ*Theory* = 22.756◦) Kelvin wake half-angles for *Fh* = 0.77 in the rectangular canal were found to compare reasonably well. However, the inflection point, which governs the value of the half-angle, was found to be in some disagreement. According to the theory, this should occur at approximately *kx* = 1.06, whereas CFD suggests the inflection point is located at *kx* = 1.46. The potential sources of this discrepancy likely pertain not only to limitations in terms of mesh size and time-step in CFD, but also, and likely more importantly, the theoretical assumptions.

One of the main conclusions drawn from this study is that the theoretical method used to predict the locus location has great potential for validation of numerical free surface problems. However, bounds of validity are yet to be determined. This forms the primary recommendation for future work. Specifically, a method should be devised to equip the practitioner with a number indicating how well the problem at hand is modelled. The answer to such a question can be found by performing a systematic experimental campaign, which should feature a variety of hull forms. These must ensure that the method is challenged with nonlinear forms to rigorously test its applicability. For example, these would include ships with high block coefficients, nonslender forms, as well as immersed transoms, where separation would be expected.

Another aspect requiring further attention is the manner in which the ship is accelerated up to the target velocity. In this study, the velocity changes linearly, inducing shocks onto the ship at the initial and target speeds. This results in longer times for convergence. To reduce the time taken to converge the result, different acceleration profiles must be compared, such as sinusoidal velocity changes. This is left as a piece of future work, alongside a turbulence dependence study and the relative impact of sinkage and trim.

**Author Contributions:** Conceptualization, M.T. and G.Z.; methodology, M.T. and G.Z.; software, M.T.; validation, M.T. and G.Z.; formal analysis, M.T. and G.Z.; investigation, M.T. and G.Z.; resources, T.T., Z.Y., A.I.; data curation, M.T. and G.Z.; writing—original draft preparation, M.T.; writing—review and editing, all authors; visualization, M.T. and G.Z.; supervision, T.T., Z.Y., A.I.; project administration, T.T., Z.Y., A.I.; All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Acknowledgments:** Results were obtained using the ARCHIE-WeSt High-Performance Computer (www.archiewest.ac.uk) based at the University of Strathclyde. The work reported in this paper is drawn from the first author's PhD thesis. The first author gratefully acknowledges the scholarship provided by the Faculty of Engineering at the University of Strathclyde, which fully supports his PhD.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


*J. Mar. Sci. Eng.* **2020**, *8*, 258


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).
