CFD Research on the Hydrodynamic Performance of Submarine Sailing near the Free Surface with Long-Crested Waves

: The simulations of submarine sailing near the free surface with long-crested waves have been conducted in this study using an in-house viscous URANS solver with an overset grid approach. First, the veriﬁcation and validation procedures were performed to evaluate the reliability, with the results showing that the generation of irregular waves is adequately accurate and the results of total resistance are in good agreement with EFD. Next, three different submerged depths ranging from 1.1D to 3.3D were selected and the corresponding conditions of submarine sailing near calm water were simulated, the results of which were then compared with each other to investigate the inﬂuence of irregular waves and submerged depths. The simulations of the model near calm water at different submerged depths demonstrated that the free surface will cause increasing resistance, lift, and bow-up moments of the model, and this inﬂuence decreases dramatically with greater submerged depths. The results of the irregular wave simulations showed that irregular waves cause considerable ﬂuctuations of hydrodynamic force and moments, and that this inﬂuence remains even at a deeper submerged depth, which can complicate the control strategies of the submarine. The response spectrum of hydrodynamic forces and moments showed slight amplitudes in the high-frequency region, and the model showed less sensitivity to high-frequency excitations. studied with different turning radii and drift angles. The speciﬁc results for the total force and moments, as well as lateral force distributions along the hull and the ﬂow properties such as ﬂow separation and leeside vortex development, are presented and compared with experimental data. Dubbioso and colleagues carried out simulations of a submarine with two different aft appendage con-ﬁgurations and analyzed the maneuvering response, including trajectories and kinematic parameters. The results showed that the X rudder conﬁguration performs better regarding turning qualities with respect to the cruciform rudder. All of the simulations conﬁrmed that CFD can be used to simulate the steady turning of submarine.


Introduction
It is necessary to evaluate the sailing performance parameters of a submarine to ensure its safe navigation before sea trials.There are many different methods to obtain the various performance parameters, such as experimental fluid dynamics (EFD), computational fluid dynamics (CFD), and extrapolation according to the empirical formulation.Previous research [1,2] conducted experiments and collected serial experimental data for SUBOFF and established a SUBOFF database, which has provided data resources for CFD validation and other submarine research.Nathan Chase [3,4] conducted simulations of self-propulsion for SUBOFF with the E1619 propeller, showing that CFD can be used to simulate the self-propulsion of submarines.Further studies [5,6] have simulated the steady turning of a submarine, with several submarine models being studied with different turning radii and drift angles.The specific results for the total force and moments, as well as lateral force distributions along the hull and the flow properties such as flow separation and leeside vortex development, are presented and compared with experimental data.Dubbioso and colleagues [6] carried out simulations of a submarine with two different aft appendage configurations and analyzed the maneuvering response, including trajectories and kinematic parameters.The results showed that the X rudder configuration performs better regarding turning qualities with respect to the cruciform rudder.All of the simulations confirmed that CFD can be used to simulate the steady turning of submarine.
In summary, the hydrodynamic and sailing performance of submarines have been studied in detail in the open water.Nevertheless, submarines will sail close to or beyond the free surface during their voyage due to some unavoidable needs, including replenishment, communication, reconnaissance, and emergent buoyantly rising.When a submarine sails close to the free surface, surface waves will be generated due to an interaction between the hull and the free surface.The dynamic pressure distribution along the hull will be affected by the surface waves further.As a result, the hydrodynamic performance of the submarine will be changed and perform differently with respect to its performance in open water.The mechanisms of interaction between the free surface and the underwater vehicle have been studied by Dawson [7], Amiri et al. [8], and Li et al. [9].The conclusions of above literature can be briefly summarized as follows.The bow, aft shoulder, and stern wave systems generated by the submerged vehicle will interact with each other, which forms the whole wave pattern.With different sailing velocities, the phase difference between several wave systems will cause an increase or decrease of the surface wave elevations and change the wave lengths.The hydrodynamic coefficients will be affected due to the change of wave elevations and wave lengths, and the effect of the free surface will decrease quickly as the submerged depth increases.In conclusion, with the submarine sailing close to the free surface, the resistance will increase and the lift and pitch moments show different features with respect to that of open water.To evaluate the influence of the free surface acting on the hydrodynamic performance of submerged vehicles, Wilson-Haffenden et al. [10] and Dawson [7] have conducted several experiments to investigate the influence of the free surface on the resistance, lift, and moments of submerged vehicles, which revealed that the free surface affects the wave-making resistance and dynamic pressure distribution along the hull, thus changing the lift and pitch moment of the vehicles.Jagadeesh et al. [11], Mansoorzadeh et al. [12], and Nematollahi et al. [13] all simulated the underwater vehicle sailing near the free surface by solving the RANS equation.They compared the simulation results with experimental data, indicating that the RANS method can be used to obtain reliable hydrodynamic coefficients.The resistance coefficients increase with decreasing submerged depths, which is mainly owing to the increase of the dynamic pressure component.In addition to the hydrodynamic performance, the sailing performance of submarines has also been discussed widely.Attiya et al. [14] studied the lift and drag performance of a finite plate near the free surface using the LES method, in which the effects of aspect ratio were investigated and the tendencies were explained by visual vortex structures.The mean drag coefficient decreased slightly with an increased aspect ratio, which was owing to the vorticial activities near the front face of the plate.Carrica et al. [15] and Wang et al. [16] conducted experiments and simulations for the self-propulsion of a submarine near the free surface.The former paid attention to the force and moment coefficients of the submarine and performance parameters of the propeller in different control strategies.the latter focused on the wake and vortex of the propeller, and compared the CFD results with EFD, thereby evaluating the reliability of CFD and pointing out the distortions regarding wake fluctuation and tip vortex strength.The research presented some consistent conclusions that the interaction between the free surface and the submarine causes the non-uniform submarine wake, affecting the performance parameters of the propeller, causing fluctuations in thrust and torque coefficients, and complicating the control systems of the underwater vehicle.
With the increase of studies focused on the interaction between a submerged vehicle and the free surface, some researchers have explored the influence of initial surface waves acting on the submarine sailing beneath the free surface.Willy [17] predicted that the force and moments acting on the submerged AUV were due to the waves based on the slender body theory, and the corresponding experiment was conducted to be compared with the simulation results.This method was found to be in good agreement with the experimental data, and the hydrodynamic forces on the model were linearly related to the wave amplitude.Milgram [18] used the strip theory to simulate the hydrodynamic characteristics of a submerged vehicle, considering the wave-induced forces and radiation components due to the vehicle motions.The results agreed well with the experiments, except for the simulations in beam seas.The fin-induced forces were determined to be important for the prediction of wave forces and radiation components.Tian and colleagues [19] numerically simulated a bare AUV sailing near the free surface with regular waves and the 6DOF motions of the model were fixed.The study concluded that the hydrodynamic coefficients show regular periodic fluctuations, and the amplitude of the fluctuations increases with an increased wave height or a decreased submerged depth and Reynolds number.The fluctuation of the hydrodynamic coefficients was explained with the periodic position of the wave trough and crest changes.Kim et al. [20] simulated the steady turning of a Joubert BB2 model near the free surface with irregular waves.They freed 6DOF motions of the submarine and explored the effects of irregular waves on the turning trajectories and kinematic parameters.The results showed that the turning maneuver is complicated due to the changes in broadband encounter frequency with irregular waves.The control performance was degraded with the lower velocity of the submarine, and the heading angles between the submarine and the wave direction had considerable influences on submarine motion.
Although there has been some previous research related to the effects of surface waves acting on the underwater vehicle, the studies were partly limited by the inherent assumption, i.e., the simplified sea state with regular waves and the theory method with the inviscid approximation.The present study focuses on the hydrodynamic performance of a submerged vehicle based on the viscous flow method, considering the effects of the initial surface wave, and the ITTC wave spectrum is discretized numerically to simulate the irregular waves, which could be used to provide support for the development of control strategies for a vehicle sailing near the free surface with realistic seaway conditions.Herein, a captive SUBOFF model sailing near the free surface with calm water and longcrested waves in different submerged depths has been simulated with an in-house viscous URANS CFD solver.An introduction to the geometry of the model, as well as simulation conditions, is presented in Section 2. Section 3 introduces the numerical methods applied in this study, including the discretization of the governing equation, the overset grid approach, and the wave-generation method.The computational domain and grid generation are described in Section 4. The verification for the numerical method has been conducted through grid/time-step convergence studies, with the numerical wave-generation technique and the results of the total resistance being validated by comparing with theoretical solutions or available experimental data in Section 5.The hydrodynamic parameters of the submarine and flow characteristics in different simulation conditions were compared to evaluate the influence of surface waves and submerged depths in Section 6. Conclusions are presented in Section 7.

Geometry of the Model
A DARPA SUBOFF model (AFF-2) was selected to be investigated in the present study.Compared with the full appendages model (AFF-8), the four identical stern appendages were simplified and only the bridge fairwater was retained, which was aimed to compare with the available experimental data.The geometry of the submarine is shown in Figure 1 and the parameters are shown in Table 1.

Simulation Conditions
The general view of the simulation conditions is represented in Figure 2, where a captive model is set below the free surface, and h refers to the submerged depth, which is measured from the axis of the main hull to the free surface.The bow and stern of the model are located at x/L = 0 and x/L = 1 in x direction, respectively.The free surface is set at z = 0 in z direction, which is the origin of the amplitude fluctuation for irregular waves.The inlet flow velocity is set according to the Froude number Fn based on hull length, where the Fn is defined as follows: where v refers to ship speed and g refers to acceleration of gravity.
To evaluate the effects of irregular waves, the cases where the model sails near the free surface with calm water were investigated first, where simulations were conducted at different submerged depths.Additionally, models sailing near the free surface with irregular waves in corresponding submerged depths were simulated, which was carried out to last 100 mean periods and total time steps were equal to N = 25,000.The simulation conditions are described in Table 2, where the flow velocity is constant for all cases, and a sea state of 4 with a significant wave height of 2.5 m and mean period of 7.2 s are selected to be the initial surface waves, which has been discretized numerically according to the ITTC wave spectrum and scaled according to the scale ratio λ = 24 and Fn = 0.336.
The present simulations were conducted based on the Tianhe-2 supercomputer at the National Supercomputer Center in Guangzhou with a parallel computational technique.The supercomputer includes 17,920 processors (Intel ® Xeon ® E5-2692 v2 Processor) and 12 cores are included for each one.The simulations used 56 cores for the simulations with calm water, and each simulation takes a physical time of about 52 h to complete.In the meantime, 155 cores were used for the simulations with long-crested waves, and each simulation takes a physical time of about 214 h to complete.

Numerical Method
The unsteady RANS equation is solved by the viscous URANS solver.The continuity and momentum conservation equations are defined as: ∂ ∂ where ρ refers to the density; p refers to the time-averaged pressure; U i (i = 1, 2, 3) refers to the time-averaged velocities component and u i (i = 1, 2, 3) refers to the fluctuation velocities component; μ refers to the dynamic viscosity; and ρu i u j refers to the Reynolds stress term.
To solve the Reynolds stress term, a two-equation eddy-viscosity SST k-ω turbulence model proposed by Menter [21] has been applied to close the RANS equation.The model is effective in predicting an adverse pressure gradient and shows low sensitivity to the inlet conditions.For solving the governing equation, the finite difference method is used to discretize the partial differential equation with temporal and spatial term discretization using the second-order differencing schemes.Meanwhile, the projection algorithm [22,23] is selected to solve the pressure-velocity coupling equation and the free surface is captured with the level-set method [24].To simplify the generation of structured grids for the complicated geometry of the ship and maintain a high quality of grid cells, an overset grid approach has been applied to generate the whole computation domain grid, using an inhouse overset code.The overall grid can be divided into different parts, such as background grids, hull grids, and appendage grids.Those grids overlap with each other and the node information is not corresponding among each other before the overset.In the process of the overset, the hole-cutting operations are conducted to dig the nodes which traverse the surface of solid boundaries and are beyond the computational domain using the cut-and-paste algorithm [25,26].Those nodes are defined as hole nodes which are not participating in the fluid computations.The grid nodes near the interface between different grid parts will be marked as fringe nodes, which are used to transmit the computing information between different grid parts.The computing information originates from nearby active nodes, which are titled as donor nodes, and the computing information among the different grid parts is connected by building the interpolation relationships between the fringe nodes and the donor nodes.More details can be found in the literature [27,28].

Irregular Wave Generation
The irregular waves are formed by various linear regular waves with the same direction, but with different wave lengths, wave amplitudes, and random phases for the longcrested waves.The wave elevation, η, is defined as the function of location x and time t: where N refers to the number of linear regular waves.For the ith regular wave, A i refers to wave amplitude; k i refers to wave number; ω i refers to the circular frequency; and φ i refers to the random phase from 0 to 2π.A specified wave energy spectrum was selected to determine the parameters of different linear regular waves.The spectral density function, S(ω), is defined as: where A i refers to the wave amplitude of the ith regular wave and ∆ω i refers to the wave frequency interval.In this study, the ITTC standard wave spectrum [29] was selected to simulate the irregular waves and the formulation is defined as: where A and B are constants ), H refers to significant wave height, and T refers to the mean period.
The wave energy spectrum is discretized by the equal-frequency method where different regular wave components have the same frequency interval.The discretization results will determine the parameters of the regular wave components.

Damping Method
To avoid the reflection from the exit boundary, a damping method was used for the wave generation.The expression of the momentum equation after damping modification is defined as: where μ(x) is the damping coefficient varying with the x-position as follows: where x > x 0 is defined as the damping region; L s refers to the length of damping region; and as refers to the wave-absorbing coefficient.

Computational Domain and Boundary Conditions
For different simulations with calm water and irregular waves, the ranges of the domain are set to be consistent in the x and y directions, and the domain extends from −1 < x/L < 3.5, −1 < y/L < 1. Herein, the size of the domain is different for two type cases in the z direction, in which −1 < z/L < 0.4 for the calm water and −1.5 < z/L < 0.4 for the irregular waves.The additional distance for the irregular waves aims to ensure the deep-water condition.The computational domain is shown in Figure 3.
The boundary conditions are set as follows: the velocity inlet and pressure outlet boundary conditions are applied on the front and end of the computational domain, respectively; the zero-pressure and zero-pressure-gradient boundary conditions are applied on the bottom and top of the computational domain, respectively; two sides are set as the zero-pressure-gradient boundary condition; and the surfaces of the model are set as a fixed no-slip boundary condition.

Grid Generation
The fully structured overset grid was applied in the present study, which could be divided into three parts before the overset: background grids, hull grids, and fairwater grids.The background grids are generated as Cartesian-type topology, and the hull and fairwater grids are set as O-type topology.For different simulations with calm water and irregular waves, the sizes of the hull grids and appendage grids are selected to be consistent.The main differences are located on the size of the background grids, as shown in Figure 4. Simulations with irregular waves require larger grid size compared to those for calm water in order to avoid strong numerical dissipation on wave propagation and to capture irregular waves accurately, where 30 grid points for the shortest wave length were used in the x direction and 75 grid points near the free surface were applied in the z direction according to the ITTC recommendation [30].The background grid size of the simulations with irregular waves will be evaluated by the generation of irregular waves as described in Section 5.1.Furthermore, Section 5.2 conducted the convergence and uncertainty analysis of simulations with calm water, which studied the influence of different hull and appendage grid sizes and demonstrated that the medium grid size can be applied for the following studies.
The hull grids and appendage grids are presented in Figure 5, where the boundary grids of the hull and fairwater were set small enough (Y+ ≤ 1) to accurately capture the flow around the boundary layer.The whole grid after overset around the model of simulations with calm water and irregular waves is presented in Figure 6, and it can be observed that the grid sizes along the x and z directions for simulations with/without irregular waves have obvious differences, where the cases with irregular waves show denser refine regions compared with those of calm water.

Verification and Validation
To evaluate the reliability of the numerical method, V&V (verification and validation) procedures proposed by The American Society of Mechanical Engineers, 2009, 2016 [31,32] have been conducted.The experimental uncertainties, numerical uncertainties, and input parameter uncertainties are used to evaluate the validation uncertainties, and the comparison between the experimental data and the numerical results shows the numerical errors.
The experimental uncertainty, UD, is assumed to be 1.00% D due to the lack of references provided, where D refers to the experimental data.The numerical uncertainty, USN, consists of the grid uncertainty, UG, the time step uncertainty, UT, and the iteration uncertainty, UI; however, the iteration uncertainty, UI, can be neglected.The input parameter uncertainties are assumed to be zero.
Due to the lack of experimental data about submarine sailing near the free surface with irregular waves, the overall V&V procedures are divided into three parts.First, the wave generation technique was validated by comparing the generation of irregular waves with the theoretical wave pattern and spectrum.Additionally, the verification study was conducted by grid/time-step convergence studies related to the simulation conditions of a submarine sailing near the free surface in calm water, which gave access to the numerical uncertainty, USN.Finally, the experimental data of the total resistance of a submarine sailing near calm water was compared with the simulation results, which provided the validation uncertainty, UV, and the numerical errors were evaluated according to the validation uncertainty, UV.

Wave Generation Technique
The validation of irregular wave generation has been conducted by comparing simulation results with theoretical values.The sea state 4 with a significant wave height of 2.5 m and a mean period of 7.2 s is selected to be validated and it has been scaled according to the scale ratio λ = 24 and Fn = 0.336.The irregular waves consist of N = 50 components, ranging from 0.455 Hz to 3.736 Hz.Meanwhile, the simulations were performed with 25,000 time steps, lasting 100 mean periods.The monitor point was set at X = 1 in the longitudinal direction, where it was located on 2L downstream of the inlet boundary.The statistical characteristics of random variables, such as the expected value (EV), standard deviation (SD), cumulative distribution function (CDF), and probability density function (PDF), are evaluated as: where Ji refers to the ith value of random variables, The time histories, EV, and SD corresponding to the wave elevation of theoretical values and simulation results are shown in Figure 7a, where the theoretical values are drawn according to different discretized linear wave components.To evaluate the deviations between the simulation results and theoretical values, the encounter energy spectra of the theory and simulation are acquired by FFT (Fast Fourier transform) analysis, which are shown in Figure 7b.The encounter circular frequency, ωe, for the heading waves is defined as: The nth order moment of the spectrum is defined as: The 0th order moment, m0, of spectrum, ωe, and S(ωe) corresponding to the peak point of the spectrum are shown in Table 3, where m0 stands for the energy of irregular waves.The differences between the CFD and theory are acceptably small, which indicates that the irregular wave generation technique is feasible and reliable for the viscous URANS solver.In the meantime, the cumulative distribution function (CDF) and probability density function (PDF) of the wave elevation are shown in Figure 7c,d

Uncertainty of the Numerical Method
The simulation of submarine sailing near calm water with Fn = 0.336 and h/D = 1.1 has been used to conduct the V&V study.The experimental data is quoted from Dawson.Due to a lack of experimental data for submarine sailing near the free surface with irregular waves, the V&V study is not related to those simulation conditions.
The numerical uncertainty USN is defined as follows: As mentioned before, the iteration uncertainty, UI, can be neglected.Grid/time-step convergence studies are carried out to obtain the grid uncertainty, UG, and the time step uncertainty, UT.The convergence study method proposed by Xing and Stern [33] has been applied and the convergence ratio, R, is defined as: where ε 32 = S f S m and ε 32 = S m S c .S f , S m , and S c refer to the solutions of fine, medium, and coarse, respectively.For different values of R, the convergence conditions are separated as: (1) 0 < R < 1, Monotonic convergence; (2) R < 0, Oscillatory convergence; (3) R > 1, Divergence.
When it achieves the Monotonic convergence, the Richardson extrapolation method is used to obtain the estimated order of accuracy, PRE, which is defined as: where r refers to the refinement ratio of the grid or time step.The distance metric P is defined as the ratio of PRE to Pth, where Pth is generally set as 2. The factor of safety method uncertainty is defined as: For the grid convergence study, three sizes of grids have been applied with the medium time step and the grid refinement ratio rG = √2.The coarse, medium, and fine grids are 1.96, 5.55, and 15.97 million grid nodes, respectively.When related to the time-step convergence study, three sizes of time steps have been applied with the medium grid and the time-step refinement ratio rG = 2.The coarse, medium, and fine time steps are 0.006, 0.003, and 0.0015, respectively.The specific number of grids for different parts is listed in Table 4.As mentioned before, the Richardson extrapolation method is used to obtain the numerical uncertainty, USN, due to the Monotonic convergence achieved.The results are shown in Table 5.The uncertainties of the grid and time step are 4.61% D and 0.42% D, respectively; the total numerical uncertainty is 4.63% D. To validate the numerical simulations, the deviations between simulation results and experimental data are supposed to be evaluated.The error, E, is defined as: where S refers to the simulation results.As mentioned before, the validation uncertainty, UV, is defined as: where UD is assumed to be 1.00% D due to the lack of any available references.The results of validation uncertainty are shown in Table 6 and the error, E, is lower than the validation uncertainty, UV, which indicates the simulation results are validated under the UV level.Accordingly, the medium grid and medium time step are selected to conduct the following numerical simulations, considering the reliable accuracy and high simulation efficiency.The comparison between CFD and EFD for the total resistance coefficient of a submarine sailing near calm water with different velocities and submerged depths is shown in Figure 8.The mean error of the total resistance coefficient between CFD and EFD is 2.33%, indicating the CFD results show a good agreement with the experimental data, where the total resistance coefficient is defined as follows: where Rt refers to the total resistance, ρ refers to the density, v refers to the flow velocity, and S refers to the wetted-surface area of the model.

Results and Discussion
In this section, the hydrodynamic performance of the model sailing near the free surface in calm water or irregular waves is presented and compared with each other.The resistance, lift, and pitch moment of the model are discussed using a statistical analysis method in the time domain and a FFT analysis in frequency domain; finally, some qualitative conclusions are drawn.

Calm Water
The resistance of the model sailing near the free surface in calm water and Fn = 0.336 at different submerged depths (h) ranging from 1.1D, 2.2D, and 3.3D is shown in Figure 9.The different components of resistance are divided into total resistance (Rt), residual resistance (Rr), and friction resistance (Rf) based on a Froude assumption.The values of Rt are 64.67N,54.68N, and 52.53N for h = 1.1D, 2.2D, and 3.3D, respectively.The values of Rt decrease 15.45% with h changing from 1.1D to 2.2D and decrease 3.93% with h changing from 2.2D to 3.3D.This indicates that Rt decreases with an increase of the submerged depth and the changes become less pronounced with an increasing submerged depth.Meanwhile, the values of Rr are 19.65N,10.11N, and 8.10N for h = 1.1D, 2.2D, and 3.3D, respectively, in which the values of Rr decrease 48.55% with h changing from 1.1D to 2.2D and decrease 19.88% with h changing from 2.2D to 3.3D.The values of Rf are slightly different with the changed submerged depth, which are 45.02N, 44.57N, and 44.43N for h = 1.1D, 2.2D, and 3.3D, respectively.The trends of Rr and Rf reveal that the total resistance is mainly dominated by the residual resistance, where Rr shows considerable changes but Rf remains relatively constant.The residual resistance appears to be more sensitive to the change of submerged depths compared with the total resistance.The free surface and wave pattern on the XOZ plane at different submerged depths ranging from 1.1D, 2.2D, and 3.3D are shown in Figure 10.The position of the crest and trough around the fairwater performs similarly for the three depths, where a crest and a trough are located before and after the fairwater, respectively.The position of the crest and trough after the fairwater changes with the submerged depths, where the crest and trough have an advanced phase for shallower submerged depths compared with the deeper one.The amplitudes of the waves decrease rapidly with an increasing submerged depth, which causes the decrease of wave-making resistance with a similar tendency.

Irregular Waves
The direct effects of irregular waves on the resistance of the model perform as the fluctuations with solution times, which indicate the first-order response caused by the wave-induced forces in fact.The time histories, EV, and SD of the total resistance (Rt) of the model sailing near the free surface with a sea state 4 and Fn = 0.336 are shown in Figure 11a-c, and the CDF and PDF of the total resistance at different submerged depths are shown in Figure 11d,e.When t = 149 s, the values of EV are 64.64N,55.85N, and 54.52N for h = 1.1D, 2.2D, and 3.3D, respectively.The values of EV decrease 13.60% with h changing from 1.1D to 2.2D and decrease 2.38% with h changing from 2.2D to 3.3D, which indicates that the influence of the free surface decreases rapidly at first and then tends to stay fairly constant with an increasing submerged depth.The steeper CDF curve, narrower peak value bandwidths of PDF, and the lower center value of the peak value bandwidth can be observed with increasing submerged depth, and the distribution of the total resistance becomes more concentrative accordingly, which leads to the decrease in EV of the total resistance and reflects the decrease of the first-order response amplitudes.Although the EV of the total resistance tends to be fairly constant, the total resistance is still considerable in the first-order response with the increase of submerged depths, which complicates the control strategies when the submarine sails near the free surface with irregular waves.To evaluate the effects of the free surface and irregular waves on the resistance components, the EV values of different resistance components at different submerged depths are shown in Figure 12.It can be observed that the Rr decreases with the increasing submerged depths and the Rf shows negligible deviations (<0.1%) for different submerged depths.When t = 149 s, the EV values of Rr are 19.92N,11.13N, and 9.86N for h = 1.1D, 2.2D, and 3.3D, respectively.The values of EV decrease 44.13% with h changing from 1.1D to 2.2D and decrease 11.41% with h changing from 2.2D to 3.3D.It can be concluded that the changes in residual resistance are dominant compared to the friction resistance.Due to the residual resistance being smaller than the friction resistance, the extent of the total resistance changing with submerged depths is weaker than the residual resistance.The free surface at different submerged depths changing with solution times is shown in Figure 13.It presents the wave pattern which is constituted by initial irregular waves and surface waves generated by the model.In the beginning of the solutions, the feature of the free surface mainly performs as the irregular waves.With the increase of the solution times, the wave generated by the model gets more and more obvious and the initial irregular waves are disturbed.When the solution time t = 2.975 s, the surface wave is combined with irregular waves and a Kelvin wave generated by the model, and the wave-making characteristics of the model are similar to those of calm water.The amplitude of the wave generated by the model decreases dramatically with an increase in submerged depth, which reveals the decrease of wave-making resistance.The EV values of different resistance components for calm water and irregular wave conditions at different submerged depths are shown in Figure 14.It can be observed that the resistance has slight deviations for different surface conditions at the same submerged depth.The added resistance which is caused by the irregular waves is not obvious, and the reason can be owing to the captive model.According to studies carried out by Sadat-Hosseini et al. [34,35] and Guo et al. [36], which simulated different models (KVLCC2 & KCS) sailing with and without a fixed degree of freedom, explored the mechanisms regarding the generation and change of added resistance.The results show that the added resistance stays relatively stable when the motions of the ship are fixed with different wave lengths, and the motions will greatly affect the added resistance.Herein, due to the restricted motions of the model, the added resistance component which originated from the motions of the model has not been considered.Furthermore, the added resistance that is related to the flow around the submarine model is the only component.It can be deduced that the added resistance related to the flow filed around the model is only somewhat relevant to the surface wave.In current simulations, the sea state 4 irregular waves consist of 50 different linear regular waves.To evaluate the characteristics of the resistance response for different component waves in the frequency domain, the FFT analysis has been conducted to acquire the resistance response at different submerged depths and the spectral responses of the resistance are shown in Figure 15.The comparison between the different submerged depths is shown in Figure 16.According to Figure 15, the amplitudes of the resistance response are different than those components of regular waves.The response spectrum has three clear, main peaks at different encounter frequencies, and the first main peak occurs at 1.064 rad/s (0.169 Hz), which is fairly near the iterative convergence frequency 1.16 rad/s (0.185 Hz) for the resistance of the model sailing near the free surface in calm water.The second main peak has a maximum amplitude among the response spectrum, in which the encounter frequency is about 5.028 rad/s (0.800 Hz).According to the encounter wave spectrum, the corresponding encounter frequency of the peak point is 5.516 rad/s (0.878 Hz).The encounter frequency at the third main peak is equal to 9.540 rad/s (1.518 Hz), which is approximately twice the value compared with the encounter frequency at the second main peak.The amplitude of the third main peak is much smaller than that of the second main peak, and the amplitude shows smaller values with an increase of encounter frequencies.The general tendency of the resistance response spectrum is similar to the incident wave spectrum, i.e., the response in the low frequency region produces considerable amplitude and the amplitude of the response remains quite small in the high frequency region.
According to the Figure 16, the frequencies of the response spectrum are fairly equal to those of the incident wave spectrum (which can be found in Figure 7).The amplitudes of the response spectrum decrease with an increasing submerged depth, in which the amplitude at the first main peak decreases 39.16% with h changing from 1.1D to 2.2D and decreases 39.55% with h changing from 2.2D to 3.3D, and the amplitude at second main peak decreases 76.96% with h changing from 1.1D to 2.2D and decreases 74.58% with h changing from 2.2D to 3.3D.The decrease is at a similar level between the several submerged depths, which is different from the trend of the EV of the total resistance as mentioned before.The general characteristics of the response spectrum are similar at different submerged depths, with the main deviations being reflected in the amplitudes of the response spectrum.

Calm Water
Similar to the chapter of resistance, the lift and pitch moment (which are related to the center of gravity) of the model sailing near the free surface in calm water and Fn = 0.336 at different submerged depths (h) ranging from 1.1D, 2.2D, and 3.3D are shown in Figure 17.The values of lift are 36.92N,7.59N, and 1.88N for h = 1.1D, 2.2D, and 3.3D, respectively.Additionally, the values of the pitch moment are 16.58N•m, 1.95 N•m, and 0.19 N•m for h = 1.1D, 2.2D, and 3.3D, respectively.The values of lift decrease 79.44% with h changing from 1.1D to 2.2D and decrease 75.23% with h changing from 2.2D to 3.3D.Accordingly, the values of the pitch moment decrease 88.24% with h changing from 1.1D to 2.2D and decrease 90.26% with h changing from 2.2D to 3.3D.The lift and pitch moments are dramatically changed with an increase in the submerged depth compared to the resistance, where the former is even over 75% and the latter is less than 20%.The lift and pitch moments show greater sensitivity to the change of the submerged depth compared with the resistance.The trend of the lift and pitch moments changing with submerged depths is different compared with the trend of resistance, in which the lift and pitch moments are getting close to zero with an increase in the submerged depth, but the resistance decreases to a fairly constant value.The total pressure acting on the surface of the model can be divided into two parts, one that is related to the distance away from the free surface is defined as the hydrostatic pressure, and the other that is caused by the non-uniform velocity field of fluid and is defined as dynamic pressure.In this section, the pressure distribution around the model is acquired to evaluate the tendencies of the lift and pitch moments changing with the submerged depths.Due to the hydrostatic pressure being only related to the distance from the free surface, the pressure gradient between the upper and lower surfaces of the model is constant for different submerged depths.This component of the total pressure will not be discussed, and the dynamic pressure is only considered for assessing the changes of the lift and pitch moments.Figure 18 shows the dynamic pressure distribution around the model at different submerged depths, which is measured on the XOZ plane.The model is normalized by the hull length in the x direction, and the dynamic pressure is also normalized, in which Cp refers to the normalized dynamic pressure coefficient according to Equation (22).
where Pr is the dynamic pressure, ρ is the density of fluid, and v is the incident flow velocity.The dynamic pressure distribution has some obvious features, such as the bow and front of the fairwater having a high-pressure region, where there are stationary points of the velocity field.Meanwhile, the top of the fairwater and aft shoulder of the model are the low-pressure regions, where there are relatively high-velocity regions of flow filed and separations of flow occur dramatically.For different submerged depths, the characteristics of the dynamic pressure distribution have some discrepancies, and the areas and amplitudes of the high/low-pressure regions are different for those depths.
To evaluate the dynamic pressure distribution on the surface of the model, the curves that are intersected between the XOZ plane and the upper/lower surfaces of the submerged model are set to monitor the dynamic pressure distribution.Figure 19 shows the dynamic pressure distribution on the curves, which is divided into the upper surface (left figure) and lower surface (right figure).From Figure 19a, it can be observed that the dynamic pressure distribution along the hull length on the upper surface for h = 1.1D has larger fluctuations compared with the other submerged depths, and the dynamic pressure distribution tends to be consistent with the increase in submerged depths.As mentioned before, the regions around the fairwater and aft shoulder of the model produce obvious low and high dynamic pressure, which will make a difference to the hydrodynamic force and moment of the model.Before the fairwater, the dynamic pressure coefficient, Cp, is about 0.1 larger than the other submerged depths for h = 1.1D and Cp is fairly similar for h = 2.2D and 3.3D.In the middle of the model, the Cp increases with an increase of the submerged depth and the increase becomes less pronounced for deeper submerged depths.The Cp of h = 1.1D is about 0.2 less than other submerged depths for 0.35 < x/L < 0.45.Around the aft shoulder of the model, the trend of Cp is similar to the one before the fairwater, and the Cp of h = 1.1D is about 0.1 larger than other submerged depths.For Figure 19b, the dynamic pressure distribution along the model on the lower surface has slight discrepancies for different submerged depths, and the main difference appears in the middle of the model (0.25 < x/L < 0.65), which has been marked with the red dash box.In the middle of the model, the dynamic pressure increases with an increase in the submerged depth and the increase becomes less pronounced for deeper submerged depths, where the Cp of h = 1.1D is about 0.03 less than other submerged depths.According to the dynamic pressure distribution on the upper and lower surfaces of the model, the pressure difference between the upper surface and lower surface will decrease with an increase of submerged depth and the trend will gradually become less obvious.In the meantime, the distribution of the pressure difference along the model will induce a decrease of the lift and bow-up moments with an increase in the submerged depth.

Irregular Waves
Similar to the chapter on resistance, the time histories, EV, and SD of the lift and pitch moments of the model sailing near the free surface with a sea state 4 and Fn = 0.336 are shown in Figures 20a-c and 3.3D, respectively.The EV of the lift decrease 78.67% with h changing from 1.1D to 2.2D and decrease 71.35% with h changing from 2.2D to 3.3D, and the EV of the pitch moment decrease 74.71% with h changing from 1.1D to 2.2D and decrease 46.47% with h changing from 2.2D to 3.3D.The trends are similar to those of the resistance, and the difference is located with the amplitudes changing with submerged depths.The changes of the lift and pitch moments are dramatically larger than those of the resistance, which demonstrates that the lift and pitch moments are more sensitive to the changes in submerged depth.Meanwhile, the steeper CDF curve, narrower peak value bandwidths of the PDF, and the lower center value of the peak value bandwidth can be observed with an increase in the submerged depth, and the distributions of the lift and pitch moments become more concentrated accordingly, which leads to the decrease in the EV of the lift and pitch moments and reflects the decrease of the first-order response amplitudes.The first-order response has considerable regions compared to the relatively small EV for the lift and pitch moments, and it can be deduced that large amplitude heave and pitch motions will be generated and the control strategies will be complicated when a submarine sails near the free surface with irregular waves.To investigate the mechanism of influence of the surface waves on the hydrodynamic performance of the submerged vehicle, the snapshots of the free surface at different solutions times of h = 1.1D are shown in Figure 23, and the hydrodynamic characteristics of the model and dynamic pressure distribution along the hull at different solution times are shown in Figure 24.It can be observed that the position of the wave crest and trough changes with the solution time due to the transient irregular waves, which are different compared with the calm water conditions shown in Figure 18.For t = 6.247 s, the surface wave coupled positively with the wave generated by the model and caused an increase in wave amplitude around the fairwater compared with that of t = 6.842 s, which resulted in the increase of wave resistance and a considerable increase of Cp on the upper surfaces, but the Cp on the lower surfaces was only affected slightly.The pressure difference between the upper surface and the lower surface around the fairwater will increase, and the pressure difference around the aft shoulder is similar for those two moments, which induces a negative pitch moment related to the coordinate at this moment.For t = 6.545 s, the dynamic pressure acting on the middle of the upper surfaces increases dramatically and the pressure difference is larger compared with that of t = 6.842 s, causing the observed decrease in lift.The values of EV for the lift and pitch moments for calm water and irregular wave conditions at different submerged depths are shown in Figure 25.The values of EV for the lift and pitch moments of the model sailing near the free surface with irregular waves are generally larger than the conditions with calm water at different submerged depths, except the pitch moment for h = 1.1D.The differences between the conditions with irregular waves and calm water decrease gradually with an increase in the submerged depth for the lift, where the pitch moment performs oppositely to the trend observed for the lift.The response spectra of the lift and pitch moments acquired by FFT analysis at difference submerged depths are shown in Figures 26 and 27, which produce some similar characteristics compared with those for the resistance, such as the encounter frequencies of the response spectrum are corresponding to those of the encounter wave spectrum.Meanwhile, the first main peak occurs at the encounter frequency 1.064 rad/s (0.169 Hz), which is near the iterative convergence frequency of the lift and pitch moments of the model sailing near the free surface in calm water.The features of the second main peak of the response spectrum produce discrepancies for the lift and pitch moments at different submerged depths; for example, it occurs at 5.021 rad/s (0.799 Hz) for lift in the three different submerged depths, but it occurs at 6.279 rad/s (0.999 Hz) for the pitch moment in h = 1.1D and 2.2D and at 5.636 rad/s (0.897 Hz) for the pitch moment in h = 3.3D.Other peaks perform small values with the increase of encounter frequencies.The response amplitudes produce larger values compared with the EV of the lift and pitch moments for different encounter frequencies, especially for the pitch moment, where the amplitude of the second main peak is several times greater compared with the corresponding EV of the pitch moment at different submerged depths.The surface waves would affect the hydrodynamic force and moments, although with larger submerged depths, which complicated the control strategies of the vehicle sailing near the free surface.
The comparisons of the response spectra between different submerged depths are shown in Figure 28, where the amplitudes of the response spectrum in different encounter frequencies decreases with increasing submerged depths.The amplitudes at the second main peak of the lift response spectrum decrease 42.54% with h changing from 1.1D to 2.2D and decrease 40.05% with h changing from 2.2D to 3.3D, and the amplitudes at the second main peak of the pitch moment response spectrum decrease 48.86% with h changing from 1.1D to 2.2D and decrease 47.15% with h changing from 2.2D to 3.3D.As mentioned above, the general characteristics of the response spectrum are similar at different submerged depths; the main deviations are reflected in the amplitude of the response spectrum.Compared with the encounter wave spectrum, the amplitudes of the response spectrum have considerable differences to those of the encounter wave spectrum in highfrequency regions, which demonstrates that the response produces nonlinear characteristics compared with the encounter wave spectrum, and the hydrodynamic force and moments acting on the model are less sensitive to high frequency excitations.

Conclusions
The numerical simulations of a SUBOFF model sailing near the free surface in calm water and long-crested waves at different submerged depths using URANS codes have been presented in this study, and the comparison between the different surface conditions at different submerged depths has also been conducted, which is aimed to investigate the influence of irregular waves on the hydrodynamic performance of the model at different submerged depths.
The results reveal that the EV of the resistance and lift decreases dramatically with an increase in the submerged depth, and the bow-up moment increases as the model gets closer to the free surface, in which the amplitudes of the surface waves decrease accordingly, indicating the influence of the free surface decreases with increased submerged depth.The trends of the hydrodynamic performance changing with submerged depths are similar to those of previously published research [8,19].
Due to the effects of irregular waves, the hydrodynamic force and moment of the model fluctuate with the solutions time.The amplitudes of the first-order response decrease with an increase in the submerged depth and produce considerable values compared with the relatively small values of EV at the deeper submerged depths, which will induce large amplitude motions and complicate the control strategies of submarine sailing near the free surface with irregular waves.In the frequency domain, the response spectrum frequencies of the hydrodynamic force and moments acquired by FFT analysis correspond to the encounter frequencies of the wave spectrum.Nonetheless, the peak characteristics of the response spectra are different to those of the encounter spectrum, where the response spectra produce several peaks in the low frequency domain and slight responses in the high frequency domain, demonstrating that the response produces nonlinear characteristics compared with the wave spectrum, and the hydrodynamic force and moments acting on the model are less sensitive to high-frequency excitations.

Figure 1 .
Figure 1.Geometry of the SUBOFF model with the fairwater.

Figure 2 .
Figure 2. General view of the simulation conditions.

Figure 3 .
Figure 3. Computational domain and boundary conditions of the model.

Figure 6 .
Figure 6.Grid around the model after overset for calm water and irregular waves conditions: (a) calm water; (b) irregular waves.
, respectively.All of the time histories, EV, SD, PDF, and CDF agree well with each other for theory values and simulation results.

Figure 7 .
Figure 7.Comparison between the CFD results and the theoretical values for the generation of longcrested waves: (a) time histories; (b) wave spectrum; (c) CDF; (d) PDF.

Figure 8 .
Figure 8.Comparison between CFD and EFD for the total resistance coefficient of a submarine sailing near calm water with different velocities and submerged depths.

Figure 9 .
Figure 9.The resistance of the model sailing near the free surface in calm water at different submerged depths.

Figure 10 .
Figure 10.The free surface at different submerged depths: (a) wave pattern; (b) wave amplitude.

Figure 12 .
Figure 12.The EV of different resistance components for the model sailing near the free surface with irregular waves at different submerged depths.

Figure 13 .
Figure 13.The snapshot of the wave pattern in the very beginning of the solution times: (a) 1.1D; (b) 2.2D; (c) 3.3D.

Figure 14 .
Figure 14.The comparison of the total resistance between calm water and the EV of irregular wave conditions at different submerged depths.

Figure 15 .
Figure 15.The spectral response of resistance for the model sailing near the free surface with irregular waves at different submerged depths: (a) 1.1D; (b) 2.2D; (c) 3.3D.

Figure 16 .
Figure 16.The comparison of the total resistance spectral response at different submerged depths.

Figure 17 .
Figure 17.The lift and pitch moments of the model sailing near the free surface in calm water at different submerged depths: (a) lift; (b) pitch moment.

Figure 18 .Figure 19 .
Figure 18.The dynamic pressure distribution around the model on the XOZ plane at different submerged depths.
and 21a-c, and the CDF and PDF of the lift and pitch moments at different submerged depths are shown in Figures 20d,e and 21d,e.The values of EV for the lift and pitch moments at different submerged depths are shown in Figure 22.When t = 149 s, the EV of the lift are 44.67N,9.53Ns and 2.73N and the EV of the pitch moment are 16.25 N•m, 4.11 N•ms and 2.20 N•m for h = 1.1D, 2.2Ds

Figure 23 .
Figure 23.The free surface at different solution times.

Figure 24 .
Figure 24.The hydrodynamic forces, moments, and dynamic pressure distribution along the hull length at different solution times: (a) hydrodynamic force and moments; (b) upper surface; (c) lower surface.

Figure 25 .
Figure 25.The comparison of the lift and pitch moments between calm water and the EV of irregular wave conditions at different submerged depths: (a) lift; (b) pitch moment.

Figure 26 .Figure 27 .Figure 28 .
Figure 26.The spectral response of the lift of the model sailing near the free surface with irregular waves at different submerged depths: (a) 1.1D; (b) 2.2D; (c) 3.3D.

Table 1 .
Main parameters of the SUBOFF model.

Table 2 .
Description of the simulation conditions.

Table 3 .
The 0th order moment, m0, of the spectrum, ωe, and S(ωe) corresponding to the peak point of the spectrum.

Table 4 .
Details for different grid nodes.

Table 5 .
Verification study for the total resistance.

Table 6 .
Validation study for the total resistance.