3.1 Potential Flow Theory
In the present study, CFD code and potential flow theory were used for seakeeping analysis of the submarine, and the results were compared to those of the model tests for verification. The seakeeping analysis based on potential flow theory was performed with Hydrostar for Experts V8.10, which has been developed by Bureau Veritas since 1991. Hydrostar is powerful 3D diffraction/radiation potential software that provides a complete solution to a first-order problem of wave diffraction and radiation and also second-order low-frequency wave loads for a floating body with or without forward speed in deep water and in finite water depth (
BV, 2019). For 3D diffraction/radiation potential theory, the hull was modeled as panels according the hull geometry. The number of panels is governed by the panel size, and the maximal panel length should be smaller than 1/6 of the minimum wave length considered for analysis. The panel model used in Hydrostar is shown in
Fig. 2.
The seakeeping analysis was performed at all wave headings from head (180º) to following (0º) seas (12 headings with 30º spacing), and a forward speed of 3 knots (1.54 m/s) was used. Wave frequencies of 0.27 rad/s to 1.32 rad/s were used with steps in wave frequencies not exceeding 0.04 rad/s. Since the viscosity effect on the roll response is important to arrive at a realistic roll angle, appropriate roll viscous damping should be included. In this study, 5% viscous roll damping was taken from a roll decay test in a published paper (
Hermanski and Kim, 2010).
After the modeling and mesh generation were completed, the response amplitude operators (RAOs) with six degrees of freedom were calculated, and the results of heave, roll, and pitch are shown in
Figs. 3,
–
5. From the results, unlike the heave and pitch RAOs of typical surface ship, those of a surfaced submarine take double peaks into account. Roll RAOs of quartering waves (30º, 60º and 120º) are higher than that of beam sea (90º). These results seem to be due to the nonlinearity near the free surface area, as will be shown later.
3.2 Computational Fluid Dynamics
The seakeeping analysis based on the regular wave simulation was performed with the commercial CFD software STAR-CCM+ V15.06.
Fig. 6 shows the computational domain that was used in this study. The computational domain consists of two regions (a tank region and overset region) and was generated with size corresponding to 0.7 times the ship length (
L) from the forward direction of the ship to 3.3 times L from the backward direction, 3.0 times Lpp in the horizontal direction, and 2.5 times Lpp in the vertical direction. The dynamic fluid body interaction method with overset grid and two sets of spring damper system were used to simulate the 6 degrees of freedom of motion.
For body-environment coupling, a pair of linear springs was attached to the bow and stern, and an elasticity coefficient of 100 N/m was used to minimize the effect on the ship’s motions. The boundary conditions for the inlet, outlet, top, bottom and side (port and starboard) were assigned as velocity inlets, and the velocity was defined as the same flow as the submarine’s target speed. For considering the free surface, the volume of fluid (VOF) method was used, and a VOF wave forcing zone was applied to the inlet (forcing zone size: 0.3L), outlet (forcing zone size: 1.5L), and side (forcing zone size: 0.25L) to eliminate disturbed flows.
The simulation conditions for regular wave cases are described in
Table 2. The regular wave simulations were performed with seven wave headings and forward speeds of 3 knots (1.54 m/s). The wavelength (
λ) was between 0.5 and 3.5 times the ship length, and the wave steepness (
H/λ) was fixed with a value of 1/50. The turbulence was modeled with the realizable
k-ε turbulence model for numerical stability and is generally used for the marine industry. In order to reduce the wave dissipation, a second-order scheme in time was applied, and a 5th-order Stokes wave was used in the simulation. The time step was selected to be able to perform 500 calculations per period based on the encounter period, and the calculation was performed so that the encounter period was repeated 10 time or more.
The mesh was made with a trimmed mesh in STAR-CCM+. Four prism layers were used to accurately measure the boundary layer flow on the hull surface, and the first layer thickness (Y+) was 10–40. The configuration of the mesh used in a regular wave simulation is shown in
Fig. 7 and
Fig. 8.
A grid convergence test was performed to find the optimal grid resolution for three grids based on Richardson extrapolation (
ITTC Resistance Committee, 2017). For this purpose, it was performed on course, medium, and fine grids for heave, pitch, and roll RAOs. The three grids are determined by increasing the size of the grid with the refinement ratio (
ri ). Based on the simulation time, a refinement ratio (
ri ) of
2 was applied due to the hardware resource limitation. The number of grids corresponding to the course, medium, and fine grids are shown in
Table 3, and the generated grids are shown in
Fig. 9.
The heave, pitch, and roll RAOs are selected for grid convergence test: ship on regular waves, head sea (
μ = 180º) and beam sea (
μ = 90º). The both cases ship speed was 3 knots (1.54 m/s) and the results of grid convergence study are summarized in
Table 4. The convergence ratio is defined by
Eq. (1).
where
Convergence conditions are defined according to the convergence rate calculated by Richardson extrapolation as follows:
Ri > 1: Grid divergence
Ri < 0: Oscillatory convergence
0 < Ri < 1: Monotonic convergence
As shown in
Table 4, RG for heave, pitch, and roll shows that monotonic convergence is achieved for all simulation variables because they are between 0 and 1. It is shown that the tendency of calculated results converge as the size of the grid is reduced. The order of accuracy
pG is defined by
Eq. (2).
Based on
Celik et al. (2008), the Grid Convergence Index (GCI) is defined using the following equations:
As shown in
Table 4, the maximum GCI is 0.029 in heave motion and head sea, and the minimum GCI is 0.008 in roll motion and beam sea, all of which are within an acceptable range. In this study, grids of 4.9 million cells were used based on the grid convergence test, and 100 grids per wavelength and 16 grids per wave height were selected. From the CFD simulation, heave, roll, and pitch responses were obtained based on the earth-fixed coordinate system.
Fig. 10 shows the time histories of the roll and pitch response. The motion RAOs were calculated from average of responses after the transient range, as shown in
Fig. 10.
Figs. 11,
–
13 show the motion RAOs obtained from the CFD simulation. The motion RAOs were normalized with the wave amplitude. Unlike the results of potential code, the maximum roll response was found in beam sea. The natural roll period was 8.3 seconds, which is a slightly smaller value compared with the target value (8.76 s). In the motion RAOs, the maximum values of heave and pitch were found in beam sea and head sea, respectively. As with the potential code, the heave and pitch RAOs from CFD simulation showed double peaks, which will be discussed later with the comparison results.