Uncertainty Analysis and Maneuver Simulation of Standard Ship Model

Abstract

Maneuver simulation of a standard ship model gives indication of numerical accuracy. In the numerical calculation of ship maneuvering, uncertainty analysis is a necessary step to ensure the accuracy of the calculation. In this study, uncertainty pair analysis is carried out in the simulation of the turning circle motion of the standard ship model ONRT in waves. According to the uncertainty analysis procedure recommended by the International Towing Tank Conference (ITTC), the change of ship resistance caused by the number of grids is studied to determine the influence of grid density on the numerical prediction. The simulation of turning motion in waves is carried out based on the uncertainty analysis. It is found that the minimum number of overset grids for this simulation is 1.4 million. The numerical results are fairly accurate compared to experimental results, and this technique provides a method with low calculated cost for this simulation.

1. Introduction

The maneuverability of a ship is directly related to its safety during navigation. Many international standards have been proposed for ship maneuvering issues, which provide regulations and safety recommendations for maneuverability indicators of different ship types and propose new research directions [1,2,3]. This indicates that ship maneuvering remains one of the most important research areas in the field of ocean engineering [4,5]. The prediction of ship maneuverability can be achieved through free-running model tests, which involve using scaled models of the hull, propellers, and rudders to conduct standard maneuvering tests in a physical maneuvering basin. The International Towing Tank Conference (ITTC) has developed corresponding specifications for ship maneuvering tests [6,7], and various classification societies have provided guidelines for maneuvering tests [8]. However, this method requires a large testing basin, precise control systems for propellers and rudders, and six-degree-of-freedom motion measurement devices for the model, making it expensive. Therefore, using mathematical models for simulating ship maneuverability has become a more focused-on approach. To evaluate the accuracy of numerical calculations in predicting ship maneuverability, ITTC provides experimental data for specified tests conducted in a water basin, which can be used for comparing with computational fluid dynamics (CFD) results, along with recommended procedures for CFD uncertainty analysis [9].

Simulation methodologies for ship maneuvering primarily encompass three approaches: direct prediction methods, system-based research methods, and CFD (Computational Fluid Dynamics) numerical simulation techniques. CFD numerical simulation methods can be divided into constrained ship model and self-propelled ship model simulations. With the advancement of computational power, CFD techniques have achieved significant breakthroughs, providing robust support for the numerical calculation of ship maneuvers under substantial interactions among the hull, propeller, and rudder. Due to their excellent intuitiveness and accuracy, CFD methods are widely employed within the industry. CFD techniques are capable of simulating various maneuvering motions for self-propelled vessels, forecasting the vessel’s trajectory and typical maneuvering parameters, such as overshoot angle, advance, and transfer, which directly reflect the ship’s maneuverability performance. Despite the complexity of simulating a self-propelled ship, this remains one of the most popular research topics due to the numerous advantages it offers.

To obtain more accurate self-propulsion maneuvering results and detailed flow field information, it is necessary to implement the coupled motion solution of the hull, propeller, and rudder. The overset grid method allows for the treatment of complex motion coupling between the ship, propeller, and rudder without causing grid deformation during the calculation process, ensuring the accuracy of the computation. Carrica and Stern [10] integrated the hydrodynamic solver CFD Ship-Iowa, which is based on the overset grid method, to achieve numerical simulation of the self-propelled ship’s maneuvering motion. They directly constructed discrete models of the hull, propeller, and rudder based on geometric bodies, using the overset grid method to handle the coupled motion of ship, propeller, and rudder during self-propulsion maneuvering operations. They conducted numerical simulations of the zigzag maneuver motion and turning circle motion for the KVLCC ship type. Shen et al. [11] use the solver naoe-FOAM-SJTU, which is developed based on OpenFOAM, and obtained results that are very close to the experimental values. Mofidi and Carrica [12] used the same solution and numerical methods to simulate the typical 10/10 zigzag maneuver test and the modified 15/1 zigzag maneuver test numerically. The predicted ship motion and maneuverability parameters were in good agreement with the experimental results, and they analyzed the detailed flow field during the self-propulsion maneuvering process. Broglia et al. [13] and Dubbioso et al. [14] respectively conducted numerical simulations of turning circle tests for twin-propeller ships with single and twin rudders. The motion of the rudder and hull was treated using a dynamic overset grid. The numerically simulated ship’s trajectory was compared with the experiment, and the turning circle trajectories under single and twin rudder conditions were compared, along with the reduction of the propeller speed, drift angle, and roll angle over time. They analyzed the changes in rudder force and lateral force on the hull and appendages throughout the turning motion, noting that under twin-propeller conditions, the presence of rudders would significantly interfere with the loads experienced by the propellers. Muscari et al. [15] also employed the same approach, using an in-house developed CFD solver with the overset grid method to handle the motion of the hull and rudder. They numerically simulated the pure yaw and turning circle maneuvering motions of twin-propeller, twin-rudder vessels.

CFD analysis of the maneuvering of a self-propelled ship requires a large computational resource that restricts its use as part of a ship design process. Researchers have been working on reducing the computational cost [16], and this could aid CFD practitioners to make their work easier. To find fewer grids for the calculation, uncertainty analysis is a good method to apply.

Uncertainty is a concept derived from metrology. Physical experiments or measurements require quantitative statements to assess the quality of the measurement results, i.e., to provide quantitative determination of the quality of the measurement results. The usefulness of physical experiments or measurement results largely depends on the magnitude of their uncertainties. In other words, measurement results must be accompanied by uncertainty statements to be complete and meaningful. Commonly used generalized probabilistic methods for uncertainty analysis include robust Bayesian theory [17], interval description [18], possibility theory [19], etc. The Dempster–Shafer theory of evidence (DSTE) proposed by Dempster provides practical and effective methods for handling cognitive uncertainty and has been widely used [20,21,22]. In CFD simulation of ship maneuvering motion, it is necessary to perform uncertainty analysis on the CFD simulation results to assess their credibility. It is generally believed that the numerical uncertainty of CFD calculations originates from three sources: truncation error, iteration error, and discretization error. In ship hydrodynamic calculations, the first one can be neglected, while the second one is easy to reduce to negligible level. Discretization error is the most important part of the error and is generated in the process of converting continuous partial differential equations into algebraic equation systems. It is primarily determined by the numerical discretization scheme and grid used for the simulation. Together with physical modelization error, they constitute the two most significant parts of the error in practical simulations of complex turbulent flows.

This study is based on the ITTC recommended procedures and focuses on CFD uncertainty analysis of the ONRT model, and is used to figure out the minimum grid number needed to simulate the standard ship movement with high precision.

2. Governing Equations

It is well known that the motion of fluid within a flow field generally needs to satisfy several conditions: the law of mass conservation, the law of momentum conservation, and the law of energy conservation. During numerical simulations, the governing equations are used to mathematically express the aforementioned conservation laws. The governing equations are represented by the continuity equation (mass conservation equation) and the Navier–Stokes (N-S) equations (momentum conservation equations), which are expressed as follows:

$${\displaystyle \frac{\partial U_i}{\partial x_i}}=0$$

(1)

$$\rho {\displaystyle \frac{\partial U_i}{\partial t}}+\rho {\displaystyle \frac{\partial \left(U_i{U}_j\right)}{\partial x_i}}=-{\displaystyle \frac{\partial P}{\partial x_i}}+\rho {\displaystyle \frac{\partial }{\partial x_j}}\left[\nu \left({\displaystyle \frac{\partial U_i}{\partial x_j}}+{\displaystyle \frac{\partial U_j}{\partial x_i}}\right)\right]+\rho g_i$$

(2)

In the equations, U_i represents the velocity component in the direction of the x_i coordinate axis; P denotes pressure; ρ signifies fluid density; and g_i symbolizes the gravity force. For the time-averaged treatment of the Navier–Stokes (N-S) equations, the Reynolds Averaged Navier–Stokes (RANS) approach is predominantly utilized. The RANS equations are capable of decomposing turbulent flows into mean flow terms and fluctuating flow terms. The transient variables within the equations are separated into time-averaged quantities and fluctuating quantities, as illustrated in the following expressions:

$$\phi =\overline{\phi }+{\phi }^{\prime }$$

(3)

In the equations: $\overline{\phi }$ represents the time-averaged quantity and ${\phi }^{\prime }$ denotes the fluctuating quantity. Consequently, the time-averaged versions of Equations (1) and (2) can be expressed as Equation (4) and Equation (5), respectively:

$${\displaystyle \frac{\partial U_i}{\partial x_i}}=0$$

(4)

$$\rho {\displaystyle \frac{\partial U_i}{\partial t}}+\rho {\displaystyle \frac{\partial \left(U_i{U}_j\right)}{\partial x_i}}=-{\displaystyle \frac{\partial P}{\partial x_i}}+\rho {\displaystyle \frac{\partial }{\partial x_j}}\left[\nu \left({\displaystyle \frac{\partial U_i}{\partial x_j}}+{\displaystyle \frac{\partial U_j}{\partial x_i}}\right)-\rho \overline{{u^{\prime }}_i{u^{\prime }}_j}\right]+\rho g_i$$

(5)

In the equations: For convenience of expression, the time-averaged quantity $\overline{\phi }$ is often denoted by $\phi$ and the Reynolds stress is represented by $\overline{{u^{\prime }}_i{u^{\prime }}_j}$, which is considered as an unknown term.

The present study employs the Finite Volume Method (FVM) for the discretization of the governing equations. The Finite Volume Method, as the name suggests, involves the integration of the governing partial differential equations within a finite volume, with the conservation equations expressed as Equation (6):

$${\displaystyle \underset{\Delta V}{\int }\frac{\partial \left(\rho \varphi \right)}{\partial t}}dV+{\displaystyle \underset{\Delta V}{\int }div}\left(\rho \overrightarrow{u}\overline{\varphi }\right)dV={\displaystyle \underset{\Delta V}{\int }div}\left(\Gamma grad\varphi \right)dV+{\displaystyle \underset{\Delta V}{\int }S}dV$$

(6)

In the formulation: $\varphi$ signifies a generic variable; $\Gamma$ represents the generalized diffusion coefficient; and S denotes the generalized source term. Using the Gauss divergence theorem, Equation (6) is transformed into the following expression:

$${\displaystyle \underset{\Delta V}{\int }\frac{\partial }{\partial t}}\left(\rho \varphi \right)dV+{\displaystyle \underset{A}{\int }\overrightarrow{n}}\left(\rho \overrightarrow{u}\overline{\varphi }\right)dA={\displaystyle \underset{A}{\int }\overrightarrow{n}}\left(\Gamma grad\varphi \right)dA+{\displaystyle \underset{\Delta V}{\int }SdV}$$

(7)

The essence of solving the discretized equations involves selecting appropriate discretization schemes for the governing equations. In accordance with the coupling principle between pressure P and velocity V, the solution methodology adopted in this article primarily utilizes the SIMPLE (Semi-Implicit Method for Pressure-Linked Equations) algorithm. The core of the SIMPLE algorithm is to iteratively guess and correct the solution using the basic grid, ultimately solving the momentum equation. In addition, there is the SIMPLEC (Semi-Implicit Method for Pressure-Linked Equations, Consistent) algorithm, which is very similar to the SIMPLE algorithm. However, in the SIMPLEC algorithm, the flux correction approach is further refined, allowing for better convergence. The PISO (Pressure-Implicit with Splitting of Operators) algorithm, in contrast to the previous two methods, introduces an additional correction term, achieving a secondary improvement and significantly enhancing the convergence speed within a single iteration step.

The Reynolds stress modelized by Equation (8):

$$-\rho \overline{{u^{\prime }}_i{u^{\prime }}_j}={\mu }_t\left({\displaystyle \frac{\partial u_i}{\partial x_j}}+{\displaystyle \frac{\partial u_j}{\partial x_i}}\right)-{\displaystyle \frac{2}{3}}\left(\rho k+{\mu }_t{\displaystyle \frac{\partial u_k}{\partial x_k}}\right){\delta }_{ij}$$

(8)

In the equation: ${\mu }_t$ represents the turbulent viscosity; $u_i$ denotes the time-averaged velocity; and k signifies the turbulent kinetic energy.

In this study, k-ω SST model [23] is selected as the turbulence model for the numerical simulations of this paper. The k-ω SST model combines the characteristics of both k-ε and k-ω, thereby employing k-ε in the free-stream regions while utilizing k-ω in the near-wall regions. The definition of the turbulent viscosity is presented as Equation (9):

$${\mu }_t={\displaystyle \frac{\rho k}{\omega }}{\displaystyle \frac{1}{\max \left(\frac{1}{{\alpha }^{\ast }},\frac{S{F}_2}{{\alpha }_1\omega }\right)}}$$

(9)

Equations (10) and (11) represent the transport equations corresponding to k and ω, respectively:

$${\displaystyle \frac{\partial \left(\rho k\right)}{\partial t}}+{\displaystyle \frac{\partial \left(\rho k{u}_i\right)}{\partial x_i}}={\displaystyle \frac{\partial }{\partial x_j}}\left[{\Gamma }_k{\displaystyle \frac{\partial k}{\partial xj}}\right]+{\tilde{G}}_k-Y_k$$

(10)

$${\displaystyle \frac{\partial \left(\rho \omega \right)}{\partial t}}+{\displaystyle \frac{\partial \left(\rho \omega u_i\right)}{\partial x_i}}={\displaystyle \frac{\partial }{\partial x_j}}\left[{\Gamma }_{\omega }{\displaystyle \frac{\partial \omega }{\partial x_j}}\right]+G_{\omega }-Y_{\omega }+D\omega$$

(11)

In equations: ${\Gamma }_k=\mu +{\displaystyle \frac{{\mu }_i}{{\sigma }_k}}$${\Gamma }_{\omega }=\mu +{\displaystyle \frac{{\mu }_i}{{\sigma }_{\omega }}}$; ${\tilde{G}}_k=\min \left(G_k,10\rho {\beta }^{\ast }k\omega \right)$; $G_k={\mu }_t{S}^2$; $G_{\omega }={\displaystyle \frac{\alpha }{{\nu }_i}}{G}_k$.

3. CFD Uncertainty Analysis Method

The error of numerical simulation is the difference between the calculated value and the true value, and it consists of two parts: modelization error and numerical error. Under certain conditions, the error value can be estimated. However, since the true value is often unknown, there is also an error in this estimation process: uncertainty is an estimate of the error. Due to the unknown true value, the error is often unknown as well, and the magnitude of the error does not change with people’s level of understanding. Uncertainty represents the degree of lack of knowledge about the error. Uncertainty is objectively present, but the results of uncertainty analysis will change with people’s level of understanding.

It is generally believed that the numerical uncertainty of CFD calculations originates from three sources: truncation error, iteration error, and discretization error. In ship hydrodynamic calculations, the first error can be neglected, while the second error is easy to reduce to a negligible level. Discretization error is generated in the process of converting continuous partial differential equations into algebraic equation systems, primarily determined by the numerical simulation grid, and constitutes the most significant part of the error with modelization in practical simulations of complex turbulent flows.

According to the ITTC recommended procedures [24], the process of CFD uncertainty analysis can be divided into two parts: verification and validation. The verification process evaluates the numerical uncertainty and estimates the magnitude and sign of the numerical error, along with its uncertainty. The verification process is essentially an assessment of whether the equations have been correctly solved. The validation process evaluates the model uncertainty and estimates the magnitude and sign of the model error, along with its uncertainty. The confirmation process is essentially an assessment of whether the correct equations have been solved, i.e., whether the mathematical model has been established correctly.

The definitions of error and uncertainty are the same as those in experimental uncertainty analysis. The difference between the numerical simulation result S and the true value T is the numerical simulation error ${\delta }_S$, which consists of the model error ${\delta }_{SM}$ and numerical error ${\delta }_{SN}$:

$${\delta }_S=S-T={\delta }_{SM}+{\delta }_{SN}$$

(12)

For specific cases, the sign and magnitude of the numerical error ${\delta }_{SN}$ can be estimated as:

$${\delta }_{SN}={\delta }_{SN}^{\ast }+{\epsilon }_{SN}$$

(13)

where ${\delta }_{SN}^{\ast }$ is the estimated value of ${\delta }_{SN}$ (including sign and magnitude) and ${\epsilon }_{SN}$ is the error in the estimation.

A more accurate extrapolated result $S_C$ can be obtained with:

$$S_C=S-{\delta }_{SN}^{\ast }$$

(14)

Verification is the process of calculating the numerical uncertainty $U_{SN}$ of the numerical simulation. If conditions permit, estimating the sign and magnitude of the numerical error ${\delta }_{SN}^{\ast }$ itself and the uncertainty $U_{S_{C}N}$ in this error estimation. For unmodified numerical simulation methods, the numerical error can be decomposed into errors ${\delta }_1$, ${\delta }_G$, ${\delta }_T$, and ${\delta }_P$, which come from the iteration count, grid size, time step, and other parameters, respectively. Therefore, the numerical uncertainty of the numerical simulation can be expressed as:

$${U_{SN}}^2={U_1}^2+{U_G}^2+{U_T}^2+{U_P}^2$$

(15)

For modified numerical simulation methods, the solution is modified, leading to the numerical benchmark $S_C$. The estimated values ${\delta }_{SN}^{\ast }$ and $U_{S_{C}N}$ of the simulated numerical error are given by the following equations:

$${\delta }_{SN}^{\ast }={\delta }_1^{\ast }+{\delta }_G^{\ast }+{\delta }_T^{\ast }+{\delta }_P^{\ast }$$

(16)

$${U_{S_{C}N}}^2={U_{I_C}}^2+{U_{G_C}}^2+{U_{T_C}}^2+{U_{P_C}}^2$$

(17)

Validation is the process of using the extrapolated result to evaluate the model uncertainty $U_{SM}$ of the numerical simulation, and, if conditions permit, estimate the sign and magnitude of the model error ${\delta }_{SM}$ itself. The error comparison E is given by the difference between the test data D and the simulation value S:

$$E=D-S={\delta }_D-\left({\delta }_{SM}+{\delta }_{SN}\right)$$

(18)

By comparing E with the validation uncertainty $U_V$, the validation can be assessed.

$${U_V}^2={U_D}^2+{U_{SN}}^2$$

(19)

If $\left|E\right|<U_V$, all the combinations of errors in D and S are less than $U_V$, and the validation at this level is achieved. If $U_V<<\left|E\right|$, the model can be improved using the sign and magnitude of $E\approx {\delta }_{SM}$. For modified numerical simulation methods, the corresponding equations are:

$$E_C=D-S_C={\delta }_D-\left({\delta }_{SM}+{\epsilon }_{SN}\right)$$

(20)

$${U_{V_C}}^2={U_{E_C}}^2+{U_{SM}}^2={U_D}^2+{U_{S_{C}N}}^2$$

(21)

4. CFD Uncertainty Analysis Results

The verification and validation methods and procedures for CFD uncertainty analysis in this study are based on the ITTC Quality Manual, and the variables’ definitions also follow the manual. The ONRT ship model used in the simulation of ship motion under calm water and wave conditions in this study is consistent. It is a full-appended model with twin propellers and twin rudders. This model is widely used for CFD validation and is listed as a standard model for the free-running model problem at the Tokyo 2015 CFD Workshop. There are abundant maneuvering test data available for this model, making it suitable for verifying the reliability of current numerical prediction methods. Optimizing this model simulation in terms of different parameters has been a popular issue in ship research, and the data are becoming more abundant [25,26]. The geometric model of the hull is shown in Figure 1, and the main dimensions of the hull are listed in

Table 1. Parameters of ONRT.

Parameters	Symbol	Modal Scale	Full Scale
Length at the Waterline (m)	$L_{WL}$	3.147	154.0
Beam (m)	$B_{WL}$	0.384	18.78
Depth (m)	$D$	0.266	14.5
Draft (m)	$T$	0.112	5.494
Wetted surface (m2)	$S_0$	1.5	NA
Block coefficient (CB)	$\nabla /\left(L_{WL}{B}_{WL}T\right)$	0.535	0.535
Displacement (kg)	$\Delta$	72.6	8.507 × 106
Vertical position of CG (m)	$KG$	0.156	NA
Buoyant to bow (m)	$LCBaftFP$	1.625	NA
Moment of inertia	$K_{xx}/B_{WL}$	0.444	0.444
	$K_{yy}/L_{WL}$, $K_{zz}/L_{WL}$	0.25	0.25
Propeller diameter (m)	$D_P$	0.1066	NA
Max steering (deg/s)		35.0	NA

.

This study is based on the STAR-CCM+ (18.02) software and primarily focuses on numerical computational analysis of typical free-running motion under wave conditions. The waves in the calculations are in the initial head sea state. Detailed calculation conditions can be found in

Table 2. Numerical conditions for simulation.

Wave Length $\lambda$ (m)	3.147
$\lambda /L_{WL}$	1.0
Wave Height $H$ (m)	0.06294
Wave Steepness $H/\lambda$	0.02
Rudder Angle (°)	35°

. In the numerical simulation, the ship navigates with full degrees of freedom at the model’s self-propulsion point. The ship’s speed is U = 1.11 m/s, corresponding to Fr = 0.20.

When conducting calculations for the design speed Fr = 0.20, the grid uncertainty analysis for this speed was first performed. The numerical uncertainty analysis used here follows the method recommended by ITTC. Three sets of grids were used for convergence verification, ranging from 0.56 × 10⁶ to 1.95 × 10⁶ grid cells. The grid scaling factor was $\sqrt{2}$, y⁺ = 11.63, and the surface grids of the ship with different sizes are shown in Figure 2.

To achieve this scaling effect, during the grid generation process, the scale of the background grid was uniformly scaled in the X, Y, and Z directions according to this ratio. Then, relative sizes were used when locally refining the grid using grid partitioning tools. This ensured that the majority of the generated grid maintained this proportional relationship. However, the automated process of generating purely unstructured grides involves complex procedures like expanding boundary layers, which may not fully meet the required proportional criteria in certain detailed areas. Nevertheless, this does not affect the overall grid uncertainty analysis. To minimize uncertainty factors, the grid uncertainty analysis is performed with simulations in which all degrees of freedom of the ship are fixed except the forward motion.

The calculation employed overset grid techniques. During computations, the computational domain is treated as an entirety, with the governing equations being discretized through the Finite Volume Method and solved within the temporal domain using a segregated solver for the unsteady solution. The numerical simulation of the six degrees of freedom motion of the ship model on the free surface is realized by the simultaneous solution of the ship body motion and surrounding flow field. In the numerical solution process, the resultant forces and moments acting on the ship model include viscous forces, pressure, and the ship model’s gravity force, all of which stem from the interaction between the ship model and the flow field. As shown by the computational flow chart below, the actual computation is a process of continuous coupling and iteration of flow field information. Initially, the surrounding flow field of the ship model is calculated to determine changes in shear stress and pressure. Subsequently, the resultant forces and moments are obtained from the external force equations, which are then input into the six degrees of freedom motion differential equations for integration to yield the ship model’s displacement. The grid node positions are then updated based on the ship model’s displacement to solve the new flow field around the ship model. Factors such as water viscosity, turbulence intensity, free surface deformation, gravity, and atmospheric pressure are considered during the solution process. Through continuous coupling and iterative computation, the external forces and the ship model’s posture are obtained in real time, ultimately enabling the simulation of the ship model’s dynamic response. The overset grid diagram is shown in Figure 3.

The main verification object of grid uncertainty analysis is the dimensionless resistance coefficients, including the total resistance coefficient ($C_T$), pressure resistance coefficient ($C_P$), and viscous resistance coefficient ($C_V$), which are defined as follows:

$$\left[\begin{array}{c}{C}_T\\ C_V\\ C_P\end{array}\right]=\left[\begin{array}{c}{\displaystyle \frac{R_T}{0.5\rho U_0^{2}T{L}_{PP}}}\\ {\displaystyle \frac{R_V}{0.5\rho U_0^{2}T{L}_{PP}}}\\ {\displaystyle \frac{R_P}{0.5\rho U_0^{2}T{L}_{PP}}}\end{array}\right]$$

(22)

where $R_T$ is the total resistance experienced by the model, $R_P$ is the pressure resistance, which is the integral of the normal forces, and $R_V$ is the viscous resistance, which is the integral of the tangential forces. ρ represents the density of water, $U_0$ is the ship’s velocity, T represents the ship’s draft, and $L_{PP}$ is the length between perpendiculars of the ship.

According to the ITTC guidelines for verification and validation (V&V) of CFD numerical calculations, the determination of convergence forms in grid uncertainty analysis is based on the convergence parameter $R_G$, whose value is determined by different density grid solutions (generally including three density grid: coarse, medium, and fine) denoted as $S_i$. It is defined as follows:

$$R_G={\displaystyle \frac{S_2-S_1}{S_3-S_2}}$$

(23)

where the subscript $S_i$, with i = 1, 2, 3, represents the results for fine grid, medium grid, and coarse grid, respectively. Different values of $R_G$ correspond to different convergence forms:

$$\begin{array}{ccc}\left(i\right)& 0<R_G<1& Monotonic convergence\\ \left(ii\right)& 0<\left|R_G\right|<1,R_G<0& Oscillatory convergence\\ \left(iii\right)& \left|R_G\right|>1& Divergence\end{array}$$

(24)

For the first case, which is consistent convergence, the grid uncertainty $U_G$ is typically evaluated using Richardson extrapolation (RE). For the second case of oscillatory convergence, the grid uncertainty is determined by the average difference between the maximum value $S_U$ and the minimum value $S_L$ of the oscillations, i.e.,

$$U_G=1/2(S_G-S_L)$$

(25)

In the third case, it is proven that convergence has not been achieved, and therefore grid uncertainty is not applicable. The resistance parameters for different grids can be seen in

Table 3. Grid uncertainty analysis of ONRT at design speed.

Grid	Symbol	Number	${\mathit{C}}_{\mathit{P}}$	${\mathit{C}}_{\mathit{V}}$	${\mathit{C}}_{\mathit{T}}$	${\mathit{C}}_{\mathit{T}}$ Error
Experiment					1.706 × 10−2
fine	S1	1.95 M	0.507 × 10−2	1.171 × 10−2	1.678 × 10−2	−1.66%
medium	S2	1.40 M	0.499 × 10−2	1.162 × 10−2	1.661 × 10−2	−2.64%
coarse	S3	0.56 M	0.437 × 10−2	1.124 × 10−2	1.597 × 10−2	−6.39%
$R_G$			0.307	0.237	0.265
$P$			3.4088	4.1559	3.834
$GC{I}_{12}$			0.87%	0.31%	0.45%
$GC{I}_{23}$			2.88%	1.27%	1.74%
Convergence			uniform	uniform	uniform

. From the table, it can be observed that as the grid is refined, the numerical predictions tend to converge consistently with the experimental values. Therefore, it is necessary to use the Richardson extrapolation (RE) method to evaluate grid uncertainty. The order of discretization $P$ is defined by the following equation:

$$P={\displaystyle \frac{In\left(1/R_G\right)}{In\left(r\right)}}$$

(26)

The Grid Convergence Index (GCI) can be obtained using the following equation:

$$GC{I}_{ij}=F_S{\displaystyle \frac{\left|e_{ij}\right|}{r^p-1}}$$

(27)

where $F_S$ is a safety parameter, which is typically set to 1.25 when using three or more grids. $e_{ij}$ represents the difference between $S_i$ and $S_j$. The Grid Convergence Index (GCI) can characterize the impact of different grids on the numerical results, and a smaller GCI value indicates low sensitivity of the numerical results to the grid.

From Table 3, it can be seen that all the resistance coefficients exhibit consistent convergence. The GCI₁₂ (difference between fine grid and medium grid) for the total resistance coefficient is only 0.45%, and the GCI₂₃ is consistently larger than the former. This indicates that after reaching a medium grid density, the numerical predictions are less affected by changes in the grid. To improve computational efficiency while maintaining accuracy, it is recommended to use the medium grid for simulation.

After using the medium grid for numerical calculations, this section provides a direct numerical validation for the free turning maneuvering motion with a rudder angle of 35 degrees. The free turning trial is the most commonly used method for assessing the maneuverability of a ship. During this trial, while the ship is cruising in a straight line, the rudder is turned to a certain angle—typically the maximum angle of 35 degrees—and held constant, resulting in the ship entering a turning motion. The entire turning motion is generally divided into three phases: the initial rudder turning phase, which extends from the commencement of the rudder turn command until the specified rudder angle is reached; the transition phase, which spans from the cessation of rudder turning until the ship enters a steady state of turn; and finally, the steady-state phase, where the ship’s motion parameters become stable and the vessel enters a new stage of equilibrium.

The numerical simulation of the free turning maneuvering trial still focuses on the fully appended ONRT ship model, which is equipped with twin propellers, twin rudders, and appendages such as struts, deadwood, and bilge keels. The ship’s parameters and principal dimensions are presented in Table 1. In the numerical simulation, the propeller rotation speed is set to the self-propulsion point of 8.97 RPS, which was determined from a self-propulsion computation, and the ship’s speed corresponds to a Froude number (Fr) of 0.20. The numerical simulation encompasses the process of the ship model completing one full cycle of free turning.

The numerical calculation of the free-running maneuver test is still carried out using a high-performance parallel computing cluster. During the numerical calculation of the free-running maneuver test, a total of 36 processes are used for parallel computing, among which 35 are for fluid field calculation and 1 is for interpolation calculation. The time step is set to $\Delta t=0.005$ s. The propeller is simulated by using the volume force method, which is equivalent to propulsion provided by the propeller.

The computed trajectory of the ship’s turning motion at a rudder angle of 35 degrees and its comparison with experimental values are shown in Figure 4.

Figure 4 demonstrates that the numerically predicted trajectory of the free-running turning motion of the ship is largely consistent with the turning circle measured in the experiment. The specific parameters of the turning circle from numerical predictions show good agreement with the experimental values, with the numerically calculated turning radius being slightly larger than that obtained from the experiment. To quantitatively analyze the precision of the current numerical prediction,

Table 4. Comparison of Gyroscopic Characteristic Parameters.

Characteristic Parameters	CFD	EFD	Error
Advance (m)	6.963	6.998	0.50%
Transfer (m)	3.872	3.880	0.21%
T90 (s)	11.582	11.570	0.10%
T180 (s)	22.653	22.410	1.08%
Tactical Diameter (m)	9.784	9.646	1.43%

presents a comparison of the characteristic parameters of the free-running turn against the experimental values.

This fully demonstrates that the current numerical methods can simulate the free-running maneuvering motion with high precision. To examine whether the attitude of the ship is safe or not, the freedom of motion was investigated. Figure 5, Figure 6 and Figure 7 present the heave, pitch, and roll of the ship versus time. It can be observed from the figures that there is significant fluctuation in the curves at the initial moment due to the presence of the rudder turn phase and the transition phase. The curves of the numerically predicted motions of the ship during free turning in waves are given. From the figure, it can be seen that the ship’s heave, pitch, and roll motions oscillated with wave frequency. In addition, because the ship’s encounter wave direction is changing all the time during the free turning, there are also low-frequency fluctuations due to the free turning maneuvering under the high-frequency motion. The maximum amplitude of heave is 0.06 m, the maximum amplitude of pitch is between −2.3° to 2.2°, and the amplitude of roll is between −3.8° and 7° during the whole turning motion. In addition, it can be seen from Figure 7 that the amplitude of the roll motion due to the waves is larger than that due to the initial maneuvering.

Figure 8 presents the free surface wave profile around the ship during the steady-state turning phase. In contrast to the towing conditions, a noticeable difference in the free surface on both sides of the ship can be observed here. Due to the presence of the turning motion, the wave height at the bow on one side of the ship is significantly different from the other side, and there is also a strong asymmetry in the wave pattern at the stern.

After conducting numerical calculations using a medium grid, the CFD results closely match the experimental results. Figure 4 illustrates the comparison between the predicted ship motion trajectory obtained from numerical simulations in waves and the corresponding experimental values. It can be observed from the figure that the current numerical predictions are in good agreement with the experimental results. Additionally, it can be seen that as the heading angle of the ship ranges from 270° to 360°, the turning curve exhibits a noticeable expansion. The numerical calculation yields a turning radius significantly smaller than the experimental value, primarily due to modifications made to the geometric model of the rudder in order to ensure sufficient interpolation between overlapping grid cells in the numerical simulation. This modification reduced the effective rudder area, resulting in a decrease in rudder effectiveness. The current numerical calculations can accurately predict the maneuvering characteristics of a ship freely turning in waves with high precision.

5. Conclusions

In this study, the standard ship model ONRT was modeled, and computational simulations were conducted. The influence of grid refinement on ship resistance was analyzed to assess the impact of grid uncertainty on numerical calculations. It was found that using a medium grid (1.4 million cells) for simulation yielded results that closely matched the standard experimental values, accurately reflecting the maneuvering performance of the actual ship during free-running motion in waves. Therefore, it is recommended to use STAR-CCM+ with overset grides for ship maneuvering simulations and the grid should not be less than 1.4 million cells.

Following the numerical methods in this paper may provide a way to make the calculating cost lower for some practitioners who do not have very good calculating resources. The future work will focus on improving the accuracy with same numbers of cell, as there is still room for improvement in the rationalization of the grid.