Numerical Simulation Study on the Effects of Course Keeping on the Roll Stability of Submarine Emergency Rising

Featured Application

The results in this research can provide technical support for the rudder control of a real submarine to suppress the excessive roll of emergency rising. The results also laid the foundation for further research of a free-running submarine model.

Abstract

A direct numerical simulation method based on Reynolds Average Navier–Stokes (RANS) equations is used to carry out numerical prediction studies of submarine emergency rising in this paper. Firstly, a numerical simulation of the nonpropelled model without rudder manipulating is accomplished as the basis of this study. The numerical prediction results are in good agreement with the experimental data, which proves the feasibility and accuracy of the direct numerical simulation method. Meanwhile, both model tests and numerical simulation results reveal the strong coupling effect between roll and yaw motions during the underwater ascending process. Based on the above observation and analysis, another two numerical simulations, whose grids are identical with the non-manipulation simulation, are conducted respectively under the condition of rudder steering, i.e., course keeping simulation and self-propulsion simulation. An optimized S surface controller based on conditional determination is designed to manipulate the rudders. As a result, the yaw angle of the latter two simulations is limited within the range of 0.2° and 0.8° respectively, proving the effectiveness of the S surface controller. Correspondingly, the maximum roll angle is reduced by 96% and 70% respectively, which demonstrates that the roll stability is significantly enhanced by improving the course keeping ability of the model. Moreover, it is also proven from the perspective of reverse verification that, the excessive yaw deviation is the root cause of emergency rising roll instability for the situation of incidence angle lower than 30°.

1. Introduction

Implementing the emergency ascent operation is the most critical measure for sailing underwater submarines to survive from emergencies, e.g., flooding, sternplane jam, even escaping from torpedo attack. The conventional operation procedure is to command maximum propeller spin rate and blow the forward ballast tanks to acquire an appropriate pitch angle. Meanwhile, sternplane control could be used to accelerate the pitching process if possible [1]. In this typical scenario, operators concentrate on getting to the surface, while ignore other parameters such as heading, and roll. Nonetheless, both model tests and full-scale trials indicate that an excessive roll, as well as serious roll oscillations, may occur when the submarine emerges through the surface, which threatens its safety. For instance, Itard [2] proved that a large roll angle of 60° appears before the free-running model getting to the surface. Moreover, Watt [3] also showed how important other variables e.g., yaw and roll etc., are to the stability of the ascent motion, henceforth referred to as “rising stability”.

In the study on emergency rising stability, previous researchers mainly focused on obtaining the hydrodynamic coefficients through physical captive model tests or numerical simulated captive model tests, then combining the 6 degrees of freedom (DOF) motion equations to carry out maneuverability prediction. Feldman et al. [4,5,6] performed a series of captive model experiments to investigate the stability and control characteristics of submarines, such as straightline towing, rotating arm test etc. Pan et al. [7,8] conducted the numerical prediction of submarine hydrodynamic coefficients and rotary derivatives using the Computational Fluid Dynamics (CFD) approach based on the standard equations of motion for submarine maneuvers [9]. In addition, Booth [10] carried out a classical stability analysis of a nonpropelled submarine by introducing nonlinear incidence terms to conventional coefficient modeling. Booth pointed out the lateral instability that is caused by excessive buoyancy, and there existed a strong couple effect between roll and yaw according to the experimental results. However, as he acknowledges, his quasi-steady stability analysis did not model the unsteady forces associated with flow separation and vortex shedding [11]. Haarhoff and Sharma [12] provided a good description of the mathematical and physical insight for the horizontal plane stability parameters. Furthermore, Watt et al. established a quasi-steady coefficient model for studying roll stability based on experimental data from the wind tunnel tests [13,14,15,16]. It was manifested that the unsteady effects could be approximated with information gleaned from steady state experiments [13]. In his rolling stability analysis, the destabilizing hydrodynamic influence of the sail was considered as the main source of roll instability. Nevertheless, this conclusion cannot explain the early onset of the instability in several diesel-electric submarine trials. Based on his previous work, a 3 DOF stability analysis had been further accomplished, confirming the existence of yaw instability, whereas there were some unexplained trial instabilities remaining. Park et al. [17] conducted a series of model tests to investigate hydrodynamic coefficients of submarines at high incidence angles, laying the foundation for analyzing rising stability.

Given the above information, the stability analysis based on quasi-steady hydrodynamic coefficients has an insufficient capacity of modeling the unsteady viscous effect associated with trailing sail vorticity [3,10,13]. Furthermore, these hydrodynamic coefficients obtained from captive model tests or empirical estimates are merely applicable to infinite depth flow fields, without considering the influence of free surface or surfacing motion on the coefficients. Thus, the prediction accuracy of a coefficient-based model needs a further verification. In addition, the wind tunnel tests are set up with static steady conditions, with ignoring the dynamic response of the model during emergency ascent. Therefore, it is difficult to predict the roll oscillations of a model on the surface, as well as some excessive roll that occurred in the full-scale trials.

Overall, the emergency rising process involves highly nonlinear, complex coupled motions, especially, as the boat becomes lighter, or it speeds up, the response will change. More importantly, the surfacing process of the model is difficult to model based on relevant hydrodynamic coefficients. However, with the development of high performance computing technology, CFD has been a powerful alternative tool that provides direct numerical solutions of submarine space maneuvers, with accounting for unsteady viscous effect [18,19]. Meanwhile, considerable details about flow characteristics can be obtained through CFD methods with lower costs and higher accuracy, which provides insight into the physics of the problem. In light of the above advantages of CFD methods, various CFD simulation capabilities for submarine maneuvers had been developed by hydrodynamics researchers [20,21,22,23], especially the Reynolds Average Navier–Stokes (RANS) methods coupled with overset grid technology [24,25]. McDonald et al. [20] provided a physics-based computational method for the trajectory prediction of the fully appended SUBOFF model with a propeller. They interactively coupled the numerical solution of the three-dimensional (3D) unsteady RANS equations, which is solved using a multiblock multigrid scheme, with a 6-DOF motion solver. Pankajakshan et al. [21] presented 6 DOF RANS simulations validated by steady state and radio controlled free swimming model experiments. Their work achieves the direct simulation of hull-propeller-rudder interaction by solving the unsteady RANS (URANS) flow about a submarine. In addition, Zaghi et al. [24] applied dynamic overlapping grids to simulate the flow around a fully appended submarine, while Martin et al. [25] adopted two approaches to simulate maneuvers of a radio-controlled submarine model, i.e., the direct simulation using dynamic overset technology and the coupled CFD/potential flow propeller solver. Regarding the study of rising stability of submarine using 6 DOF RANS solver, Mark et al. [26] took the lead to employ CFD method to investigate the rising stability problem, and developed an unsteady 6 DOF RANS capability for simulating submarine rising maneuvers. He established a submarine model with fixed appendages and no propeller, while the hull-propeller-rudder interactions were modeled through some approximation methods. Meanwhile, the submarine model adopted in his study was an envelope model without any permeable regions or flooding holes in the superstructure. Thereby, Mark was unable to take the influence of the drainage process on surface stability into consideration, such that his study simply predicted the behavior of submarine before surfacing. To sum up, to date, it appears that no one has applied the CFD method to study the rising stability of an ascending submarine with moving rudders and a rotating propeller. In the present work, a 6 DOF URANS method coupled with the overset grid technique is developed for the direct numerical simulation of self-propulsion submarine emergency rising.

On the other hand, model tests continue to be the most reliable means to validate the accuracy of CFD-based simulations due to the insufficiency of turbulence modeling ability of CFD methods. With regard to experiments that investigate the roll stability of emergency rising, Zhang et al. [27] validated the coupling effect between yaw and roll motions by implementing a set of model tests of submarine emergency rising. Zhang pointed out that the excessive yaw deviation is the key factor resulting in excessive roll during buoyant ascent. Zhang believed that the excessive roll could be avoided effectively by minimizing yaw deviation and then limiting the increase of drift angle. However, due to the limitations of the experimental equipment capabilities, Zhang et al. did not perform the model test with rudder manipulation to verify his inference, let alone self-propelled model tests.

In order to verify the assumptions about yaw instability and further explore the influence of propeller on roll stability, the author employs the CFD method to carry out the numerical simulation of submarine emergency rising maneuvers. A self-propulsion simulation with moving rudders and a rotating propeller is accomplished by using Star-CCM+ solver, combined with overset grid technique and S surface control method. The effects of rudder and propeller on roll stability are analyzed to provide technical support for submarine emergency rising maneuvers. Therefore, the structure of this paper is organized as follows: In Section 1, a brief background about the model test is presented firstly, followed by the preliminary analysis of the strong coupling effect between roll and yaw motions, which lays the foundation of this study. Section 2 describes the numerical simulation method in detail established in this paper, including modeling analysis and hydrostatic attribute calculation, meshing methodology for the small gaps, superimposed rotation of moving objects, etc. In Section 3.1, an equivalent simulation with the same state as the model test is conducted, while the accuracy of the numerical simulation is verified by comparing it with the experimental data in Section 3.2. Section 3.3 introduces an S surface controller based on the conditional determination mechanism firstly. Then, using an identical grids model, a nonpropelled model simulation with rudder steering is carried out to analyze the influence of course keeping on the roll stability. Section 3.4 further accomplishes the self-propulsion simulation to preliminarily explore the effect of the propeller on its course keeping capability and roll stability. Finally, the Section 4 summarizes the innovative conclusions of this paper, and proposes the next research plans.

2. Methods: Experiment and Numerical Simulation

2.1. The Model Tests Method

According to the public data of diesel-electric submarine at home and abroad, a medium-sized diesel-electric submarine usually has a relatively high sail. This is hydrodynamically destabilizing because, as the submarine begins to roll, the rolling moment increases in the same direction, thereby rolling the submarine further. That is to say, submarines are prone to experience excessive roll, which greatly threatens the safety of submarine equipment and personnel. In order to investigate the emergency rising characteristics of these diesel-electric submarines, the Science and Technology on Autonomous Underwater Vehicle Laboratory of Harbin Engineering University has conducted early research in the area of model test of submarine emergency rising. Referring to the public data of some mother ships, researchers design a specific high-sail submarine model for test purposes, meanwhile, a series of model tests were organized in the comprehensive test water pool. Based on the previous model, tests without any tailplanes deflections, this laboratory further implemented the emergency rising test under the conditions of fixed sternplane angles. Please refer to the reference [27] for detailed information about the experiment, which is beyond the scope of this paper. Here is a brief introduction to these model tests.

The model test is conducted in a 50 × 30 × 10 m water pool, and the layout of the experiment facilities is shown in Figure 1a, which mainly consists of a submarine model and its data acquisition system, a towing sled and its support bracket, and camera devices with lighting, etc. Before the start of a test, the submarine model was locked on the towing sled by an electromagnetic clamping mechanism and hoisted to the top of the support bracket in the pool. At the beginning of the tests, the variable-frequency motor is activated to accelerate the towing sled along the track. When the model speed is substantially stabilized at the target value, the electromagnetic release device is activated so that the submarine can be launched immediately and smoothly. During the rising process, the model posture is recorded by the embedded PC104 bus-based measurement system and then transmitted to the host computer via zero buoyancy cable. Meanwhile, a set of camera devices are employed to record snapshot videos for observing. As a result, an emergency rising model test is achieved by means of releasing the model with reserved buoyancy instantaneously.

Based on the above experimental principle, the researchers have carried out a large number of model tests for combinational conditions of various factors, including different blown position, speed, rudder angle, etc., and established a sufficient experimental database. In order to highlight the issues that this paper focuses on, 5 representative model tests as shown in

Table 1. The working conditions list of 5 representative model tests.

Test Number	d0 (m)	u0 (m/s)	${\mathit{x}}_{\mathit{b}\mathit{l}\mathit{o}\mathit{w}\mathit{n}}$ (mm)	${\mathit{\delta }}_{\mathit{s}}$ (°)
1	8	1.15	362	0
2	8	1.15	362	3
3	8	0.95	362	9
4	8	0.75	568	3
5	8	0.75	568	3

, are selected to illustrate the irrelevant influence of moderate incidence angle and further verify the accuracy of the numerical results as well. Some video snapshots of the model test 1 from Table 1 are presented in Figure 1b–d, which occurs excessive roll before the model broaches the surface.

2.2. Theory of Numerical Simulation Method

2.2.1. Description of Submarine Motion

The earth-fixed inertial coordinate system O_E—ξηζ and the body-fixed coordinate system O_B—xyz are defined respectively to describe the motion of submarine. As shown in Figure 2, each coordinate axis is determined according to the right-handed system. The state variable y is used to describe the speed and position of the submarine:

$$\mathbf{y}=\mathbf{y}(u,v,w,p,q,r,{\xi }_0,{\eta }_0,{\zeta }_0,\varphi ,\theta ,\psi )$$

(1)

Among them, $(u,v,w)$ and $(p,q,r)$ are the translational and rotational velocity of the submarine with respect to the x, y, z body-axes direction respectively, while $({\xi }_0,{\eta }_0,{\zeta }_0)$ is the position of the submarine with respect to the $\xi ,\eta ,\zeta$ direction of the inertial coordinate system. $(\varphi ,\theta ,\psi )$ is the Euler angles orientating the body with respect to the inertial coordinate system, and are respectively the roll, pitch and yaw angles. In order to determine the state variable y, it is necessary to solve the equation of motion for both the submarine and fluid regions simultaneously, which will be discussed in detail in Section 2.2.2 and Section 2.2.3.

Furthermore, referring to Watt’s research ideas, the hydrodynamic angles used in the stability analysis is defined as the Equation (2):

$$\{\begin{array}{l}\beta ={\tan }^{-1}(-v/u)\\ \Theta ={\tan }^{-1}(\sqrt{{(-v)}^2+{(-w)}^2}/u)\end{array}$$

(2)

2.2.2. Solid Body Equations of Motion

Based on Newton’s second law, applying the centroid motion theorem and the momentum theorem relative to the centroid, the 6 DOF solid body equations of motion (EOM) for the submarine with respect to the body axes can be expressed as Equation (3):

$$\{\begin{array}{l}\frac{d\mathit{B}}{dt}+\mathit{\Omega }\times \mathit{B}=\mathit{F}\\ \frac{d\mathit{K}}{dt}+\mathit{\Omega }\times \mathit{K}+\mathit{U}\times \mathit{B}=\mathit{M}\end{array}$$

(3)

The two right-hand-side terms of these equations, F and M, describe the sum of all external forces and moments applied on the submarine, which including the hydrodynamic forces, the static weight and buoyancy forces, and where applicable the appendage control forces and thrust forces. Besides, Ω = (p, q, r) and U = (u, v, w) are the velocity vector of the submarine as stated above, while B = (mu, mv, mw) and K = (I_xxp, I_yyq, I_zzr) are the momentum and the moment of momentum for the submarine respectively. Since the specific solid body equations of motion employed in this work have been derived in literature [9], we will not repeat it here.

2.2.3. Fluid Equations of Motion

The three-dimensional viscous incompressible flow field around the emergency rising submarine is determined by solving the URANS equations. The derive process of fluid equations of motion based on the local velocities results in apparent body forces appears on the right-hand side (RHS) of the transport equations. Therefore, the conservation equations for mass and momentum based on local velocities are, respectively [28]:

$$\frac{\partial \rho }{\partial t}+\frac{\partial \rho U_j}{\partial x_j}=0$$

(4)

$$\frac{\partial \rho U_i}{\partial t}+\frac{\partial \rho U_i{U}_j}{\partial x_j}=-\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}\right)\right]-\frac{\partial }{\partial x_j}\rho \overline{u_i^{\prime }{u}_j^{\prime }}+{\mathit{F}}_B$$

(5)

$${\mathit{F}}_B=-\rho \left(\frac{d^2{\mathit{R}}_0}{d{t}^2}+\frac{d\mathit{\Omega }}{dt}\times {\mathit{r}}_1+\mathit{\Omega }\times (\mathit{\Omega }\times {\mathit{r}}_1)+2\mathit{\Omega }\times U_i\right)$$

(6)

Among the above equations, the instantaneous fluid velocity is decomposed into the mean components U_i, U_j and fluctuating components u_i, u_j, where the subscripts identify corresponding components along i and j directions respectively. x_i and x_j are position coordinates of fluid particles along i and j directions. Additionally, ρ is the density, μ is viscosity, and P is the mean component of the instantaneous pressure. Furthermore, F_B is the apparent body force as a result of the translating–rotating frame of reference, and these terms in Equation (6), in order from left to right, account for apparent linear, angular, centripetal, and Coriolis accelerations sequentially. In this force, derived in [29], R₀ is the position vector locating the submarine with respect to the inertial coordinate system, Ω is the angular velocity vector of the body coordinate system, and r₁ is a vector that locates a fluid particle in body axes. Besides, the linear acceleration term can alternatively be written as follows when referenced to the body axes:

$$\frac{d^2{\mathit{R}}_0}{d{t}^2}=\widetilde{U^{\prime }}+\mathit{\Omega }\times \tilde{U}$$

(7)

where $\tilde{U}=u{\mathit{e}}_x+v{\mathit{e}}_y+w{\mathit{e}}_{\mathrm{z}}$, ${\tilde{U}}^{\prime }=\dot{u}{\mathit{e}}_x+\dot{v}{\mathit{e}}_y+\dot{w}{\mathit{e}}_{\mathrm{z}}$ are the velocity and accelerate velocity vector of the submarine respectively.

For turbulent simulations, the k-ω based Shear Stress Transport (SST k-ω) turbulence model has been selected to model the Reynolds stress term, $\rho {\overline{u}}_i{\overline{u}}_j$, and equations are shown as following forms:

$$\begin{array}{l}\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({\Gamma }_k\frac{\partial k}{\partial x_j}\right)+{\tilde{G}}_k-Y_k+S_k\\ \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({\Gamma }_{\omega }\frac{\partial \omega }{\partial x_j}\right)+{\tilde{G}}_{\omega }-Y_{\omega }+S_{\omega }+D_{\omega }\end{array}$$

(8)

In the above formula, k is the turbulent kinetic energy, ω is the specific dissipation rate of turbulence, Γ_k, Γ_ω are the effective diffusions of k and ω, respectively. Y_k, Y_ω are the dissipation of k and ω; ${\tilde{G}}_k$, ${\tilde{G}}_{\omega }$ are the generator terms of k and ω; D_ω is a cross-diffusion term; S_k, S_ω are the custom source terms.

According to the previous numerical investigations, it has been proven that the maximum prediction deviation between several turbulence models, such as SST k-ω, Realizable k-ε, does not exceed 5%, which has a negligible effect on the analyses of the emergency rising stability. Therefore, the comparative studies on turbulence models are ignored in this paper, and only the optimal turbulence model is selected for calculation, i.e., SST k-ω, which is also to establish a constant reference benchmark. The SST k-ω model has been shown to predict the onset of flow separation under adverse pressure gradients [28,29] without computation intensively. Moreover, the SST k-ω model solves the boundary layer with the aid of wall functions [28], which model the wall velocity profile when the near-wall y+ values are larger than 12, thereby avoiding fully resolving the viscous sublayer with a very large number of mesh elements.

In the numerical solution process, the continuum equation, momentum equation and turbulence equation are discretized by second-order upwind scheme, and the pressure-velocity coupling is solved iteratively using Semi-Implicit Method for Pressure-Linked Equations (SIMPLE) algorithm. The hydrodynamic forces and moments applied on the hull at different times are evaluated by integrating pressure and shear data along the surface. Combined with the 6 DOF equation of the submarine, the time-varying curves of the emergency rising motion are solved. In the meantime, the colorful cloud map of flow field information such as velocity, pressure and vorticity can be captured.

2.3. Description of the Research Model

This section focuses on the analysis of the static properties of the submarine model used in the numerical simulation. The self-built submarine model is shown in Figure 3. The main hull consists of a blunt bow, an axisymmetric parallel middle body, a parabolic stern, while the appendages including stern control surfaces with “+” form, the sail and its sailplanes. In order to study qualitatively the effects of propeller thrust and torque on the roll stability, a 7-blade highly skewed propeller is self-designed by referring to relative data of INSEAN E1619 propeller, which is commonly used to the wake behavior analysis of a generic submarine propeller [30]. Notice that the bow, stern, and deck casing section of the model have permeable regions with flooding holes, rather than a conventional envelope body. In other words, the present work takes the influence of flooding regions on the surfacing motion into consideration as truly as possible. Furthermore, since this paper focuses on the response characteristics of emergency rising motion after completing the blowing process, it is necessary to analyze the hydrostatic properties of the model with residual buoyancy.

With regard to the numerical simulation of emergency rising, there are three following assumptions to be made:

• Assuming that the time required for the submarine to blow the ballast tank is very short, the numerical simulation does not model the blowing process of pressure air. At the beginning of the calculation, the model obtains buoyancy force instantaneously, and the variable external force produced by the blowing process does not applied on the submarine.
• The model is released with a fixed sternplane angle at the desired speed, it will be not imposed the control force of sternplanes, except the propulsion forces and the rudder control force.
• Assuming there is no cavitation effect for the propeller because of its low revolution speed.

Based on the above assumptions, the weight of the submarine after blowing off ballast tank is calculated as follows:

$$m=m_{\nabla \downarrow }-m_{blown}-m_{casing}$$

(9)

where, $m_{\nabla \downarrow }$ is the weight of water that the external hydrodynamic envelope of the submarine displaces, $m_{blown}$ is the blown ballast mass and $m_{casing}$ is the displacement of the permeable region.

The position of the center of gravity of the submarine is calculated by:

$$G=\frac{m_{\Delta \downarrow }{G}_{\Delta \downarrow }-m_{blown}{G}_{blown}-m_{casing}{G}_{casing}}{m}$$

(10)

where, $G_{\nabla \downarrow }$ is the center of gravity of the submarine envelope displacement, $G_{blown}$ is the center of gravity of buoyancy load, $G_{casing}$ is the center of gravity of the permeable region.

The submarine moments of inertia are calculated by:

$$I=I_{\nabla \downarrow }-I_{blown}-I_{casing}$$

(11)

where, $I_{\nabla \downarrow }$ is the moment of inertia of the submarine envelope displacement, $I_{blown}$ is the moment of inertia of the buoyancy load, and $I_{casing}$ is the moment of inertia of the permeable region.

The origin of the body coordinate system O_B is set to the center of gravity of the submarine. Similar to the adding weight method, the model has residual buoyancy under water by reducing the weight of the model. Then calculate the center of gravity of the model that removes bow, stern, deck permeable regions, and blown ballast load, so does the center of buoyancy. Therefore, the mass centroid of the model in air varies with the change of longitudinal position of the blown ballast. By convention, the center of gravity $(x_G,y_G,z_G)$ of the model with residual buoyancy are defined with respect to (w.r.t.) the midship and model central axis, so is the position of buoyancy and blown ballast load. For instance, the mass centroid coordinates of blown ballast load in

Table 2. Main computational parameters of buoyant model and self-designed propeller.

Parameters	Value	Parameters	Value
L/m	5.09	$I_{xx}$/kg·m2	34.348
D/m	0.58	$I_{yy}$/kg·m2	1001.61
h/m	0.02	$I_{zz}$/kg·m2	1000.32
m/kg	862.80	$\Delta \downarrow$/kg	1107
mblown/kg	44.886	$(x_{blown},y_{blown},z_{blown})$/m	(0.362, 0, 0.01534)
Douter/m	0.3365	$(x_G,y_G,z_G)$/m	(0.14676, 0, −0.02824)
Dinner/m	0.0664	$(x_B,y_B,z_B)$/m	(0.16863, 0, −0.00823)

are (0.36258, 0, 0.015343), i.e., 0.36258 m forward from the midship and 0.015343 m upward from the central axis, respectively. Last but not the least, the other hydrostatic parameters (w.r.t. body axes) of the buoyant model that completing the blowing process are listed in Table 2. That is to say, the author only provides the final parameters, and ignores the pre-process of blown ballast load as expressed in Equations (9)–(11).

Note that the calculation parameters are given based on the model scale, the simulation results actually can be translated to a real submarine’s prediction results by the similitude rule, because the model tests are implemented in accordance with the similarity criteria. Specifically, the Re number of model tests is higher than the critical Re number, and for the mother ship and self-designed model, the flow similarity of flooding holes, the St number similarity, the Fr number similarity and geometry similarity are all satisfied as well. Therefore, the dimensionless parameters of a real submarine can be obtained by nondimensionalizing the corresponding parameters under the model scale. For instance, the blown mass and longitudinal position can be nondimensionalized with respect to L and $\nabla \downarrow$ respectively. Similarly, the predicted attitudes of the model reflect a real submarine’s motion, except that the time axis of attitude curves needs to be scaled based on the root-mean-square of the scale ratio.

2.4. Meshing Methodology

In the present work, all simulations are accomplished in the commercial software Star-CCM+ with employing multi-level overset mesh technique, and the computational domain is decomposed into several overset regions and background region by region-based meshing methodology. This multiple overset region consists of a main hull, rudder and propeller three grid blocks as shown in Figure 4b, in which each moving component has its own grid to deal with complex motion problems. Specifically speaking, the hexahedral grids are generated using the trimmer mesher for the background region, while polyhedral grids are generated using the polyhedral mesher for the main hull subregion because polyhedral grids can achieve a reasonable level of accuracy with the relative smaller number of grids and shorter computation time [31]. Besides, it is worth noting that trimmer mesh also is applied to the sternplanes and propeller subregions to deal with some small gaps better. Furthermore, prism layer mesher is utilized to gridding boundary layers of the no-slip wall for the employment of wall function. Prism layer mesh are generated such that average near wall Y+ values could be larger than 30, as required by the k-ω turbulence model.

In the study of naval ship computational fluid dynamics, the polyhedral meshes have been widely used in the computation of flow around objects with complex geometry due to its strong adaptability and high computational efficiency. Furthermore, the trimmed mesh of the hexahedral type has a relatively higher accuracy and is more suitable for the background mesh generation of multiphase flow with a free surface. In this paper, three meshing schemes, i.e., polyhedral mesh, trimmed mesh and hybrid mesh (polyhedral mesh for the overset area and trimmed mesh for background area), are analyzed for grid independence. Taking the pitch angle that is relatively stable in the emergency rising process as the comparison benchmark parameter, the grid number and calculation error of the three grid schemes are given in

Table 3. Mesh independence analysis summary.

Mesh Type	Background Cells (M)	Overset Cells (M)	Total Cells (M)	θmax (°)	Error (%)	θmin (°)	Error (%)
Model test	—	—	—	12.78	0.00	−3.10	0.00
Polyhedral mesh	1.156	3.241	4.397	11.92	6.73	−1.04	66.5
Trimmed mesh	0.864	3.795	4.659	11.35	11.2	−2.52	18.7
Hybrid mesh	0.864	3.241	4.105	11.83	7.43	−3.54	14.2

.

It can be seen from Table 3, that the polyhedral mesh is not suitable for the calculation of multiphase flow in this paper because of its large deviation in predicting the pitch angle of the submarine after surfacing. To the contrary, both the trimmed mesh and the hybrid mesh scheme have higher prediction accuracy before and after surfacing, and the hybrid mesh has the minimum perdition error for submarine underwater pitch angle. Therefore, the third hybrid mesh scheme is adopted in this paper.

The interface between each two overlapping regions is set to overlap boundary for flow information interpolation. Considering the calculation accuracy in “cutting hole” and “interpolation” processes, the grid dimensions of different grid blocks in overlapping areas should be similar to offer enough donor cells for interpolation. Therefore, the grid around submarine transitions from a refined prism layer to a coarse polyhedral mesh away from submarine to form an overset region as shown in Figure 4a,c. At the same time, several refinement regions are applied to the background region by volume control method. So that the grid dimensions of the background region increase gradually from the overlapping boundary to the far-field boundary. In addition, in order to capture the free surface of the two-phase flow, the grids at the free surface are refined through anisotropy methodology, with densifying along normal direction of the free surface. Eventually, the resulting meshed regions are obtained as shown in Figure 4.

In order to achieve self-propulsion simulation, an overset mesh approach with multiple regions overlapping are used to model the superposition motion including rudder and propeller rotation. Since overset grids require a minimum of three active cell layers for resolving a gap between two wall boundaries, there are some small gaps between rudders and stabilizing fins that require special treatment. Therefore, when coupling the overset mesh and the background mesh, we activate prism layer shrinkage to morph the cells of the prism layers in these small gaps, so that the above requirement is met. As shown in local enlarged view in Figure 4b, there are at least 8 layers prism cell in small gap, which improves the accuracy of the overlapped grids interpolation.

2.5. Boundary and Initial Conditions

The Star CCM provides solutions to the unsteady RANS equations discretized through the finite volume method (FVM), meanwhile, the volume of fluid (VOF) method is applied to capture the free surface interface between air and water. All boundary settings and initial conditions are depicted in

Table 4. The boundary settings and initial conditions.

Items	Boundary Condition	Initial Condition
Submarine	No-slip wall	Specified for different simulations
Behind side	Pressure outlet	The hydrostatic pressure of flat VOF wave
Remaining sides	Velocity inlet	The velocity of flat VOF wave
Water	—	The volume fraction of heavy fluid of flat VOF wave
Air	—	The volume fraction of light fluid of flat VOF wave

. Besides, a snapshot of the self-propulsion simulation with the rotating rudders and propeller is illustrated in Figure 5. This scenario is set up to observe the movement of rudder and propeller conveniently, and generate a simulation animation if necessary.

Based on the research ideas of verification and validation, the accuracy and feasibility of numerical simulation under the equivalent working conditions of model test are verified firstly. Then, using the same numerical simulation method and grids, another two simulations with rudder steering or propeller rotating is further carried out respectively, to verify the viewpoint of this paper. In summary, the numerical simulation conditions involved in this piece of work are listed in

Table 5. The working conditions of numerical simulation present in this work.

Simulation Condition	Rudder Steering	Propeller Rotating	d0 (m)	u0 (m/s)	${\mathit{x}}_{\mathit{b}\mathit{l}\mathit{o}\mathit{w}\mathit{n}}$ (mm)	${\mathit{\delta }}_{\mathit{s}}$ (°)
A	×	×	8	1.15	362	0
B	√	×	8	1.15	362	0
C	√	√	8	1.15	362	0

, with being divided into three categories, i.e., A, B, C. They are all simulated at the longitudinal blown ballast position of 362 mm where an excessive roll is prone to occur.

Taking self-propulsion simulation as an example, at the beginning of the simulation, the submarine is set to straight motion with a fixed depth of 8 m (the height of keel to the free surface), while the propeller revolutions per minute (rpm) is manually adjusted to obtain overall axial force equilibrium. The other 5 DOF motions are released when the initial velocity maintains around the target speed of 1.15 m/s at an appropriate propeller spin rates, while propeller speed remains constant in the remaining phase. After released, the submarine will gradually rise to the surface due to the buoyancy greater than the gravity, and continue to advance for a while on water surface under the action of propeller trust. The computation time step is 0.01 s, with 20 times iterations for each time step. Simulation results indicate average Y+ values increase from 30 at the beginning and reach to its peak value 80 when the submarine is moving at maximum axial velocity of 1.5 m/s.

3. Results and Discussions

3.1. The Coupling Effect of Yaw and Roll

In view of the fact that this paper mainly uses CFD as a means to carry out emergency rising simulations that cannot be achieved by the model test at present, this section firstly describes the coupling effect between yaw and roll found in the model test by analyzing the experimental data in detail. By comparing the maximum roll angle of several different test conditions whose incidence angle are almost the same, it is explained that the incidence angle is not the key factor to determine whether the model excessive roll or not, at least not at moderate incidence angles.

Based on the experimental data of test condition (1) given in

Table 6. The working conditions of numerical simulation present in this work.

Test Condition	Θmax (°)	θmax (°)	ϕmax (°)	tout (s)	ψout (°)	βout (°)	θout (°)	ϕout (°)	Vout (m/s)	Tϕ-ϕ (s)
(1).8-362-0-1.15	19.24	12.98	53.49	16.92	14.11	15.53	5.52	40.23	0.66	5.12
(2).8-362-3-1.15	19.53	13.54	50.87	16.13	13.61	13.54	8.93	29.43	0.87	5.13
(3).8-362-9-0.95	18.98	20.10	15.42	15.21	7.43	6.36	19.89	9.19	1.02	5.68
(4).8-568-3-0.75	20.37	26.13	24.16	14.33	8.94	8.90	25.19	17.66	1.20	5.63
(5).8-568-3-0.75	19.68	26.46	9.57	14.07	3.32	2.51	26.09	4.49	1.22	5.98

, the typical curves of attitudes for non-propelled emergency rising without rudder steering are depicted in Figure 6a. The code in the first column of Table 6, taking 8-362-0-1.15 for example, indicates the test condition of d₀ = 8 m, $x_{blown}$ = 362 mm, ${\delta }_s$ = 0° and u₀ = 1.15 m/s. Apparently, the model occurs excessive roll and obvious yaw deviation simultaneously, while the pitch curve does not have a plateau stage. Meanwhile, there is an obvious positive correlation among roll, yaw and drift angles before surfacing. On one hand, as the model occurs roll, the v component of local vertical ascent velocity increases rapidly while the w component decreases correspondingly, such that the drift angle increases with the increase of the roll angle, and further results in an evident yaw deviation. On the other hand, the incipient increase of drift angle may be induced by the growth of yaw deviation, as shown at around t = 10 s, due to the lack of course keeping capability. The yaw deviation emerges at the earlier beginning and continues accumulating until the model rises to the surface. During the underwater rising process, the drift angle grows with the increases of yaw deviation, further causing excessive roll. Thus, it is hard to tell who induces whom exactly. In other words, there are complex coupling effects between yaw and roll motions, considering the drift angle as the medium.

The comparison curves of maximum roll angle for working conditions listed in Table 6 are drawn in Figure 6b, where the x coordinate is labeled the test number 1–5 from left to right. It can be seen from this subgraph that there is an obvious linear correlation trend among the three curves of roll angle (magenta dashed line), yaw angle (red solid line) and drift angle (blue dot dash line) on the premise that the incidence angle (black dotted line) is approximately 20°. Along with yaw deviation decreasing from 15° to smaller than 5°, the maximum incidence angle gradually decreases from 60°to lower than 30°. Meanwhile, the curves of yaw and drift angle approximately coincide, which demonstrates that it is yawing motion that causes the drift angle to change synchronously.

The postures statistics of the test scenarios with the maximum incidence angle of about 20° among the model tests are listed in Table 6. It is remarkable that under the premise of incidence angle approximately equal, the maximum roll angle ϕ_max of the test conditions (1) and (2) is larger than 50°, with occurring the excessive roll, while there is no excessive roll occurs for the test conditions (3)–(5) because ϕ_max < 30°. According to the theory of fluid hydrodynamics, it has been validated that the separated vortices emitted from leeward side of the main hull and sail are symmetrical when incidence angle ranges from 5° to 30°. Thus, it may be known, the incidence angle is not the determinant factor of leading to excessive roll when the emergency rising is conducted at a certain speed with obtaining medium magnitude incidence angle. In contrast, before the model broaches the water surface, the maximum yaw deviation as well as maximum drift angle shows a significant positive correlation trend with the maximum roll angle as shown in Figure 6b. In other words, the strong coupling effect between roll and yaw motion in submarine emergency rising is verified from the perspective of model test.

Based on the above analysis of test results, it is further deduced that under the condition of moderate incidence angle (Θ_max < 30), the larger yaw deviation is the root cause of the excessive roll for the emergency rising at an initial speed. In order to validate this inference in detail, a numerical simulation that is equivalent to model test conditions is carried out and discussed in Section 3.2. Once verifying the validity of numerical simulation, another two emergency rising simulations under the condition of course keeping would be further conducted in Section 3.3 and Section 3.4 respectively, with making up for the deficiency of the model test.

3.2. The Validation Accuracy

In order to verify the validity and accuracy of the numerical method, a nonpropelled emergency rising simulation without any automatic steering devices is carried out first to reproduce the test results, as accurately as possible. The first working condition (1) in Table 1 is selected as the representative case of the numerical simulation due to its excessive yaw and roll simultaneously. For clarity, the initial conditions for simulation A are repeated as the following items:

• The blown ballast mass is 44.886 kg, and locates its mass centroid with 362 mm forward from midship and 15.34 mm upward from central axis of model.
• The model is launched at speed of 1.15 m/s from the depth of 8 m, while the sternplane angle is fixed as 0°.
• The opening states of flooding holes remain unchanged for all simulations, as does the metacentric height of 0.02 m.

Comparing the simulation results with the experimental data to obtain the contrast curves as shown in Figure 7a–f. Firstly, it should be pointed out that the data of model tests and numerical simulation are aligned to the fifth second. As can be noticed from Figure 7f, at the beginning of the test, the towing sled is gradually accelerated from 0 m/s to about 1.15 m/s. The time histories of ship speed are converged to the target value in about 5 s. Due to the influence of the speed control accuracy of the towing sled, the actual speed of the model fluctuates slightly around the target speed of 1.15 m/s during the uniform motion phase. As shown in the figure, after the model is released, the axial velocity curves of the numerical simulation coincide with the model test. It can be seen that the effect of the slight velocity deviation around releasing is negligible.

From the pitch curve of Figure 7a, it can be seen that the pitch curve of the model test has similar trends with that of numerical simulation. On one hand, both of them are difficult to maintain constant pitch state after reaching the maximum pitch angle. This is because the pitching moment provided by the buoyancy load is relatively small when the longitudinal position of blown ballast is close to the center of gravity of the ship, i.e., 362 mm. With the increase of roll angle, as shown in Figure 7c, the longitudinal hydrodynamic moment that trims bow down applied on the sail increases gradually, with further coupling the restoring moment provided by metacentric height, so that the pitch angle decreases after reaching its peak value. On the other hand, the deviation of the predicted peak value of pitch is lower than experimental results. This may be ascribed to the numerical dissipation of the turbulence model in solving the flow around the hull, especially the numerical solution during the rapid bow-lifting phase. The same phenomenon also exists in the prediction of incidence angle as illustrated in Figure 7b.

As the model emerges through the surface, the gravity center of submarine rises sharply while its displacement decreases as well, so that the model trim decreases rapidly and reaches its reverse maximum value of about 3°. At the same time, the roll angle increases drastically to 50°, resulting in the severe excessive roll. An obvious excessive rolling can be observed in the movie of the Simulation A, provided in the “Electronic supplementary material” (Video S1), where a significant yaw deviation is also shown under the condition of no rudder manipulating. More importantly, according to Figure 7c, the roll angle curves of numerical simulations and model tests also show a tendency of coincidence. Moreover, the peak-to-peak periods of rolling after surfacing are nearly identical, which indicates that both of them have the same metacentric height. In other words, it verifies the equivalence of static properties between the simulation model and the test model. To sum up, when the submarine blows off the ballast tank near its midship, it is prone to occurring excessive roll and threatening the safety of submarine and its crew.

Furthermore, according to the transverse and vertical comparison of the drift angle, roll angle and drift angle, the maximum underwater yaw angle and the drift angle corresponding to the twice roll instability states both are larger than 5°. At the same time, the peak values of the three curves show an obvious monotonous relationship, with the variation trend of parameters coinciding with the model test results. That is to say, the coupling effect between yaw and roll motion is verified by numerical simulation again. Last but not the least, the initial speed loss agrees well with the model test measurements as shown in Figure 7f, which demonstrates the completely equivalent and high reliability of numerical simulation for emergency rising of the nonpropelled model.

Overall, for the state parameters, e.g., pitch angle, incidence angle, axial velocity, etc., which presents good repeatability among multiple repeated tests, numerical simulation and model tests both have similar trends. As time goes on, the deviation of parameter curve is further reduced, even almost synchronized in time. Therefore, the prediction deviation of the longitudinal parameters in the bow-lifting phase is acceptable, with producing no essential impact on the rising stability analysis. Similarly, for the state parameters, e.g., roll angle, drift angle, yaw angle etc., which are susceptible to disturbance and exhibit dispersion within certain ranges of angle, the predicted results and test data also have similar trends. It is even possible to find out the test case that almost coincides with the numerical results, as shown in the above instance. Obviously, the roll, yaw, and drift angle are almost coincident. In other words, the numerical simulation has a high accuracy for predicting the lateral motion parameters, and it plays an extremely important role in analyzing the roll stability.

To sum up, according to the above analysis, it is demonstrated that the numerical simulation method established in this paper can substantially reproduce the model test results with high accuracy and reliability. Therefore, the numerical simulation can be used as an important auxiliary tool to study the emergency rising stability under the condition of rudder steering, which cannot be accomplished in model test up to now. Certainly, the self-propulsion simulation with rudder steering would provide theoretical basis and technical support for the further model test.

3.3. The Course Keeping Simulation Based on S Surface Controller

3.3.1. Description of S Surface Controller

The sigmoid surface functional controller (S surface controller) is a well-proven control algorithm with good effects. At present, it is widely used in various complex motion control of underwater vehicles [32,33,34]. The Sigmoid surface function is a nonlinear fitting to the output of the conventional fuzzy controller. When designing fuzzy controller, the control is usually loose on both two sides and dense in the middle, i.e., the control is relatively rough when the deviation is large, while the control is relatively fine when the deviation is small, which is consistent with the trend of sigmoid function. Therefore, the sigmoid function embodies the idea of fuzzy control [34]. In this paper, the model of the S plane controller for rudder steering is expressed as Equation (12):

$$\{\begin{array}{l}e=\psi -{\psi }_d\\ \dot{e}=(\psi -{\psi }_d)/\Delta t_i\\ o=\frac{2}{1+{\mathit{e}}^{-k_e\cdot e-k_v\cdot \dot{e}}}-1\\ {\delta }_{ri}={\delta }_{r\max }\cdot o\\ {\omega }_{ri}={\delta }_{ri}-{\delta }_{ri-1}/\Delta t_i\end{array}$$

(12)

In Equation (12), $\psi$ and ${\psi }_d$ are the current yaw angle and desired yaw angle. $e$ and $\dot{e}$ are the normalized yaw deviation and the rate of yaw deviation respectively. $k_e$ and $k_v$ are the control parameters corresponding to the deviation and the rate of deviation, which reflect the weight of the yaw angle deviation and the rate of yaw angle in the control. Both of them are parameters that need to adjust for the rudder angle controller. $o$ is the normalized output value of control system, i.e., normalized rudder angle, ${\delta }_{r\max }$ is the maximum allowable output rudder angle, ${\delta }_{ri}$ is the actual desired rudder angle, ${\omega }_{ri}$ is the rate of rudder steering acquired according to time differentiation of rudder angle in every time step. In numerical calculations, the rudder rotates to the target rudder angle according to the angular rate.

From the formula of the control model (12), the S surface controller combines the idea of fuzzy control with the structure of PD control, which not only simplifies the design of the controller, but also guarantee the control effect. There are only two control parameters in S surface controller model, i.e., $k_e$, $k_v$. These two parameters can be adjusted with referring to the concept of a proportional-differential (PD) controller. The parameters are manually adjusted so that the angular velocity of the rudder can meet the requirements of control system. Eventually, the yaw angle deviation can be limited within a small range to achieve good course keeping ability. $k_e$ and $k_v$ are larger, the response sensitivity to small deviations is higher, certainly, it is prone to causing the system oscillation as well. If the overshoot is large, it can be appropriate to reduce $k_e$ and increase $k_v$. Conversely, if the convergence is slow, it can be increased $k_e$ and decreased $k_v$ appropriately. In this paper, $k_e$ and $k_v$ are set to $k_e$ = 1, $k_v$ = 0.5 respectively.

Furthermore, in order to simulate the rudder steering characteristics of an actual submarine as truly as possible, an S surface controller based on conditional decision mechanism is proposed to steer the rudder. Specifically, when the deviation between the actual response rudder angle and the desired rudder angle is lower than 0.05°, the steering gear is not actuated, thereby avoiding the frequent manipulation of rudder in each time step. This ensures that the numerical simulation of rudder steering characteristics is more similar to the manipulation law of an actual submarine. The pseudo-code of decision mechanism is summarized as the framework shown in Figure 8.

In the above algorithm, ${\delta }_{ro}$ is the observed actual rudder angle, ${\delta }_{ri}$ is the desired rudder angle calculated by the controller at time step i. ${\omega }_{ri}$ is the angular velocity of rudder, and the index is set as an indicator to judge whether manipulating rudder or not when the rudder angle deviation locates within the range of 0.05°~1°. Because this rudder angle deviation will fall within the range of 0.05°~1° for twice, i.e., in the course of manipulating rudder and deviation accumulation without steering, respectively, which needs to be distinguished.

3.3.2. The Course Keeping Simulation

Based on the simulation A in Section 3.3.2, only a S surface controller is added into the simulation B through user defined function (UDF) to further simulate the emergency rising without propeller. The simulation results are compared with those of 3.2, as shown in Figure 9.

Firstly, in order to verify the effectiveness of the S-plane control method, the manipulation curves of the rudder are described in Figure 9f. It can be seen from the enlarged local curves that the rudder manipulation is performed only when the deviation larger than 0.5°, which is more in line with the actual steering mechanism, and further achieves the tracking control of the desired rudder angle. The yaw angle curve, as shown in Figure 9e, also indicates the effectiveness of the rudder steering. A local view of the rudder in the animation of Simulation B (see the “Electronic supplementary material”, Video S2) reveals the rotation of the rudder, where both roll and yaw angles are significantly reduced by manipulating rudder.

As shown in Figure 9a, the parametric comparison curves between two numerical simulation results (steering and no steering) and the experimental data (no steering) are presented. Firstly, compare the numerical results under two different conditions according to whether manipulating the rudder or not. As shown in Figure 9a–b, during the period of 5–15 s, the two predicted pitch curves as well as incidence curves substantially coincide, which shows that the numerical calculation has good computational stability. After t = 15 s, the pitch angle deviation between two numerical simulations increases gradually, which exactly reflects the effect of steering rudder on the trim. In other words, the manipulation of the rudder keeps the trim of the model before surfacing is relatively stable, which improves the longitudinal stability of the model to some extent. Obviously, this further leads to a nearly increase one time in the minimum pitch angle after the model surfacing, which is also understandable. In summary, the numerical simulation has high reliability for solving rudder steering process.

From the yaw angle curve of Figure 9c, it can be seen that the underwater yaw angle before surfacing is limited to 0.2°, which indicates that the model has excellent course keeping ability, and this further proves the effectiveness of the S surface controller. To be more precise, after t = 15 s, the yaw angle deviation of the model decreases gradually under the state of steering rudder, while the ones increases further when the rudder is not manipulated. Corresponding to this, according to curves of roll angle, the maximum roll angle in the rudder steering state is only 2°. In contrast, when the rudder is not manipulated, the maximum roll angle of the model approaches to 50° after surfacing. Undoubtedly, there occurs an obvious excessive roll during the underwater process of the emergency rising. In conclusion, for two numerical simulations with identical computational parameters and meshing, merely changing the steering state of the rudder, the roll stability of the model is significantly enhanced with reducing the maximum roll angle by 96%, which further validates the effectiveness of numerical simulation for the rudder manipulation.

Furthermore, by comparing the yaw angle and drift angle curves, it can be seen that under the steering rudder state, the underwater drift angle before surfacing is lower than 5°, and the variation trend of drift angle also present obvious linear correspondence with the yaw angle curve. Thereby it is demonstrated that the increase of drift angle is caused by the variation of yaw angle. In other words, the aim of improving the course keeping ability of the model is to enhance the roll stability through the way of limiting the growth of drift angle. This further proves that, from the perspective of reverse verification, excessive yaw deviation is the root cause of roll instability for emergency rising at medium incidence angle.

To sum up, in the small longitudinal blown position of 362 mm forward from the model center of gravity, the yaw angle can be controlled within the deviation range of 0.2° by manipulating the rudder, which is helpful to limit the growth of the model’s underwater drift angle, and then greatly reduce the maximum roll angle after surfacing. Under these conditions, the model can not only obtain good course keeping ability, but also avoid excessive pitch and excessive roll simultaneously. Consequently, the model will achieve the best emergency rising stability and safety.

3.4. The Self-Propulsion Simulation

In order to validate the feasibility of direct numerical simulation method for emergency rising simulation of self-propelled submarine, and further analyze the effect of steering rudder on the roll stability in self-propulsion simulation, in this section a self-designed propeller model is coupled with the main hull to conducts emergency rising simulation. The simulation results of the nonpropelled model with steering rudder are compared with the simulation results of the 2 different steering states in the self-propulsion simulation, and the curves shown in Figure 10a–h are obtained.

Compared with the conventional virtual disk method, this paper adopts a propeller model with real geometry, taking into account the interaction among the propeller, the main hull and sternplanes, which improves the accuracy of the simulation results. Figure 10g–h shows the parametric curves during the process of propeller revolution adjustment and initial velocity matching respectively, before the model is released. For obtaining the same initial advancing state as the non-propeller model test, at the beginning of the simulation, the model is set to move at the desired speed in the x-axis direction. Then the propeller thrust is balanced with the hull resistance by manually adjusting the propeller rpm for several times. According to Figure 10g, the propeller revolutions is gradually adjusted from 240 rpm to 210 rpm, so that the thrust and resistance achieve equilibrium when the initial velocity is maintained at 1.15 m/s. After the model is released, the increase of drag is attributed to the pressure difference between bow and stern caused by the model pitching. When the initial velocity is almost stable, the other five degrees of freedom motions are released to imitate the emergency rising, meanwhile, the propeller speed remains constant during the remaining simulation.

Firstly, according to the curves of axial velocity in Figure 10f, except for the peak value, the axial velocity of the self-propeller model is about 20% higher than that of the nonpropelled model during the whole process of emergency buoyancy. Obviously, it is ascribed to the propeller thrust. Furthermore, it can be seen from the enlarged velocity curves that, in the nonpropelled simulation, the velocity decreases slightly in the initial stage after release, while the velocity increases directly for the self-propulsion simulation. Similarly, due to the large axial velocity of the self-propelled model, the maximum pitch angle and incidence angle of the model are reduced by 23% and 8%, respectively, compared with the peak value of the nonpropelled model, as shown in Figure 10a,b. In addition, note that the influence of propeller mass and centroid on the center of gravity of the main hull is not considered in this paper. Therefore, it is acceptable that the model trims stern down in the stable advancing on the free surface, because this does not affect the discussion on course keeping ability and the qualitative analysis about the influence of the propeller.

Since the propeller is a left rotary propeller, the model always rolls to the starboard during the underwater rising process due to the propeller torque. For comparison, the roll and yaw angles of the nonpropelled simulation are taken as opposite values. It can be seen from roll angle curves in Figure 10c that, the maximum roll angle is only about 13° under the condition of steering rudder, which is significantly decreased by 70% than the roll angle in the simulation without steering rudder. Meanwhile, by comparing curves in Figure 10c for roll angle and Figure 10d for drift angle, we can see that they have similar trends. In terms of the above 3 various states, the roll angle increases with the rising trend of drift angle, i.e., the drift angle is a direct reflection of the roll angle. By manipulating the rudder, the drift angle is limited, which in turn weakens the roll angel. This shows that by keeping the course of the self-propelled model, the roll stability is significantly improved. Two animations of simulation C showing the real-time changes of attitude curves (see the “Electronic supplementary material”, Video S3 and S4), which intuitively illustrates the beneficial effect of manipulating rudder on reducing roll amplitudes in self-propulsion simulation.

On the other hand, the roll angle of self-propulsion simulation increases by about 11° compared with that of the nonpropelled simulation, which is mainly attributed to the effect of propeller torque. Correspondingly, the course keeping ability of self-propelled model is slightly deteriorated as illustrated in the Figure 10e. It is also mainly is attributed to the reaction influence of roll motion caused by the propeller torque. In other words, the coupling effect between roll and yaw motion makes the model yaw angle slightly increase around the time of broaching the free surface. However, the final heading is also kept in the initial direction, which verifies the effectiveness of the S surface controller in the self-propulsion simulation.

In summary, in the constant 362 mm blown position, for the self-propulsion simulation, when the rudder is not manipulated, the model is prone to excessive roll due to the influence of propeller torque. In the meantime, because of the influence of propeller torque and smaller pitch angle, the course keeping ability of the self-propelled model is slightly deteriorated. However, by manipulating the rudder, it is still possible to limit underwater yaw angle within a certain range, e.g., 1°, and further restrain the growth of drift angle, thereby significantly reducing the roll angle. In other words, it is reasonable to improve the roll stability of the model by enhancing the course keeping ability in the self-propulsion simulation. It is proven once again that under the condition of moderate angle of attack, the yaw instability is the root source of the emergency rising roll instability.

4. Conclusions

This paper presents a direct numerical simulation method based on RANS equation for submarine emergency rising emulation. Numerical simulations for 3 various conditions, i.e., 2 nonpropelled model simulations based on whether rudder steering or not, and self-propulsion simulation with rudder steering, are performed using Star-CCM+ solver. Here are some conclusions that can be drawn from the above discussion:

First, the predicted attitudes of the nonpropelled submarine without rudder steering agree well with the experimental results, which proves that the present CFD method coupled with overset grid technology is suitable for the direct numerical simulation of emergency rising. In addition, both numerical simulation and model tests indicate the strong coupling effect between roll and yaw motions, which serves as the theoretical foundation for the following studies.

Second, during the rudder steering simulations, the rudder rotational rate is controlled by a S surface controller based on conditional determination mechanism to achieve excellent course keeping capability. Consequently, the heading angle is limited within the range of 0.2°, while the maximum roll angle is reduced by 96%, which validates that the roll stability is significantly enhanced by improving the course keeping ability. More importantly, it is also proved that the excessive yaw deviation is the root cause of emergency rising roll instability at moderate incidence angles. Concisely, the model achieves optimal emergency rising stability under the condition of Simulation B.

Third, in the self-propulsion simulations, the time histories of propeller revolutions and initial velocity are converged to the target value in about 5 s, and then release other 5 DOF motions to implement emergency rising. Due to the influence of propeller torque, the course keeping capability of the self-propelled model is slightly deteriorated. However, by manipulating the rudder, it is still possible to limit underwater yaw angle within 0.8°, and further restrain the growth of drift angle, thereby significantly reducing the roll angle by 70%. Moreover, the coupling effect between roll and yaw is demonstrated through the submarine self-propulsion simulation again.

Future work will focus on self-propulsion emergency rising simulation with sternplanes manipulating as well as the full-scale model simulation. There are some difficulties in acquiring the metacentric height curves of the model surfacing process. Furthermore, more work will be done to analyze the influence of scale effect on the roll stability of submarine emergency rising.