Investigation of the ship–seabed interaction with a high-fidelity CFD approach

Ship–seabed interaction is highly critical to ship safety and ship performance when ship operates over shallow water with uneven seabed. However, the prediction method for ship reaction is still yet to be fully developed. In this paper, a high-fidelity transient numerical simulation method with sliding mesh is developed based on computational fluid dynamics to model a vessel passing a step bank, demonstrating to be a computational cost economical solution. The comprehensive numerical model is validated and verified against benchmarked experimental and numerical studies, which is proved to be highly accurate in predicting force characteristics and wave development. Meanwhile, during the research the impulse effect generated by the step bank was found to have striking effects on the wave elevation, wake development and vessel sinkage. Five regions on the behaviour of vessel sinkage are defined in the present work. According to the result, the vessel will encounter the extreme sinkage after a relatively long distance (12 Lpp\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${L}_{\mathrm{pp}}$$\end{document} at Fh2=0.519\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${F}_{h2}=0.519$$\end{document}) passing the step bank instead of immediately.


Introduction
When ships operate along waterways or in coastal regions, while in close proximity to the shore, the ship has been known to experience a significant motion as squat motion [1]. This motion could cause ships unexpectedly grounded, damaging the hull. On the other hand, under harsh sea conditions, ships struggling to manoeuvre in shallow waters could also cause ship grounding. And squat effect is known to play an adverse role. However, with most ships designed to sail in deep water, the seabed effect onto the ships is often ignored. This becomes critical when the ships enter shallow water areas.
During the operation in the shallow water, the ship starts to feel the seabed, encounter larger hydrodynamic forces, lose control and lose energy efficiency [2]. To date, the ship-seabed interaction effect is still lack of understanding, such as what is the effect of uneven seabed to ship propulsion performance, what is the effect to ship safety and whether ship-seabed interaction would actually cause ship grounding. These questions need to be carefully clarified. To perform such research and to generate the best practice to guide ship operation in real seas, either experimental study or numerical investigation needs to be conducted. Due to the facility limitation, most of tests use the false bottoms/step bank, which often are deemed to be too short. The length effect of the false bottom/step bank could not be extensively studied by the experimental method. A high-fidelity time resolved numerical simulation tool is desired to demonstrate this fundamental fluid dynamic problem. However, due to the significant computational time and expense using large number of mesh modelling the seabed bathmetry and the ship together, a computational cost economical method has not yet been developed. Under this framework, the present study is aimed to develop a novel transient numerical simulation method based on computational fluid dynamics (CFD) using sliding mesh to enhance the ability to resolve the highfidelity fluid features with a minimum requirement on mesh and hence the computational cost. A classic benchmark case of ship passing a step bank is studied using the developed numerical model to demonstrate the capability. With such approach, predicting the dynamics of the ship within restricted water area can be performed with cost economical and time efficient way, which would support the design of channels, development of safe port operating procedures and assist ship operation trainings [3].
The case study of step bank simulation is also aimed to reveal the interaction between the ship and the seabed. In the early stage of squat effect study, the work conducted by Constantine [1] undertook a study of squat effect in response to an industry-based request. Both experimental and numerical investigations were conducted. In their work, three distinct flow regimes for ships operating in restricted water, namely sub-critical, critical and supercritical, are revealed. In the sub-critical regime, the flow is steady and can be treated as a simple steady state by employing Bernoulli and continuity equations. When entering into the critical range, the flow turns to be unsteady piling up the solitary wave. Finally, in the super-critical regime, the "steady states" are achieved again. Kijima and Nakiri [4] proposed a numerical method to predict the manoeuvrability of ship in both deep and shallow waters at the initial stage of design. Gourlay [5] used the slender-body method to predict the squat effect. Following the theoretical works, a review was carried out by Gourlay [5] summarising the concerns on sinkage and trim of modern container ship in shallow water, two potential flow methods are discussed with some comparisons to the model test results. Apart from the theoretical development, numerous experimental studies have been conducted to investigate the hydrodynamic performance of ship model in shallow water over the last few decades [6][7][8][9]. The capabilities of different ship manoeuvring simulation methods in shallow water were later benchmarked by SIMMAN [10].
On the other hand, the development of CFD has improved modern ship design and simulation. In the last decade, the continued technological advances offer ever-increasing computational power, in which CFD methods are rapidly gaining popularity for simulating ship performance including the performance in shallow water. Saha et al. [11] used the CFD method to improve the hull form in shallow water and pointed out that the sinkage should be considered as a hydrodynamic design constraint in shallow water. Good discussions of the shallow water effects on resistance, forces and moments, form factor, as well as local flow field are provided by Toxopeus [12] based on a serious of CFD calculations on a KVLCC2 tank. Meanwhile, CFD is one of the numerical tools to perform the ship manoeuvring simulation. Apart from the works mentioned above, research on ship manoeuvring using CFD method was developed by the group, Carrica et al., at the University of Iowa. Carrica et al. [13] predicted the heave and pitch motions of the DTMB 5512 model in head waves using the CFDShip-Iowa V4 with the implementation of overset grids. The method was further expanded and used for simulating ship self-propulsion [14], manoeuvring [15] as well as dynamic stability [16]. Based on the works listed above, Carrica et al. [17] also experimentally and numerically investigate the 20/5 zigzag manoeuvre for the container ship KCS in shallow water.
However, most of the investigations on the ship related simulation are still conducted based on uniform water depth condition. This is mainly because of that the tank tests are usually conducted by steady towing of the ship model, either with constant deep water depth or with constant shallow water depth. And it is still highly challenging to perform CFD simulation over uneven seabed to obtain the time resolved solution. Therefore, the study in this paper is initiated to tackle this challenge. Figure 1 presents the sketch of the unsteady problem with some of the principal definitions. L t is the overall length of the bottom. h 1 and h 2 are the water depth in deep and shallow water region, respectively. The vessel is travelling with a constant speed U towards the shallow water region, where d l is the horizontal distance between the vessel's centre of gravity and step bank. To research the unsteady squat motion, Duffy [3] conducted an experimental investigation similar to the setup shown in Fig. 1. This is also being referred as the false bottom test in the towing tank, which is used to investigate ship moving into shallow water region from deep water region. During this event, ship will experience a sudden change in the water depth, which will excite the ship motion in a sudden manner and gradually decays away.

Introduction of the case study
It is a challenge to perform both the experimental study and the numerical investigation. First, to conduct the experimental study, the towing tank facilities are often subjected to the limited configuration in terms of setting the step bank length and height. And due to the long decay period, the ship model often cannot reach the quasi-steady condition before reaching the end of the step bank. In addition, setting different water depths are extremely costly. In terms of the numerical investigation, to simulate such case with the current methodology, the whole length of towing tank needs to be modelled within a computational domain, which is highly expensive computationally.
Giving an example, ITTC [18] recommends that for simulations in the presence of incident waves, the inlet boundary should be located 1-2 L pp away from the hull, whereas the outlet should be positioned 3-5 L pp downstream to avoid any wave reflection from the boundary walls. Therefore, for a traditional ship in head wave simulation with a constant water depth, the computational domain is around 7.5 × 5 L pp (length × width) [19]. However, as an unsteady phenomenon, the computational domain for ship passing over the step bank requires to be significantly expanded to simulate this dynamic process. Thus, the computational cost is significantly increased by the expanded domain. To simulate the ship-bank interaction, the computational domain requires a 33 × 2.3 L pp (length × width) in the present study to obtain the quasi-steady result in both deep and shallow water region, which is more than two times larger than the traditional ship in head wave simulation.
In this paper, a benchmark ship model, KCS, is used to conduct the following studies with a scale factor, 1:31.599, without rudder appended. The principal dimensions of the KCS model [20] are presented in Table 1. A speed of 18 knots (1.646 m/s for the model) is selected to perform the simulation based on the previous tests performed by Enger et al. [21] to validate the simulation accordingly.

Computational methodology and mesh generation
In the present study, an incompressible unsteady Reynolds Averaged Naiver Stokes (URANS) method is employed using the Finite Volume Method CFD package, Star-CCM+ [22]. The volume of fluid (VOF) method is applied to perform the simulation capturing the characteristics of air-water interface. The k-omega two-equation shear stress transport (SST) [23] is solved using a segregated iterative solution method based on SIMPLE-algorithm together with the VOF equation and the conservation equations.

The use of two regions
Unlike the traditional CFD simulation for ships, to provide a cost-effective and mesh efficient solution. Our present numerical model is proposed to have two sub-regions: (1) the seabed region and (2) the ship region. The region (2) moves along with the vessel, whereas the region (1) is fixed on the earth coordinate system as a stationary region. In between, there is a sliding mesh interface between the region (1) and region (2) to exchange the flow information.
In this way, the region (2) is actively extracting the seabed condition surrounding the ships. With such a method, the size of the region (2) can be greatly reduced to only the domain around the ships. While a thin region, region (1) seabed region, can be extended to a large area with small increment on the total mesh number. With the ship motion, the region (2) perform translational motion over the region (1), "scanning" through the seabed condition. Two different coordinate systems are applied, the local coordinate system moving with the ship and the global coordinate system fixed to the earth for the seabed. To replicate the benchmark testing conditions Enger et al. [21] and Gourlay [5], the ship is fixed to perform a steady towing with no free trim and sink allowed. The principle dimensions for each region in this case study are shown in Fig. 2.

Definition of the boundary conditions
The boundary conditions are illustrated in Fig. 2. It is noted that a velocity inlet boundary condition was assigned to the frontal surface, the top surface of region (1) and both sides of both region (1) and (2), with a flat wave condition to define the water level and the wind-wave-current condition. Zero wind and current speed and zero wave height relative to the global coordinate system are used in the simulation to simulate the calm water condition in the towing tank. The negative x-direction boundary is employed as the pressure outlet condition, using the flat wave pressure condition. The hull of the vessel, as well as the seabed, were modelled as a non-slip wall boundary. The interface, where the two regions are connected, is set as an internal interface with uniform mesh control to guarantee the mesh quality. The internal interface takes the responsibility to convey data between two regions [22]. The flow filed was solved as well as the hydrodynamic forces and moments acting the vessel hull.

Mesh generation
Following the setup of the computational domain and boundaries, the mesh was subsequently generated using the Star-CCM+ trimmer meshing technology [22]. Figures 3  and 4 show the mesh around the vessel, where free surface refinements are used to capture the Kelvin wave around the hull. Based on prior experience and the recommendation in Star-CCM+ [22], a minimum of 20 cells was used in the vertical direction, where the free surface wave is.
To guarantee the precision in the data communication between these two regions, the internal interface used in the present study requires an identical and high-resolution mesh at the boundary. Thus, a local refinement is set to guarantee that the interface mesh is finely and uniformly resolved (see Fig. 4).

Mesh and time step sensitivity study
The numerical set-up for the mesh sensitivity study is provided, as shown in Table 2. The M1 case is presented as the coarsest mesh in the sensitivity test. In this case, the computational domain consists a total of 8.4 million elements. Further cases are illustrated in Table 2 ranged from 8.4 million to 21 million elements. As can be seen in Table 2, the M1, M2, and M3 generally increase all the mesh refinements in the computational domain including wall boundaries, wake of the vessel and free surface. The M3-W1 and M3-W2 particularly focus on the refinement of the bottom wall boundary. Due to the present study was carried based on finite water depth, especially focused on the non-uniform seabed effect, the bottom wall boundary has a significant effect on the hydrodynamics of the vessel. Thus, the boundary layer at the bottom wall needs to be further refined to provide accurate simulation results.
The current overall computational domain is separated by a step bank with two different water depths: deep water section and shallow water section. The deep water section has a water depth ( h 1 ) of 12 times the vessel draft (T). And the shallow water section water depth ( h 2 ) is 3 times of the vessel draft in this mesh sensitivity study.
The mesh sensitivity study, based on the effects of mesh refinement, provides results for time-averaged nondimensional total resistance coefficient ( − C T ) and heave force derived midship sinkage ( − S ) in both deep and shallow water sections. C T and S are obtained by the equations as follows:  where F D is the total resistance force, F H is the heave force, is the fresh water density, g is the gravity of earth, A s is the area of wetted surface and A w is the waterplane area.
To obtain a quasi-steady result of − C T and − S ∕L pp , both parameters were obtained by averaging the results from d l = −6.5 ∼ −1.5L pp in the deep water region and d l = 8.5 ∼ 13.5L pp in the shallow water region ( d l = 0L pp indicates the position of the step bank) to avoid the unsteady change at the step bank. It should be noted that all the cases for the mesh sensitivity study used the same time step ( Δt ) of 0.02 s. This is determined based on the related procedures and guidelines of ITTC [18] for resistance computations in calm water (where Δt = 0.005 ∼ 0.01L pp ∕U , U is the vessel speed).
As can be seen in Table 3, the M3 case is considered to be fine enough for the mesh convergence of . The transient section starts as when the vessel is reaching to the step bank ( d l ∕L pp = −1.5 , with d l ∕L pp = 0 being the moment when the centre of gravity of the vessel reached the step bank) and ends when the vessel reaches to a quasi-steady state in the shallow water section ( d l ∕L pp = 5.5 ). By considering the unsteady problem, there is a relatively large difference between the M3 and M3-W1 case for both C T and S∕L pp . However, by further refining the mesh at the bottom wall from the case M3-W1 to M3-W2, the results are very close indicating the mesh convergence has been reached.
Apart from the mesh sensitivity study, the time step sensitivity study was also carried out for the convergence evaluation. The mesh model chosen for initial time step sensitivity study is the case M3-W1, which was then repeated for two further cases with different time steps. As can be seen in Table 4, the M3-W1 case is considered to be fine enough. Therefore, M3-W1 has been chosen for further validation and verification of the numerical model against previous experimental and numerical results.

Comparison with experiments and previous studies
Due to the lack of research on this unsteady problem, the present numerical model is extremely hard to validate and verify. To do the comparison with experiments and previous studies, the present section is separated as two main part: (1) comparison of the deep water section results against the previous experimental and numerical results; and (2) comparision of the shallow water section results against the previous theoretical results.

Comparision of the deep water section results
Experimental data available from Perić et al. [21] is the main reference for validating the present numerical simulation in deep water. The resulting − C T and − S ∕L pp from the present numerical simulation and the experimental measurement from Enger et al. [21] are presented in Table 5. It is noted that the presented total resistance coefficients are well agreed with the experimental results. However, there is a relatively large discrepancy for the sinkage. A 15.9% discrepancy is observed. Therefore, a further comparision study is performed comparing the current numerical results with the CFD results provided by Perić et al. [21] along with their experimental study. As shown in Table 6, a discrepancy of only 6.74% has been observed. Based on both validation and verification studies, the results from the present numerical model show a good agreement in the deep water section.

Comparision on the shallow water section results
Due to the lack of the experimental results available in the shallow water, the present numerical results of non-dimensional sinkage at a range of shallow water depth ( h 2 = 3-5 T) are compared with the numerical predictions based on slender-body theory carried out by Gourlay [5]. As shown in Fig. 7, the results from the present numerical model show a good agreement against previous numerical predictions in the shallow water section.

Total resistance drag force and sinkage
A selection of total resistance and sinkage results are presented for a KCS vessel travelling over a simplified step bank, described in Sect. 3.1. As can be seen in Fig. 8, for all F h2 , a significant increase of the total resistance is observed when the vessel passing over the step bank. Two peaks (P2 and P3) occur when the vessel travelled over a step bank. It is noted that C T at P3 is higher than its value at P2, especially for a larger F h2 value ( C T at P3 is 5% higher than it at P2 when F h2 = 0.519 ). After the vessel passed the step bank, the total resistance dropped rapidly to a value closing to C T at deep water. However, it is observed that, C T increases along with F h2 in the shallow water region.
Apart from C T , the midship sinkage is illustrated in Fig. 9. As seen in the figure, S∕L pp increased from a   Gourlay et al. [6] Present result Fig. 7 Midship sinkage from present simulations as a function of depth Froude number ( F h2 ) compared with the theoretical predictions provided by Gourlay [5] steady-state value when the ship travelled over the step bank. For all F h2 simulated in the present work, S∕L pp are very simlair before the LCG of the vessel reaching to the step bank. After P2, the midship sinkage starts to increase rapidly. It is noted that after the stern of the vessel passed the step bank (P3, indicating the ship has fully entered into the shallow water region), the midship sinkage fluctuates. Four peaks (named A, B, C, and D in Fig. 9) are observed in the current study. To better describe the unsteady progress, five regions are manually separated and defined (named I, II, III, IV and V in Fig. 9). The first peak value of S∕L pp is the lowest peak value in the shallow water region. The midship sinkage then starts to drop until the vessel enters into region II. In region II, S∕L pp is increased significantly till reaching to point B.
It is noted that,S∕L pp at point B is much higher than it at point A. This indicates a dangerous sinkage on the vessel. In region III, S∕L pp at point C is very similar to the value of point B. A huge decrease in S∕L pp is observed at point C′ in region III indicating a recover of the sinkage on the vessel. However, when the vessel reaches region IV, a very large peak of S∕L pp is observed at point D. The value of S∕L pp at point D is higher than any peaks appeared before (6% higher than the value at point B). It is noted that point D is located in the downstream area more than 10 times of L pp from the step bank at F h2 = 0.519. For a smaller F h2 , the peaks in each region is occurred slightly early than a larger F h2 . This novel finding indicates that the vessel will not encounter the extreme sinkage immediately after passing the step bank. The sinkage fluctuates during the travel of the vessel. As mentioned before, due to the limitation of the test facility, the experiments carried out by Duffy [17] only performed 5 L pp after the step bank. However, in the present study, it is showed that a more critical peak has occurred when the vessel travelled 12 L pp after the step bank. The present numerical results show that, a severe sinkage can occur even when the vessel has passed over the step bank after a relatively long distance.

Correlation of wave elevation, wake and sinkage
To have a general visual appreciation of the wave pattern around the ship, the surface wave at F h2 = 0.519 is plotted along with the longitudinal wave cuts in Fig. 10. The unsteady wave patterns may reveal the key factor which induced squat when a vessel passes over a step bank. As can be seen in Fig. 10(I), before the vessel enters into the shallow water region, a clear Kelvin wake is created by the vessel consisting both divergent and transverse waves (named "original waves"). Once the vessel entered into the shallow water region, as shown in Fig. 10(II), the free water surface in front of the vessel is strongly piled up due to the shallow water effect. Similar observations can be found in the longitudinal wave cut along the vessel in Fig. 10. The rise of the free surface can generate a new series of waves (named "shallow water waves") which is different from the original waves. These waves make the vessel behaves like riding on a incoming wave. At P2, the shallow water waves only affect the free surface around the bow of the vessel. When the whole vessel enters into the shallow water region, at P3, the original divergent waves start to be affected by the shallow water waves. Deep troughs occur around the vessel and the divergent waves become flatter and wider. This indicates that the impulse generated by the step bank strats to affect the waves around the ship. However, due to the vessel only travelled a very short distance after the step bank, the transverse waves have not been affected yet. When the vessel reaches point B, the troughs around the vessel become deeper and wider. The impulse effect on the waves around the vessel can be clearly found in the longitudinal wave cuts around the vessel. In addition, the wake of the vessel is influenced by the impulse. The impulse effect weakens the troughs in the transverse waves at point B. Point C and C′, corresponding to the third peak and bottom in the time history of the sinkage, are plotted in Fig. 10(V) and (VI), respectively. As can be seen in the figures, the divergent waves are stronger and larger when the vessel reaches the peak of the sinkage, compared with the waves observed at C′ point. This induces the change in the sinkage. At point C, the shallow water waves already passed over the bow of the vessel. Thus, the wave elevation in front of the vessel's bow drops back to the value close to the original level at P1. In addition, the wake area of the vessel still has not been affected until point C. Marginal difference can be found at point C′. When the vessel reaches point D, the largest peak observed in the present sinkage time history, clear changes can be found in the wake area. As can be found in the longitudinal wave cuts, significant differences can be found on both wavelength and height in the wake area immediately after the vessel. The impulse effect starts to appear in the wake region. It excites the sinkage of the vessel. Additionally, the divergent waves are stronger and larger than them at point C. This phenomenon also enhances the sinkage of the vessel. Apart from the phenomena mentioned above, the divergent waves are separated into two series in the wake area of the vessel, which can be clearly founded in Fig. 10(VII). The corresponding animation for this figure is provided in Electronic Annex I and Electronic Annex II.

Water depth effects on the unsteady process
As can be seen from Fig. 11, all the wave patterns and elevation around the KCS vessel is the same at P1. Once the vessel reached to P2, the wave patterns around the bow of the vessels start to change due to the difference of F h2 (see Fig. 12). With a larger F h2 , a higher wave elevation around the bow of the vessel is observed. However, at P2, the wave patterns in the wake region are still the same under all F h2 . Compared with Fig. 11, it is noted that the wave elevation if front of the ship has been piled up due to the shallow water effect. When the stern of the vessel reached the step bank, not only the wave in front of the vessel but also the wave around the vessel also have significant changes (see Fig. 13). Additionally, it can be observed that the Kelvin wake changes along with the change of F h2 . With a higher F h2 [e.g., F h2 = 0.519, in Fig. 13(I)], the Kelvin wake becomes wider than the other lower F h2 cases. The wave patterns at point B (when the second largest S∕L pp occured), is plotted in Fig. 14. The wake region of the vessel starts to be affected by the change of F h2 . Compared with Fig. 13, the Kelvin wake becomes stronger at a higher F h2 . Strong impulse effect on the waves are observed from the troughs in the Kelvin wake. These troughs become deeper and wider, leading to an increase in both resistance and sinkage. However, for a relatively low F h2 ( F h2 = 0.450 ∼ 0.402 ), the Kelvin wake does not affected significantly. In the present study, the largest sinkage occurred at point D for all F h2 . Figure 15 shows the wave patterns around the vessel at point D. It can be clearly seen that the Kelvin wake becomes stronger than in Fig. 13(I), (II). Deeper troughs have also been observed around the vessel. In addition, the wave in the wake region also varied significantly at a high F h2 , especially at F h = 0.519 . However, as shown in Fig. 15, the wave elevations in front of the vessel are nearly the same for all F h2 simulated in the current study.

Conclusions
This paper introduces the use of sliding mesh to reduce the computational resource dependant when predicting a vessel's reaction passing an uneven seabed. The numerical model is well validated and verified against previous experimental and numerical results. A benchmark process on simulating a vessel passing over a step bank has been demonstrated in the present work. Observations on the interactions between a ship and waterway boundaries are examined in detail to reveal some insights of fluid physics due to the unsteady phenomena. On the basis of the wave characteristics, the mechanism of ship squat is revealed. The present numerical results show that, a severe sinkage can occur even when the vessel has passed over the step bank after a relatively long distance (12 L pp at F h2 = 0.519 ). This novel finding indicates that the vessel will not encounter the extreme sinkage immediately after passing the step bank. Five regions on the behaviour of vessel sinkage are defined in the present work. This also offers a benchmark for future experimental studies on investigating the unsteady ship-seabed interactions.
Correlation has been studied and demonstrated. The impulse effect generated by the step bank found to have the most striking effect on the wave elevation, the wake development and the ship sinkage. The shallow water waves generated by the ship-seabed interaction make the vessel behaves like riding on a wave introducing a strong effect on the Kelvin wake. It is worth noting that a higher water depth Froude number in the shallow water region may have a stronger nonlinear effect on the sinkage. To improve the accuracy of the present numerical method, a further study considering a wider range of F h2 and a larger domain with fully released vessel motion is needed to examine their effects in the numerical model properly.