The Investigation of a Sliding Mesh Model for Hydrodynamic Analysis of a SUBOFF Model in Turbulent Flow Fields

Abstract

A computational fluid dynamics (CFD)-based simulation using a finite volume code for a full-appendage DARPA (Defense Advanced Research Projects Agency) SUBOFF model was investigated with a sliding mesh model in a multi-zone fluid domain. Unsteady Reynolds Averaged Navier–Stokes (URANS) equations were coupled with a Menter’s shear stress transport (SST) k-ω turbulence closure based on the Boussinesq approximation. In order to simulate unsteady motions and capture unsteady interactions, the sliding mesh model was employed to simulate flows in the fluid domain that contains multiple moving zones. The pressure-based solver, semi-implicit method for the pressure linked equations-consistent (SIMPLEC) algorithm was employed for incompressible flows based on the predictor-corrector approach in a segregated manner. After the grid independence test, the numerical simulation was validated by comparison with the published experimental data and other numerical results. In this study, the capability of the CFD simulation with the sliding mesh model was well demonstrated to conduct the straight-line towing tests by analyzing hydrodynamic characteristics, viz. resistance, vorticity, frictional coefficients, and pressure coefficients.

1. Introduction

The relationship between the configuration of a submerged body and fluid friction has been determined by model analysis and ship model resistance tests. Regarding the analysis of the resistance force of a submerged body, Hoerner [1] indicated that when a streamlined body was attached to an additional object, this additional object interfered with the main body and affected the stress distribution in its surroundings. The induction from experiments revealed that when this additional object was placed away from the main body’s maximum diameter, the interference was minimized; by contrast, when this object was placed on the maximum diameter of the main body, the interference was maximized. Arentzen and Mandel [2] indicated that the resistance resulting from the sail of a submerged body is approximately 15%–30% of the resistance experienced by a bare hull. Friedman [3] subsequently obtained a similar result from an experiment, revealing that the resistance caused by a sail is approximately 30% of the total resistance of a submerged body; moreover, the resistance from its appendages accounts for 20%–60% of the overall resistance. For the Defense Advanced Research Projects Agency (DARPA) SUBOFF model, the appendages would cause a mean increase of about 16% in the total resistance Shariati and Mousavizadegan [4]. Joubert [5] stated that future design trends for a submerged body should focus on mitigating resistance, enhancing operating stability, minimizing noise, and optimizing cruising ability and layout. The optimization of shape is a common indicator for evaluating the performance of a submerged body. The nose shape, length, diameter, and aft body shape of the hull and the sail, location, and height of the sail must undergo design optimization. Lin, et al. [6] employed a planar motion mechanism to conduct static and dynamic maneuvering tests for a submarine and obtained the hydrodynamic coefficients of the submerged body.

With the rapid development of computer technology, the feasibility of using the computational fluid dynamics (CFD) model to implement the hydrodynamic simulations of the submerged body and ships gradually becomes higher. The advantage of employing the CFD model is the ability to directly solve nonlinear equations and rapidly determine the designed vehicle shape. However, to render simulation conditions of the CFD model as close to reality as possible, a grid independence test must be conducted, and independent analysis of a fluid domain must be designed, and the convergence must be determined [7]. Additionally, simulation results must be examined and compared with relevant experimental data to establish a complete analytical system. The International Towing Tank Conference also proposed conceptual instructions on the theory and methods for using CFD to conduct a ship model test [8].

Fretes, et al. [9] conducted a CFD simulation and revealed that appendages caused changes in the stress distribution of a bare hull. Similarly, the resistance experienced by the submerged body was affected. According to the simulation results, the appendage resistance was 78.81% of that of the bare hull. Specifically, the sail, fore hydroplane, aft hydroplane, and rudder comprised 51.57%, 18.87%, 4.22%, and 4.15% of resistance, respectively, which was similar to the results of another experiment. To mitigate the influence of appendages on the resistance of a submerged body, Gorski and Coleman [10] used the CFD method and Reynolds-averaged Navier–Stokes (RANS) equations to conduct the simulation and optimal design of sail shape. These results were consistent with those obtained from wind tunnel and tank experiments. Vaz, et al. [11] compared simulations of a DARPA SUBOFF model for two configurations, i.e., one bare hull and one fully appended hull, under different inflow angles by investigating different turbulence models with a RANS approach. However, their simulated results were not realized in the flow conditions of higher Reynolds numbers.

In terms of the turbulence model, Karim, et al. [12] used the SST k-ω turbulence model to conduct a simulation analysis of the viscosity resistance of an axially symmetric underwater vehicle. The study employed both structural and non-structural grids to conduct analyses, and the results were in agreement with the measurements in the experiment. Moreover, the study indicated that simulating the motion of a submerged body using non-structural grids was superior to using structural grids. Gross, et al. [13] employed the CFD method to determine the turbulence model of a Large Eddy Simulation (LES) to simulate the streamline distribution around a SUBOFF bare hull in a steady flow field, where the maximum inlet angle was set at 30°. An experiment using a small submerged body was conducted to verify the simulation. Moonesun, et al. [14] conducted a similar simulation analysis and revealed that using numerical simulation to estimate the resistance of a submerged body in an infinite flow field is more feasible than using conventional empirical equations. By using the CFX software, Hu, et al. [15] employed the k-ε turbulence model to simulate the translational hydrodynamic coefficients, and applied the k-ω model to analyze the rotational and other coupled hydrodynamic coefficients.

There are two basic numerical techniques used for solving flow problems of moving domains, i.e., Moving Reference Frame (MRF) method [16,17] and Dynamic Mesh (DM) method [18,19]. On the other hand, Nguyen, et al. [20] utilized the MRF method to conduct virtual captive model tests for a full-scale submarine. Pan, et al. [21] used the DM method in numerical simulation to implement the submarine’s straight-line towing test and the pure heaving test. In addition, Wu, et al. [22] employed the Multi-block Hybrid Dynamic Grids (MHDG) method, including a multi-block mesh topology, moving mesh, and dynamic-layer method, to simulate AUV underwater docking with a discretized propeller. Generally speaking, the MRF method is adopted if the fluid domain moves based on the assumption of rigid motion. In contrast, the DM method is utilized if the fluid domain deforms as it moves. If unsteady interaction between the stationary and moving objects cannot be neglected in the multiple moving reference frames, the sliding mesh method is usually suggested to calculate the relative motions of the stationary and translational components [23]. In the case of steady-state maneuvering motions, only the velocity-based hydrodynamic coefficients can be determined generally. To determine the acceleration-based hydrodynamic coefficients, the Planar Motion Mechanism (PMM) experiments need to be performed [6]. Another approach to performing the PMM test is suggested by using the overlapping grid method in the CFD model [24]. However, the overlapping grid method is computationally expensive when the overlapped fractions have to be calculated afresh at each time step [25].

The purpose of this study is to explore the feasibility of the CFD model based on the software ANSYS Fluent 2019 R1 to evaluate the straight-line towing test by conducting the Boussinesq-based RANS solver with the SST k-ω turbulence closure and the sliding mesh method in the multi-grid fluid domain. The software SolidWorks 2017 was employed for the computer aided design (CAD) model. In order to lower mesh cell counts while maintaining solution accuracy, Fluent Meshing was adopted for mesh generation. The theoretical background and the numerical method would be introduced in the following sections. Subsequently, hydrodynamic characteristics of the full-appendage SUBOFF model were simulated and verified by the David Taylor Research Centre (DTRC) experimental data. Finally, the quantitative results were summarized to present the variations along the top meridian line of the full-appendage SUBOFF model in different flow fields.

2. Theoretical Background

2.1. Geometry

In this study, the full-appendage SUBOFF model (Configuration 8) [26] was considered in the numerical simulation. The dimension of the designed CAD model can be illustrated in Figure 1. The SUBOFF model is an axisymmetric hull with a total length of 4.356 m and a maximum diameter (D) of 0.508 m. The SUBOFF model has a sail above the hull, its leading edge is located at 0.924 m (1.820D) from the bow and the trailing edge is 1.293 m (2.545D) away, so the total length of the sail is 0.368 m (0.724D). There are four identical appendages placed on the stern, with angles of 0°, 90°, 180°, and 270° (vertical and horizontal control planes).

2.2. Coordinate System

Both the earth-fixed and the body-fixed coordinate systems in this study were defined by Feldman [27]. The velocity, angular velocity, corresponding force, and moment in each direction are as presented in Figure 2 and

Table 1. Parameters of the six degree of freedom (DOF) motion modes of the submerged body.

Mode	Force/Moment	Velocity/Angular Velocity	Displacement/Angle
Surge	X	u	x
Sway	Y	v	y
Heave	Z	w	z
Roll	K	p	ϕ
Pitch	M	q	θ
Yaw	N	r	ψ

, respectively. The hypothesis indicates that the fluid is incompressible, irrotational, and the free surface effect is ignored.

2.3. RANS Equations

In this study, the Reynolds numbers of the turbulent flows in the numerical simulation are all above ${10}^7$. RANS equations were adopted and can be written in Cartesian tensor form by substituting decomposed form of variables into continuity and momentum equations and taking time average as below:

$$\frac{\partial \rho }{\partial t}+\frac{\partial }{\partial x_i}\left(\rho u_i\right)=0,i, j =1,2,3$$

(1)

$$\frac{\partial }{\partial t}\left(\rho u_i\right)+\frac{\partial }{\partial x_j}\left(\rho u_i{u}_j\right)=-\frac{\partial p}{\partial x_i}+\frac{\partial }{\partial x_j}\left[\mu \left(\frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i}-\frac{2}{3}{\delta }_{ij}\frac{\partial u_i}{\partial x_i}\right)\right]+\frac{\partial }{\partial x_j}\left(-\rho \overline{u_i^{\prime }{u}_j^{\prime }}\right)$$

(2)

where $i$ and $j$ are tensor indices $\rho$ is the volume-fraction-averaged density; $u (=\overline{u}+u^{\prime }$) is the flow velocity and can be decomposed into the mean and fluctuating components; $\delta$ is the Kronecker delta; $p$ is the pressure term; and $\mu$ is the dynamic viscosity.

The additional terms in Equation (2) represent the effects of turbulence and are called Reynolds stresses. The Reynolds stresses $\left(-\rho \overline{u_i^{\prime }{u}_j^{\prime }}\right)$ based on the Boussinesq hypothesis [28] can be related to the mean velocity gradient, which is given as:

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

(3)

where ${\mu }_t$ denotes the turbulent viscosity, and $k$ represents the kinetic energy.

2.4. Turbulence Model

The Boussinesq hypothesis can be applied to the SST k-ω turbulence model [29], which is a two-equation eddy-viscosity model to improve the preciseness and computational stability of such calculations at near-wall regions. The advantages of the SST k-ω model are its capability of treating adverse pressure gradients and separating flow. It is able to produce turbulence levels which are a bit too large in regions with large normal strain, like stagnation regions and regions with strong acceleration. Meanwhile, the SST k-ω model is designed to be applied throughout the boundary layer provided that the grid resolutions in the near-wall are fine enough. Since model equations do not contain undefined terms at the wall, they are able to integrate to the wall without using wall functions. In that case, the accuracy increases accompanied with the decrease of the cell size or the non-dimensional wall distance (${\mathrm{y}}^{+}$) value [7]. On the other hand, the SST k-ω turbulence model used in this study belongs to the low Reynolds number turbulence model, which is suitable for processing in the near-wall treatment, the y^+ setting is thus chosen to approach 1.

More specifically, the SST k-ω model adopts a blending function to gradually transition from the standard k–ω model near the wall to a high Reynolds number k–ε model in the outer part of the boundary layer. It contains a modified turbulent viscosity formulation to calculate the transport effects of the principal turbulent shear stress. SST k-ω model can be expressed as follows:

$$\frac{\partial }{\partial t}\left(\rho k\right)+\frac{\partial }{\partial x_i}\left(\rho k{u}_i\right)=\frac{\partial }{\partial x_j}\left({\mathsf{\Gamma }}_k\frac{\partial k}{\partial x_j}\right)+G_k-Y_k+S_k$$

(4)

$$\frac{\partial }{\partial t}\left(\rho \omega \right)+\frac{\partial }{\partial x_i}\left(\rho \omega u_i\right)=\frac{\partial }{\partial x_j}\left({\mathsf{\Gamma }}_{\omega }\frac{\partial \omega }{\partial x_j}\right)+G_{\omega }-Y_{\omega }+D_{\omega }+S_k$$

(5)

$G_k$: the generation of turbulent kinetic energy $k$ due to the mean velocity gradient;

$G_{\omega }$: the generation of specific dissipation rate $\omega$ due to the mean velocity gradient;

${\mathsf{\Gamma }}_k$ and ${\mathsf{\Gamma }}_{\omega }$: the effective diffusion coefficients of $k$ and $\omega$;

$Y_k$ and $Y_{\omega }$: the dissipation of $k$ and $\omega$ due to turbulence;

$D_{\omega }$: the cross-diffusion term;

$S_k$ and $S_k$: undefined source terms.

2.5. Sliding Mesh Method

The sliding mesh method is one of several common approaches adopted to compute the unsteady flow field when a time-accurate solution is necessary for a moving object in the fluid. Specifically, the unsteady solution obtained from the sliding mesh method is in a time-periodic manner, which repeats with a period related to the speeds of the moving domains. For a moving grid with boundary motion, the integral form conservation equation can be described as:

$$\frac{d}{dt}{{\displaystyle \int }}_V\rho \varphi dV+{{\displaystyle \int }}_{\partial V}\rho \varphi \left(\overrightarrow{u}-{\overrightarrow{u}}_m\right)·d\overrightarrow{A}={{\displaystyle \int }}_{\partial V}\mathsf{\Gamma }\nabla \varphi ·d\overrightarrow{A}+{{\displaystyle \int }}_V{S}_{\varphi }dV$$

(6)

$$\frac{d}{dt}{{\displaystyle \int }}_V\rho \varphi dV=\frac{\left[{\left(\rho \varphi \right)}^{n+1}-{\left(\rho \varphi \right)}^n\right]V}{\Delta t}$$

(7)

where $\varphi$ is defined as an universal scalar; $V$ is an arbitrary control volume; $\stackrel{\rightharpoonup }{u}$ is the flow velocity vector; ${\overrightarrow{u}}_m$ is the moving mesh velocity; $\mathsf{\Gamma }$ is the diffusion coefficient; and $S_{\varphi }$ is the source term of $\varphi$. In addition, $n$ and $n+1$ represent the respective quantity at the current and the next time step.

The sliding mesh method is considered as a special case of a general dynamic mesh motion where all the meshes are allowed to move rigidly in a dynamic cell zone [30]. The advantage of the sliding mesh method is that neither the curved nor flat parts of the sliding interfaces between two cell zones introduce any local anomalies in the flow field, and flow variables vary smoothly across the mesh discontinuity [25]. Although not formally conservative, the sliding interface is effectively transparent and yields negligible mass inconsistency.

3. Numerical Solution

3.1. Solver

The finite volume method (FVM) employed in this study obtains a discrete equation based on a conservation equation of an integral form by dividing the computational area into fine grids and discretizing the flow field [31]. It is beneficial for the FVM method to allow for formulating unstructured meshes.

The present study utilized the semi-implicit method for pressure-linked equations-consistent (SIMPLEC) proposed by Patankar and Spalding [32]. SIMPLEC is an improved extension of the semi-implicit method for pressure-linked equations (SIMPLE) [33]. For complicated flows like turbulence, SIMPLEC is able to improve convergence only when it is limited by the pressure–velocity coupling.

Furthermore, the second-order upwind scheme [34] was adopted as the discretization scheme for convective terms of momentum equations and turbulence quantities transport equations, whereas the second-order central difference scheme was used for the diffusive fluxes. Therefore, it can be effective in estimating the flow field of a submerged body by adopting the numerical setup [35,36] mentioned above.

3.2. Fluid Domain

According to [21], the computational domain was designed to be a cylinder flow field. It is ensured that the length of the flow field was appropriate for capturing flow field characteristics, such as the following wake effect behind the submerged body at any speed. The length of the numerical flow field was set at 23.96 m (approximately 5.5 times the vehicle’s length, 5.5L). Moonesun, et al. [37] indicated that the width and height of the flow field can be set approximately seven times the vehicle’s diameter (7D) for the CFD modeling of the submerged body. In other words, provided that the height and width of the flow field were both set at 19 times the model’s diameter (19D), the simulation results of the numerical flow field would be further guaranteed to be free from the influence of the computational domain boundary. More details of the computational domain can be illustrated in Figure 3.

Figure 3 exhibited that the pressure outlets (B.C.1 and B.C.4) were set as the outlet boundary conditions at both ends of the fluid domain. Meanwhile, the free slip wall (B.C.2) was determined on the outside boundary of the external fluid region. For considering the effect of viscous flows on the rigid body, the boundary condition of the no-slip wall (B.C.3) was applied to the submerged body surface in the internal fluid region. Between the external and the internal regions, there is a non-conformal and sliding interface.

3.3. Mesh Generation

For the mesh generation of the SUBOFF model and the flow field, both the curvature size function and the proximity size function [38] were considered in this study. The present study adopted unstructured grids and structured grid for grid partition to properly fit the shape of the submerged body, for example, at the stern of the vehicle where the curvature varies greatly and the junction of the submerged body and appendages such as the rudder and sail, and use a structured grid in the flow field away from the vehicle to save computing resources. For grid refinement, common methods such as body sizing and face sizing can conduct additional refinement on areas with greater variations in model features or curvatures by appointing the mean length of the grid to a particular body, plane, or side. The results of multiple attempts revealed that body sizing tends to result in an excessive number of grids. Therefore, face sizing was adopted for grid refinement in the following simulation. To determine the flow field features near the wall, that is, the flow field at the boundary layer, the inflation function was applied to the CFD model.

According to [18], the multi-zone method can automatically decompose the geometric model into blocks, generate a structured hexahedral mesh in the area suitable for dividing the hexahedral mesh, and in other areas generate unstructured tetrahedral meshes. In Figure 4a, there are three regions in the computational domain [39]: one is the interior region (region B) that surrounds the unstructured tetrahedral meshes of the object (region C) without relative motion with the object, the other one is the exterior region (region A) that excludes the computational domain around SUBOFF. More details of mesh properties can be introduced in

Table 2. Mesh properties in the multi-zone fluid domain.

Zone	Sizing	Element Size (m)	Type
A	body sizing	0.112	hexahedral
B	body sizing	0.065	tetrahedral
C	face sizing	0.037	tetrahedral

. Furthermore, Figure 4b shows the interface between region A and region B. In the schematic diagram, two regions of cells meet at the non-conformal interface. Specifically, fluxes across the grid interface would be calculated using the faces resulting from the intersection of the two interface zones, e.g., faces 9-2 and 2-10, instead of the interface zone faces, e.g., face IX-X themselves.

3.4. Grid Convergence Index

Numerical uncertainty is often induced by the errors in the CFD model and comprises the effects of discretization, round-off, and iterative convergence errors on the accuracy [40]. The grid refinement is considered as the most common solution to numerical uncertainties by using grid convergence index (GCI) [41,42]. Basically, the GCI relates to the comparison of discrete solutions at different grid resolutions and gives an indication of how much further grid refinement would change the numerical solution [43]. The GCI estimated for the grid independence test of the resistance over two grid solutions can be expressed as:

$$GCI=\frac{f_s}{s_R^p}\left|\frac{R_{i+1}-R_i}{R_i}\right|·100\% ,i=1, 2, \dots , N$$

(8)

where $f_s$ is defined as a safety factor and implies 95% confidence for the uncertainty estimate [42], $R$ is the resistance, $f_s=1.25$ is suggested for employment over these grids [44]; $N$ is the total number of grid sizes for GCI tests, $s_R$ is the grid refinement ratio, and $p$ is the observed order of convergence which can be formulated as:

$$p=\frac{ln\left|\left(R_{i+2}-R_{i+1}\right)/\left(R_{i+1}-R_i\right)\right|}{ln\left|s_R\right|},i=1, 2, \dots , N-2$$

(9)

3.5. Courant–Friedrichs–Lewy Condition

The Courant–Friedrichs–Lewy condition is principally used for numerical convergence. More specifically, temporal discretization of governing equations is accomplished by an explicit time-marching algorithm [45]. Therefore, the determination of computational stability and convergence is achieved by revealing the change of the dimensionless parameter, the Courant number, to understand the correlation between grid size and time step [46] as below:

$$C_{courant}=\frac{u·\Delta t}{\Delta x}$$

(10)

where $u$ is denoted as the towing velocity; $\Delta t$ is the time step; and $\Delta x$ is the grid size.

4. Results and Discussion

4.1. Simulation Condition

In order to ensure the accuracy of computational results, the fully appended DARPA SUBOFF Model 5470 (Configuration 8) stated in the resistance experiment data of the DTRC was used as a test object. Relevant experimental conditions are as presented in

Table 3. Experimental conditions for resistance tests.

Speed (knot)	Speed (m/s)	Reynolds Number	Resistance of DTRC Experimental Data (N)
5.93	3.050	1.33$\times$ 107	102.3
10	5.144	2.23$\times$107	283.8
11.85	6.096	2.65 $\times$107	389.2
13.92	7.160	3.11$\times$ 107	526.6
16	8.230	3.58 $\times$107	675.6
17.79	9.151	3.98 $\times$ 107	821.1
Model Configuration No.8 (Fully Appended)

[47].

4.2. Grid Independence Analysis

Before the grid independence analysis, the iterative convergence criteria adopted for all numerical simulations is a maximum normalized residual ($L_{\infty }$) [48] of the order of 10⁻⁶. For a given Courant number, viz. the grid size is fixed, the time step decreases in accordance with the increase of the towing velocity which is related to the Reynolds number. In order to extract stable data for further processing, the mean values of time histories between 10 and 30 s would be presented.

During the test, the mean length of the grid cells around the submerged body was gradually reduced by adjusting the face sizing and body sizing setting to increase the number of grids. In order to investigate the effect of grid size on the resistance computation, only the base cell size was changed for obtaining geometrically similar grids. In addition, the first cell height was prescribed the same value for all grids when the near-wall treatment was considered. This grid refinement method [38] can place fine grids on the surface of the submerged body. Similarly, fine grids were placed on the sail, stern, and other appendages where the curvature varied remarkably. Then, by adjusting the growth rate of the grid, their skewness did not increase enough to affect the simulation results. With respect to the resistance test, the forces and moments in all directions except for the components in the longitudinal direction were small. Therefore, only the GCI values analyzed in the longitudinal direction as a function of the Reynolds number were compared. Specifically speaking, the GCI curves not only present the convergence of grid size, but also the influence of time step.

In Figure 5, GCI₁₂ represents the analysis results of $C_{courant}$ = 2.66 and 2.22, and GCI₂₃ represents the analysis results of $C_{courant}$ = 2.22 and 1.85. The corresponding Courant–Friedrichs–Lewy conditions were presented in

Table 4. The respective Courant–Friedrichs–Lewy conditions for GCI analysis.

${\mathit{C}}_{\mathit{c}\mathit{o}\mathit{u}\mathit{r}\mathit{a}\mathit{n}\mathit{t}}$	Reynolds Number ( $\mathit{R}\mathit{e}$)	$\mathbf{\Delta }\mathit{x}$ (m)	$\mathbf{\Delta }\mathit{t}$ (s)
1.85	1.33$\times$107	$1.65\times {10}^{-2}$	$1\times {10}^{-2}$
	2.23$\times$ 107		$5.93\times {10}^{-3}$
	2.65 $\times$107		$5\times {10}^{-3}$
	3.11$\times$ 107		$4.26\times {10}^{-3}$
	3.58 $\times$107		$3.7\times {10}^{-3}$
	3.98$\times$ 107		$3.33\times {10}^{-3}$
2.22	1.33$\times$107	$1.37\times {10}^{-2}$	$1\times {10}^{-2}$
	2.23$\times$ 107		$5.93\times {10}^{-3}$
	2.65 $\times$107		$5\times {10}^{-3}$
	3.11$\times$ 107		$4.26\times {10}^{-3}$
	3.58 $\times$107		$3.7\times {10}^{-3}$
	3.98$\times$ 107		$3.33\times {10}^{-3}$
2.66	1.33 $\times$107	$1.15\times {10}^{-2}$	$1\times {10}^{-2}$
	2.23$\times$ 107		$5.93\times {10}^{-3}$
	2.65 $\times$107		$5\times {10}^{-3}$
	3.11$\times$ 107		$4.26\times {10}^{-3}$
	3.58 $\times$107		$3.7\times {10}^{-3}$
	3.98$\times$ 107		$3.33\times {10}^{-3}$

. It is obvious that the GCI values decrease accompanied by successive grid refinement for a given $Re$ number. On the other hand, the variation of the time step appears to have little influence on the GCI curve with the increase of the Reynolds number. However, the results at high Reynolds numbers are slightly inconsistent with the ones at low Reynolds numbers. It is due to the fact that the SST k-ω turbulence model may produce a bit too large turbulence levels in regions with large normal strain, like stagnation regions and regions with strong acceleration, especially at a high Reynolds number. Since all GCI values are below 0.07%, it is verified that the dependency of numerical results on the grid size has been reduced and the solution achieves the grid-independent solution.

4.3. Resistance Analysis

Figure 6 shows the simulated results of the resistance from $C_{courant}=1.85, 2.22, 2.66$ in comparison with the DTRC experimental data and the published numerical data [17,49]. According to the GCI analysis, the grid-independent solution can be essentially achieved among these three Courant numbers. Besides, the simulation results appear to agree generally well with the DTRC experimental data and the other two published numerical data. It is shown that the divergence between the numerical simulation ($C_{courant}=2.22$ and $2.66$) and the experiment would gradually become large with the increase of Reynolds number. Nevertheless, the simulation results of $C_{courant}=1.85$ are still found to be quite close to the DTRC experimental data and even have better predictions than the other two published numerical results. Özden et al. [17] established a RANS solver coupled with the SST k-ω turbulence model to implement the resistance simulation of the DARPA SUBOFF model (Configuration 8) in the incompressible flows based on the SIMPLE algorithm. In addition, Sezen et al. [49] used a CFD model, in which the RANS equations were coupled with the k–ε turbulence model and the pressure field was solved based on the SIMPLE algorithm, to simulate the resistance analysis of the DARPA SUBOFF model with full appendages. For both these cases, the moving reference frame (MRF) method was only adopted for open water and self-propulsion analysis. In comparison, the present numerical simulation utilized the sliding mesh model to simulate the straight-line towing test. Accordingly, it seems to be appropriate for comparing our numerical simulation with the resistance experiment by setting $C_{courant}=1.85$ and the sliding mesh model in the following hydrodynamic analysis.

4.4. Analysis of Turbulent Flow-Field Characteristics

Figure 7a–f present the vorticity distribution around the SUBOFF hull surface for different towing speeds. It can be seen from Figure 7a–f that as the speed increases, the intensities of vortices also increase significantly. It is exhibited that the flow pattern becomes complicated as the flow with embedded vortices approaching the SUBOFF stern. According to Chase and Carrica [50], the flow field is encountering an adverse pressure gradient when the hull geometry is sharpened. As a result, the boundary layer will thicken rapidly accompanied by the decrease of the axial flow velocity. Moreover, the vortices separated from the sail are transferred to the stern and amplified when the boundary layer is evidently thickened. Since all these vortex structures interact with each other, the non-uniformity of the inflow field behind the rudders would be aggravated significantly for the large towing speed.

Figure 8a–f present the normalized velocity contours around the SUBOFF model under different towing speeds. In order to highlight the flow variation around the SUBOFF model, the velocity contours have been normalized as U/U₀, where U is the flow velocity relative to the towing velocity U₀. It is found that the boundary layer develops from the stagnation point. Subsequently, the boundary layer and the upstream flow interact with the leading edge of the sail, which causes a local velocity gradient in front of the sail. The adverse velocity gradient occurs near the stern due to the contraction of the hull shape and contributes to the flow separation of the boundary layer. When the boundary layer flow from the hull interacts with the rudders, it causes rapid rises of local velocities in front of the rudders. Eventually, the shear layers shed from the hull and the rudders are propagated downstream into the wakes. Moreover, the extent of radiation to the surroundings also increased as the towing speed increased.

Figure 9 shows the friction coefficients $C_f$ along the top meridian line by comparing simulation results for different Reynolds numbers Re with the experimental data from Huang, et al. [51], and Qiu, et al. [52] for $Re=1.2\times {10}^7$. The friction coefficients $C_f$ is defined as $C_f={\tau }_w/\frac{1}{2}\rho U_0^2$, where ${\tau }_w$ is frictional stress of the hull surface, $\rho$ is the mass density of the fluid, and $U_0$ is the given towing velocity. By discussing the influence of the Reynolds number $Re$ on the friction coefficient $C_f$, the results indicate that $C_f$ decreases with the increase of $Re$. Due to the interaction of boundary layer flows with leading edges of appendages, the local peaks of $C_f$ occur in front of the sail and the rudders. The decrease of $C_f$ from the top of the sail is accompanied by the generation of the tip flow, which moves downstream in the form of the tip vortex. In addition, boundary separation interacting with shear layers from the hull and the rudders may cause a decrease of $C_f$.

Figure 10a–f present the normalized velocity distribution on the cross-section of the sail at x/L = 0.255 of the SUBOFF model under different towing speeds. It reveals that adding a sail to the hull notably affects the resistance and flow field of the submerged body [9]. The thickening of the hull boundary layer is due to the adverse pressure gradient on the sail, which eventually leads to flow separation and wake formation. The stern wake produced by the sail affected the overall wake distribution. Drastic changes in flow velocity generated a horseshoe vortex at the front edge of the sail, which produced additional fluid-borne noise and wake. Therefore, the shape and local flow field of a sail must be carefully examined, and the interaction between the sail and the hull must be considered.

A detailed observation of the normalized velocity contours around the rudders at x/L = 0.896 of the SUBOFF model under different speeds has been illustrated in Figure 11a–f. Due to the boundary layer interacting with the rudders, the hole in the center of each wake profile is generated. It is found that the wake pattern for each case is very similar, with a reduced axial velocity around the rudders. Such non-uniform flow may cause notable noise when it meets the propeller. Therefore, propeller designers must consider the wake distribution at the stern rather than evaluating its performance solely on the basis of propeller properties.

In addition to velocity distribution in the flow field, the pressure distribution on the submerged body is a substantial factor in the design process. Although the change of pressure caused by flow velocity change cannot be compared to the pressure generated by seawater as the submerged body dives deeper underwater, the distribution of high and low pressure on the submerged body serves as an additional reference for evaluating the shape design. Figure 12a–f present the distribution of pressure coefficients $C_P$ on the SUBOFF model under different speeds. The pressure coefficient is defined as $C_P=\left(P-P_0\right)/\frac{1}{2}\rho U_0^2$, where P is dynamic pressure and $P_0$ is static pressure.Figure 12a–f reveal that pressure distribution around the submerged body did not change notably with the towing speed, and the areas with the highest pressure were the bow, front edge of the sail, and front edge of the stern appendage. Moreover, low-pressure areas were observed where the flow accelerated, such as the bow, sail, and stern, whereas no notable pressure changes were observed amidships.

Furthermore, Figure 13 exhibits comparisons of simulated and experimental pressure coefficients $C_P$ along the top meridian line of the SUBOFF model with full appendages under different $Re$ conditions. From the comparison of SUBOFF without appendages and full appendages in Figure 13, it can be seen that the $C_p$ value is related to the geometric shape, and the wavering of the $C_p$ value is most obvious before the position of the sail (x/L = 0.2). In front of the after body (x/L = 0.8), the $C_p$ value rises relatively as the hull diameter decreases and also decreases in front of the rudder. The experimental data from [51] were measured at $Re=1.2\times {10}^7$. It is obviously found that the Reynolds number $Re$ has almost a tiny effect on the variation of the pressure coefficient $C_P$. In general, the simulation results are almost consistent with the experimental data from [51]. Meanwhile, the simulation results reveal that $C_P$ has a maximum value at the stagnation point (x/L = 0), and then decreases rapidly before it reaches the trip wire (x/L = 0.04). In the rear of the tripwire, $C_P$ has a local peak at the top meridian in front of the sail (x/L = 0.21) and then has a sharp decrease. Since the wake from sail affects the pressure before x/L = 0.4, $C_P$ varies constantly only in the region 0.4 < x/L < 0.7. Apparently, the adverse pressure gradient occurs in the stern (0.7 < x/L < 0.88) and then reaches its local maximum in front of the rudders (x/L = 0.88). Definitely, the variation of $C_P$ in the stern (0.7 < x/L < 0.88) is caused by the interaction of boundary layers with the rudders.

5. Conclusions

In this study, the sliding mesh model in the multi-zone fluid domain was integrated with the Boussinesq-based URANS solver and the SST k-ω turbulence closure into the CFD model. After the grid-independent tests, the straight-line towing test of a full-appendage DARPA SUBOFF model was conducted by analyzing the resistances, turbulent flow fields, vorticity, frictional coefficients, and pressure coefficients. Since the simulation results of $C_{courant}=1.85$ are mostly in good agreement with the DTRC experimental data, the hydrodynamic characteristics have been analyzed qualitatively and quantitatively. When the multi-zone method was combined with the sliding mesh model, it ensures the good accuracy of the unsteady motion. It is helpful to enhance the computation efficiency and reduce the numerical error due to the difference in the calculation domain. By discussing different towing speeds, the results indicate that $C_f$ decreases with the increase of $Re$. It is also shown that the local peaks of $C_f$ occur in front of the sail and the rudders as a result of the interaction of boundary layer flows with leading edges of appendages. In contrast, it appears that there is less variation of the pressure coefficient $C_P$ with increasing $Re$.