Numerical Simulations for Determining Resistance, Trim, and Sinkage Characteristics of a Tugboat

8.16. × 105

Fig. 1. Geometry of the tugboat. (a) Lines plan, (b) CAD plan.

physical insights into the hydrodynamic behavior of these specialized vessels, thereby supporting future design and analysis efforts.

2.1. Geometry and simulation cases

The vessel under investigation is a conventional azimuth stern drive (ASD) tug, designed for operations requiring high maneuverability and bollard pull capacity. The hull form design was initially published by Brandner (1995), and scaled model experiments were subsequently conducted by Zheng et al. (2023).

The main particulars of the full-scale and model-scale tugboat are presented in Table 1. A three-dimensional (3D) model of the tugboat was created using its twodimensional (2D) lines plan. The geometry of the model is presented in Fig. 1. Numerical simulations with a model-scale tugboat were performed at seven different flow speeds. Table 2 summarizes all cases. Here, V represents the velocity of the tugboat and Reynolds number (Re) and Froude Number (F r) calculated from the following equations: Re =

Here ρ and µ are the density and dynamic viscosity of water, respectively, while g is the acceleration due to gravity.

2.2. Mathematical model

The viscous flow around the tugboat hull was modeled by solving the three-dimensional, incompressible, Unsteady Reynolds-Averaged Navier–Stokes (URANS) equations. The continuity and momentum equations are given in Cartesian tensor form as: ∂ (ρui ) =0 ∂xi

∂ (ρui ) ∂ ∂ p̄ + ρui u j + = − ∂t ∂x j ∂x j     ∂u j ∂ ∂ui 0 0 + µ + − ρui u j (4) ∂x j ∂x j ∂xi where ui and p̄ are the mean velocity and pressure components, ρ is the fluid density, µ is the dynamic viscosity, and ρu0i u0j represents the Reynolds stress tensor (Ferziger et al., 2019). The Reynolds stress terms were closed using the two-equation Shear Stress Transport (SST) k-ω turbulence model (Menter, 1994). The model solves transport equations for the turbulent kinetic energy k and the specific dissipation rate ω:   ∂k ∂ ∂k ∂k + uj = pk − β ∗ kw + (5) (υ + σk υT ) ∂t ∂x j ∂x j ∂x j   ∂w ∂w ∂ ∂w + uj = αS2 − βw2 + (υ + σw υT ) ∂t ∂x j ∂x j ∂x j + 2 (1 − F1 ) σw2

where pk is the production rate of turbulent kinetic energy, υT is the turbulent eddy viscosity, and F1 is

a blending function that activates the standard k-ω model near the wall and transitions to the k- model in the free stream. The constants α, β, β ∗ , σk , σw are model coefficients (Menter, 1994). This formulation is particularly suited for capturing complex flow separation and adverse pressure gradients present in the stern and skeg regions of the tugboat (ITTC, 2011). The interface between the water and air was captured using the Volume of Fluid (VOF) method (Hirt & Nichols, 1981). The method introduces a scalar volume fraction, α, defined as the fraction of the cell volume occupied by the primary phase (water). The transport equation for α is: ∂αi ∂αi + uj =0 ∂t ∂x j

The volume fraction transport equation was solved using the High-Resolution Interface Capturing (HRIC) scheme (Muzaferija, 1998) to maintain a sharp interface while mitigating numerical diffusion. The fluid properties in each cell are calculated as a weighted average based on the volume fraction: ρ = αρwater + (1 − α) ρair

To predict the vessel’s dynamic response, the hull was allowed to move freely in heave and pitch (trim) using a Dynamic Fluid Body Interaction (DFBI) model (CD-ADAPCO, 2017). Pressure-velocity coupling was achieved using the SIMPLE (Semi-Implicit Method for Pressure-Linked Equations) algorithm of Patankar and Spalding (1983). The governing equations were discretized using a cell-centered finite volume method (FVM). Second-order upwind schemes were employed for the convection terms of momentum, turbulence, and volume fraction equations, while a second-order implicit scheme was used for temporal discretization. All numerical simulations were performed using the commercial CFD software Siemens Star-CCM+.

2.3. Computational domain

The computational domain was designed to accurately capture the far-field flow while ensuring minimal influence from the boundaries on the near-hull solution. The computational domain was designed to accurately capture the far-field flow while ensuring minimal influence from the boundaries on the near-hull solution. The extents of the domain were selected following established guidelines for marine CFD applications, particularly those recommended by the International Towing Tank Conference (ITTC,

2011). These guidelines aim to minimize numerical wave reflection and blockage effects, ensuring that boundary conditions do not artificially influence key results such as resistance, trim, sinkage, and wave patterns. The chosen dimensions are consistent with common practice for simulations of displacement hulls at model scale (e.g., Wilson et al., 2001; Carrica et al., 2007). As shown in Fig. 2, the domain extends 2 × LW L upstream of the bow and 3 × LW L downstream of the stern. Side boundary is located 2 × LW L from the centerplane, while the top and bottom boundaries are positioned 1 × LW L above and 2 × LW L below from the still waterplane, respectively. The boundary conditions were assigned as follows: a free-stream velocity condition was applied at the inlet (upstream), side, top, and bottom boundaries; a pressure outlet condition was specified at the downstream boundary; and a no-slip wall condition was enforced on the hull surface. The symmetry plane was defined along the vessel’s centerline. The volume mesh was generated using the unstructured, hexahedral-dominant trimmed-cell mesher within Siemens Star-CCM+. The background cell size was set to 0.32 × LW L . A systematic, multi-level refinement strategy was then implemented to resolve the critical flow features. The hull surface was discretized with a target cell size of 0.005 × LW L to adequately capture its curvature. Given the hydrodynamic importance of the stern appendages, the skeg surface received an additional local refinement, reducing the target size to 0.00125 × LW L to accurately resolve its sharp geometry and the associated complex flow separation. To sharply capture the air-water interface, a box-shaped control volume encompassing the expected free-surface deformation was defined, extending from the upstream to the downstream and laterally to the domain sides. Within this region, the

mesh was refined in the vertical direction to a target cell size of 0.00125 × LW L .To capture the wake field, a triangular volumetric refinement region was created. This region originates at the bow, extends

1.2. × LW L downstream from the stern, and employs a

target cell size of 0.01 × LW L . For boundary layer resolution, a prism layer mesh consisting of seven layers with a geometric stretch factor of 1.2 was applied on the wetted hull surface. Crucially, for each simulation, the total thickness of this prism layer stack was adjusted to achieve an average target wall unit of y+ ≈ 50 at the first cell centroid across the hull. This value is within the recommended range (30–100) for the use of wall functions with the “All y+ wall treatment” in StarCCM+. This approach represents a computationally efficient and validated methodology for predicting integral forces and moments on ship hulls (ITTC, 2011), providing a practical compromise that avoids the prohibitive cost of maintaining a y+ < 5 across the entire hull under dynamic motion. The final mesh for the high-fidelity case contained approximately 3.0 million cells, which is shown in Fig. 3. The resulting wall y+ distributions for the extreme Froude numbers, confirming the attainment of the target, are presented in Fig. 4.

2.4. Computational details

All numerical simulations were performed using the commercial CFD software Siemens Star-CCM+ (version 2022.1). The computations were conducted on a personal computer equipped with a 6-core processor and 32 GB of RAM. The simulations were run in parallel, utilizing all available computational cores. For the production cases (fine mesh, 1t = 0.010 s), each simulation required approximately 96 hours of

Fig. 3. Volume mesh: (a) Perspective view, (b) Mesh on the hull surface, (c) yz-plane cut at amidships, (d) xz-plane cut at center plane, (e). xy-plane cut on the still water plane.

Fig. 4. Wall y+ distribution on the hull: (a) Fr = 0.106, (b) Fr = 0.375.

wall-clock time to achieve a fully developed, statistically steady solution for resistance, trim, and sinkage. Fig. 5 shows convergence history of CT at Fr = 0.375.

3.1. Verification methodology

RG = The verification study was conducted to estimate the numerical uncertainties stemming from spatial (grid) and temporal discretization. The procedure adhered to established best practices in computational fluid dynamics for ship hydrodynamics, including the guidelines of the International Towing Tank Con-

The convergence is classified as: • Monotonic Convergence: 0 < RG < 1 • Oscillatory Convergence: RG < 0 • Divergence: 1 < RG

Fig. 5. Convergence history of resistance coefficient CT at Fr = 0.375.

For cases exhibiting monotonic convergence, the generalized Richardson extrapolation method is used to estimate the observed order of convergence, pa , 21 extrapolated value, φext , factor of safety, F S, and numerical uncertainty, UG , based on ITTC (2024). In this method, the apparent order of accuracy, pa , can be calculated iteratively from the following equations (Celik et al., 2008): pa =

where r21 and r32 are the refinement ratios of consecutive meshes or time steps. ε21 and ε32 are defined as: (14)

Finally, the approximate relative error, the extrapolated relative error, and the GCI of the fine grid are calculated from (Celik et al., 2008): e21 a =

Where UG = GCI 21 f ine , and F S is the factor of safety, which is calculated based on Xing and Stern (2010): 

2.45. − 0.85P → 0 < P ≤ 1

FS = (20) 16.4P − 14.8 → 1<P Here, P is the ratio of the theoretical order of accuracy to the apparent order of accuracy. When oscillatory convergence is achieved, numerical uncertainty, UG , is determined based on Stern et al. (2001): UG =

3.2. Mesh and time-step convergence results

The total resistance coefficient, CT a primary parameter of interest, is defined as: RT 0.5ρSV 2

where RT is the total resistance, and S is the wetted surface area of the hull, as given in Table 1. First, mesh convergence was evaluated using three systematically refined grids: coarse (721,782 cells), medium (1,485,666 cells), and fine (3,068,731 cells). Table 3 shows the results for grid convergence analysis. According to the results, monotonic convergence

Coarse 721,782 0.1027 0.739 21.4 Medium 1,485,666 0.1043 0.757 20.1 Fine 3,068,731 0.1054 0.769 21.0 Convergence condition Monotonic Monotonic Oscillatory pa 1.574 1.701 N/A FS 1.781 1.727 N/A UG (%) 4.0 5.3 3.1

Table 4. Time-step convergence results at Fr = 0.288. Simulation 1t (s)

Coarse 0.020 0.1051 0.774 Medium 0.015 0.1056 0.771 Fine 0.010 0.1054 0.769 Convergence condition Oscilatory Monotonic pa N/A 2.21 FS N/A 3.301 UT (%) 0.2 0.6

is achieved for CT and trim, while oscillatory convergence was achieved for sinkage. The uncertainty related to grid size in CT , trim and Sinkage are calculated as 4.0%, 5.3% and 3.1%, respectively. Consequently, the fine mesh was selected for all subsequent simulations to ensure the highest practical spatial accuracy. Subsequently, time-step convergence was investigated to ensure temporal independence. The selection of time-step sizes was guided by the ITTC (2011) recommended practice, which suggests a nondimensional time step (1t ∗ ) in the range of 0.005 to 0.01, defined as 1t ∗ = 1t (U/L), where U is the flow velocity and L is the characteristic length (taken as LW L ). For the convergence study case at Fr = 0.288 (V = 1.632 m/s and LW L = 3.272 m), this corresponds to a dimensional time-step range of approximately 0.010 s to 0.020 s. Accordingly, simulations were performed with three different time-step sizes within this range: 1t = 0.020 s, 0.015 s, and 0.010 s. Table 4 shows the results for time-step convergence analysis, where UT is the numerical uncertainty due to the time-step. The results show that the uncertainty related to time step in CT , trim and Sinkage are below 1%. Based on this analysis, a time step of 1t = 0.010 s was selected for all subsequent cases.

4.1. Resistance, trim and sinkage

The calm-water hydrodynamic performance of the tugboat, characterized by its total resistance coefficient CT , trim angle, and sinkage at the center of gravity, was numerically investigated across the

Froude number range of 0.106 to 0.375. The results from the fine-mesh, fine-time-step simulations are compared with the experimental data published by Zheng et al. (2023) in Fig. 6. Fig. 6(a) presents the variation of the total resistance coefficient CT with Froude number. The numerical predictions closely follow the experimental trend, successfully capturing the characteristic resistance curve of a conventional displacement hull. The agreement is quantified by a mean absolute percentage error (MAPE) of 2.9% across the speed range, with the largest individual discrepancies of 5.8% occurring at the lowest speed (Fr = 0.106), and 6.9% occurring at the highest speed (Fr = 0.375). The comparison of trim angle is shown in Fig. 6(b). The numerical model reproduces the trend measured experimentally, with the trim angle increasing with Froude number as expected. The computational results show consistent alignment with the physical measurements, with a MAPE of 6.0% for trim. The maximum error was 9.3% at the highest speed (Fr = 0.375), where wave-making effects and dynamic motions are most pronounced. Sinkage at the center of gravity is presented in Fig. 6(c). The sinkage increases with Froude number due to dynamic pressure changes along the hull. The computational fluid dynamics (CFD) results show good agreement with the experimental measurements, with a MAPE of 7.2% for sinkage. It should be noted that the relative error at the lowest speed (Fr = 0.106) is 27.3%; however, this corresponds to an absolute discrepancy of less than 0.5 mm in model scale, reflecting the challenge of accurately measuring and predicting very small sinkage values near the static condition. Excluding this lowest speed point, the MAPE for sinkage reduces to 3.8%. It is noted that the comparison error for trim and sinkage does not vary monotonically with Froude number but shows fluctuation. This behavior is physically consistent and reflects the changing dominance of underlying hydrodynamic mechanisms. At lower speeds (Fr < 0.2), performance is governed by viscous pressure fields and boundary layer development, which are well-captured by the RANS model. At higher speeds, dynamic wave-making pressures— characterized by nonlinear interactions between the bluff bow, skeg, and stern wave systems—become the primary driver for trim and sinkage. The tugboat’s specific geometry, with its high displacement-tolength ratio and prominent stern appendages, makes this transition particularly sensitive, leading to local fluctuations in error as the flow regime shifts. Overall, the numerical simulations demonstrate a consistent and reliable ability to predict the key calm-water performance metrics of the tugboat. The

Fig. 6. Comparison of the experimental (EFD) and numerical (CFD) results for (a) resistance coefficient, (b) trim, and (c) sinkage. Experimental results are from Zheng et al., 2023.

agreement with the experimental data from Zheng et al. (2023) validates the selected numerical methodology, including the URANS approach with the SST k-ω turbulence model, the VOF free-surface treatment, and the DFBI motion model, as a practical tool for the hydrodynamic assessment of this type of vessel.

4.2. Flow field around tugboat

A detailed analysis of the flow field around the tugboat hull was conducted to elucidate the hydrodynamic mechanisms governing resistance, trim, and sinkage across the investigated Froude number range. Three key characteristics were the main parameters of interest: the evolution of the free-surface wave pattern, the hull pressure and skin-friction distributions, and the wake structure at the propeller plane. Together, these elements are critical for a comprehensive understanding of ASD tugboat performance in calm water.

Fig. 7 presents the computed free-surface elevations around the hull for selected Froude numbers. At the lowest speed condition (Fr = 0.106), the wave system is relatively weak and symmetric, consisting of a modest divergent bow wave and a gently recovering stern wave. The free surface remains largely undisturbed along the midship region, indicating that viscous resistance dominates the total resistance at this speed. As the Froude number increases, wave-making effects become increasingly significant. Distinct bow and shoulder wave systems develop, accompanied by a pronounced wave trough near the bow and an amplified stern wave depression. At higher speeds (Fr ≥ 0.288), strong divergent and transverse wave components are clearly observed, with increased wave steepness and downstream propagation. These features are characteristic of displacement-type hulls operating at moderate to high Froude numbers and directly correspond to the observed rise in the total resistance coefficient reported in Section 4.1.

Fig. 7. Free surface deformations around the tugboat at various Froude numbers.

Fig. 8. Dynamic pressure distribution on the tugboat: (a) Fr = 0.106, (b) Fr = 0.26, (c) Fr = 0.375.

Fig. 9. Skin friction distribution on the tugboat: (a) Fr = 0.106, (b) Fr = 0.26, (c) Fr = 0.375.

The dynamic pressure distributions on the hull surface for representative speeds are illustrated in Fig. 8. At low Froude number, the pressure field is largely governed by hydrostatic effects, with elevated pressures concentrated around the bow stagnation region and a relatively smooth pressure gradient along the hull. As speed increases, the influence of dynamic pressure becomes more pronounced. High-pressure regions intensify at the forebody due to increased stagnation effects, while extended low-pressure zones develop along the aft shoulder and stern regions. Notably, the skeg and stern appendage areas experience localized pressure peaks, particularly at Fr = 0.26 and Fr = 0.375, reflecting strong flow acceleration and deceleration in these regions. These pressure asymmetries contribute significantly to the longitudinal force balance and play an important role in the predicted trim behavior of the tugboat. Fig. 9 shows the distribution of skin-friction coefficient on the hull surface. At lower speeds, the skin-friction distribution is relatively uniform, with higher values concentrated near the bow and forward bottom region due to thinner boundary layers and higher local velocities. With increasing Froude number, the magnitude of skin friction increases over the entire wetted surface, especially along the flat bottom and stern regions where flow acceleration and turbulence intensity are enhanced. Elevated skin-friction levels are also observed near the skeg and stern areas, indicating strong viscous effects associated with flow curvature and possible local separation. The wake structure downstream of the hull is examined through normalized streamwise velocity contours at the propeller plane (x/LW L = 0.11, where x is measured from stern), as shown in Fig. 10. At Fr = 0.106, the wake is relatively narrow, with a modest velocity deficit centered along the hull centerline. As the speed increases, the wake expands both vertically and laterally, and the velocity deficit becomes more pronounced. At higher Froude numbers, complex wake patterns emerge, characterized by strong velocity gradients, low-speed regions behind the skeg, and accelerated flow zones near the outer portions of the propeller disk. These wake features are of particular relevance for ASD tugboats, as they directly influence inflow conditions of the propulsors, propulsion efficiency, and maneuvering performance. Overall, the flow field analysis demonstrates that the numerical model successfully captures the essential hydrodynamic features associated with tugboat operation across a wide speed range. The predicted free-surface patterns, pressure and shear distributions, and wake characteristics are physically consistent and provide valuable insight into the resistance mechanisms and dynamic behavior of the vessel.

Fig. 10. Normalized streamwise velocity distributions at propeller plane: (a) Fr = 0.106, (b) Fr = 0.26, (c) Fr = 0.375.

This detailed flow-field information complements the global performance metrics presented earlier and further confirms the robustness of the adopted URANS–VOF–DFBI modeling framework for tugboat hydrodynamic analysis.