Numerical Investigation of Internal Solitary Wave Forces on Submarines in Continuously Stratified Fluids

A numerical model of internal solitary waves in continuously stratified fluids is developed by introducing a density transport equation to the three-dimensional Navier–Stokes equation and adopting the fully nonlinear models of the Dubreil-Jacotin-Long equation to obtain the initial field of the ISW. The corresponding turbulence model has also been modified to ensure that it considers the variable density field. Comparisons between numerical simulation results and experimental results show that the total resistance, the nondimensional pressure coefficient, and the nondimensional friction coefficient for the standard submarine model proposed by the Defense Advanced Research Projects Agency under different flow field conditions are highly consistent with the experimental results. The model established is used to numerically analyse the forces and moments of the standard submarine model encountering ISWs at different submergence depths. The influence of the rotation centre position on the moment is discussed, and the position range of the optimal rotation centre is proposed.


Introduction
Internal solitary waves (ISWs) are a kind of fluctuation that occurs in stratified seawater and often appears in the actual ocean environment. ISWs carry much energy during propagation to induce strong convergence [1,2] and divergence of water bodies and sudden strong currents, which have an important impact on the ecological environment of the sea [3], the exploitation of oil and gas resources [4], and the safety of submarines [5,6]. When a submarine encounters ISWs, it will be dragged downwards and bumped violently. Submarines will be damaged because of rapid sinking exceeding the maximum safe submergence depth, as was the case with the KRI Nanggala-402 submarine disaster [7]. To evaluate the influence of ISWs on a submarine, hydrodynamic forces due to ISWs must be investigated.
At present, most investigations into the interaction between ISWs and a submarine are focused on the standard submarine model proposed by the Defense Advanced Research Projects Agency (DARPA), which is called the DARPA SUBOFF submarine model and can be represented by the bare hull (AFF-1) and fully appended hull (AFF-8) with stern appendages and a sail [8]. For example, Fu et al. [9] applied Fluent software to study the influence of ISWs on the bare hull and calculated the load when the submerged body encounters ISWs. Chen et al. [10] applied Fluent software to develop an ISW numerical tank and simulated the interaction between an ISW and self-propelled SUBOFF. Guan et al. [11] applied Gerris software to study the forces of the fully appended hull caused by ISWs. Huang et al. [12] simulated the free motion of SUBOFF under the action of ISWs by using Fluent software. However, most of the above studies used the two-layer fluid system to investigate the interaction between ISWs and SUBOFF, and the continuous stratification in fluid density was not considered.
In fact, according to the measured data from the Asian Seas International Acoustics Experiment (ASIAEX), a certain continuous stratification in density occurs along the water depth in the South China Sea, which is a prerequisite for the frequent occurrence of ISWs [13,14]. It is of great significance to investigate the interaction between ISWs and submarines in continuously stratified fluids to reflect the effect of the actual ocean environment. Therefore, a numerical model of ISWs in continuously stratified fluid will be developed and applied to investigate the forces of ISWs on a submarine.
The outline of the paper is described as follows. First, in Section 2, the development of a numerical model of ISWs in continuously stratified fluid is presented. Then, the verification of the model is given in Section 3. Subsequently, in Section 4, the forces of ISWs on SUBOFF are obtained through numerical simulation, and the effects of submergence depth are discussed. Finally, conclusions are drawn in Section 5.

Governing Equations
In order to reflect the characteristics of continuously stratified fluid in the ocean, we introduce the convection-diffusion equation of density to a fully 3D Navier-Stokes equation. The Boussinesq approximation is introduced because the density variation in the actual marine environment is small. The governing equations of the model are: where U = (u i , u j , u k ) is the velocity vector, ∇ is the gradient operator, t is the time, Q is the source term, ρ 0 is the reference density, ρ is the density field, p _rgh is a modified pressure field, g is the gravitational acceleration vector, and X is the position vector. v Eff is the effective kinematic viscosity, defined as v Eff = µ Eff/ ρ 0 , where µ Eff is the effective dynamic viscosity, including the molecular viscosity (µL) and turbulent viscosity (µ t ). D is the diffusion coefficient, and its value is the same as the effective dynamic viscosity (µ Eff ). Ω is the Coriolis parameter, which is twice the speed of rotation around the vertical unit vector e 3 = (0,0,1). r is the damping coefficient, and the specific formula [15] is as follows: where x s is the starting point of the sponge layer, x e is the endpoint of the sponge layer, the length of the sponge layer is twice the wavelength, a c and M are undetermined coefficients, and a c = 20 and M = 2 are recommended [15]. The turbulent viscosity in the above equation needs to be solved by the turbulence model. As an effective turbulence model, the two-equation k-ε model is widely used to solve turbulence. For example, Small and Hornby [16] use the k-ε model to close the equations when studying the propagation and evolution of ISWs. However, the twoequation k-ε model is generally suitable for high Reynolds number conditions, but for low Reynolds number conditions, it will cause excessive dissipation [17], and it cannot reflect the proper behaviour of the turbulent boundary layers caused by the adverse pressure gradients [18]. Huang et al. [19] and Menter et al. [20] suggested using the k-ω SST model to obtain substantially more accurate results. Considering the characteristics of ISWs, such as large amplitude, small wave speed, and small wave energy attenuation, this paper uses the k-ω SST model to solve the turbulence viscosity.
Notably, incompressible turbulence models in OpenFOAM are actually designed for strict incompressible flows with constant density, and they neglect the density variation in the flow field. For the numerical model of an ISW in continuously stratified fluid, the k-ω SST model must be modified to consider the variable density field. The modified specific equations can be expressed as: where k is the turbulent kinetic energy, ω is the specific dissipation rate, P k is the production term of k, given by P k = τ R : ∇U, P * k is related to the production term of turbulence kinetic energy P k in the k equation, v t is the turbulent kinematic viscosity, S t is the mean rate of the flow strain, given by S t = 0.5(∇U + ∇U T ), the model constants are assigned the values β * = 0.09, a 1 = 0.31, c 1 = 10 and C µ = 0.09, F 1 and F 2 are blending functions, and the values of σ k , σ ω , C γ and C β are blended using the equation Φ = F 1 Φ 1 + (1 − F 1 ) Φ 2 , in which Φ 1 and Φ 2 are given in Table 1.

Numerical Methods
The governing equations are numerically discretised based on the finite volume method. The finite volume method requires that Equations (2) and (3) are satisfied over the control volume around the centre point in integral form. In the discretisation process of the finite volume method, Gauss's theorem is applied to convert the volume integral into a surface integral, and then the obtained surface integral form is transformed into the algebraic form of the summation of fluxes over the faces of the control volume (Li et al., 2021). Equations (2) and (3) are solved by the PIMPLE algorithm [21], which combines the pressure implicit with splitting of operator (PISO) algorithm [22] and semi-implicit method for pressure-linked equations (SIMPLE) algorithm [23].

Initialised Field of ISW Generation
Fully nonlinear models of the DJL equation are applied to generate ISWs [24,25]. The DJL equation is as follows: where η is the isopycnal displacement, H is the water depth, c is the propagation speed, N is the definition of the buoyancy frequency, and z is the vertical position.
where ρ 0 (z) is the reference density, and g is the gravitational acceleration. η and c are obtained by solving the above DJL equation, and then the internal solitary wave-induced velocity field can be obtained through the relationship ψ = ηc, where ψ is the stream function. The DJLES open-source package developed by Dunphy et al. [26] is used to solve the DJL equations. Then, the initial velocity field calculated by DJLES is input into OpenFOAM for the start of the numerical simulation.

Model Verification
The capability of the numerical model for the simulation of ISWs was verified in Li et al. [27,28]. In this section, the simulated results of the flow resistance of a submarine model are compared with experimental measurements.

Computational Settings for SUBOFF Simulation
The standard submarine model proposed by the Defense Advanced Research Projects Agency (DARPA), called DARPA SUBOFF, will be used to investigate ISW-submarine interactions. Therefore, the SUBOFF resistance in the complex flow field is verified through the towing tank experiment of Liu and Huang [29]. The SUBOFF models of the bare hull (AFF-1) and fully appended hull (AFF-8) with stern appendages and a sail are shown in Figures 1 and 2, respectively. The main parameters of the SUBOFF model are shown in Table 2.
where ρ0(z) is the reference density, and g is the gravitational acceleration. η and c are obtained by solving the above DJL equation, and then the internal solitary wave-induced velocity field can be obtained through the relationship ψ = ηc, where ψ is the stream function. The DJLES open-source package developed by Dunphy et al. [26] is used to solve the DJL equations. Then, the initial velocity field calculated by DJLES is input into OpenFOAM for the start of the numerical simulation.

Model Verification
The capability of the numerical model for the simulation of ISWs was verified in Li et al. [27,28]. In this section, the simulated results of the flow resistance of a submarine model are compared with experimental measurements.

About the SUBOFF Model
The standard submarine model proposed by the Defense Advanced Research Projects Agency (DARPA), called DARPA SUBOFF, will be used to investigate ISW-submarine interactions. Therefore, the SUBOFF resistance in the complex flow field is verified through the towing tank experiment of Liu and Huang [29]. The SUBOFF models of the bare hull (AFF-1) and fully appended hull (AFF-8) with stern appendages and a sail are shown in Figures 1 and 2, respectively. The main parameters of the SUBOFF model are shown in Table 2.

Computational Domain and Grids
The computational domain and parameter settings are consistent with the experimental settings of Liu and Huang [29]. The grids around the SUBOFF should accurately depict the surface of the SUBOFF hull and describe the boundary layer. The boundary layer is generally divided into a viscous sublayer, buffer layer and log-law layer according  The computational domain and parameter settings are consistent with the experimental settings of Liu and Huang [29]. The grids around the SUBOFF should accurately depict the surface of the SUBOFF hull and describe the boundary layer. The boundary layer is generally divided into a viscous sublayer, buffer layer and log-law layer according to the distance along the wall normal. Among these layers, the flow of the log-law layer is in a fully developed turbulent state, and the velocity distribution is close to the logarithmic law. For the study of the interaction between ISWs and SUBOFF, we consider the full development of turbulence, so the first layer of grid control points should be located in the log-law region of the boundary layer. The boundary layer region division can be determined by the dimensionless value y + , and the expression is as follows: According to Ferziger et al. [30], 1 < y + < 60 corresponds to the log-law region. Fivelevel local mesh refinement is used to characterise the SUBOFF surface and ensure that the first layer of grid control points is located in the log-law region, and fifteen layers of grids on the SUBOFF surface are used as a transition, with an increasing ratio of 1.1. Figure 3 shows the y + distribution on the SUBOFF surface; y + is concentrated in the range of 25-55. In the SUBOFF surface refinement process, complex parts (such as the sail and stern appendages) are further refined using local mesh refinement, as shown in Figure 4.

Computational Domain and Grids
The computational domain and parameter settings are consistent with the experimental settings of Liu and Huang [29]. The grids around the SUBOFF should accurately depict the surface of the SUBOFF hull and describe the boundary layer. The boundary layer is generally divided into a viscous sublayer, buffer layer and log-law layer according to the distance along the wall normal. Among these layers, the flow of the log-law layer is in a fully developed turbulent state, and the velocity distribution is close to the logarithmic law. For the study of the interaction between ISWs and SUBOFF, we consider the full development of turbulence, so the first layer of grid control points should be located in the log-law region of the boundary layer. The boundary layer region division can be determined by the dimensionless value y + , and the expression is as follows: According to Ferziger et al. [30], 1 < y + < 60 corresponds to the log-law region. Fivelevel local mesh refinement is used to characterise the SUBOFF surface and ensure that the first layer of grid control points is located in the log-law region, and fifteen layers of grids on the SUBOFF surface are used as a transition, with an increasing ratio of 1.1. Figure  3 shows the y + distribution on the SUBOFF surface; y + is concentrated in the range of 25-55. In the SUBOFF surface refinement process, complex parts (such as the sail and stern appendages) are further refined using local mesh refinement, as shown in Figure 4.    Figure 5 shows the simulated and measured total resistance of the bare hull (SUB-OFF-AFF-1) and fully appended hull (SUBOFF-AFF-8). The numerical simulation results are consistent with the experimental results, and as the flow velocity increases, the total resistance of SUBOFF gradually increases. According to the relative errors in Tables 3 and  4, the maximum relative errors of resistance for SUBOFF-AFF-1 and SUBOFF-AFF-8 are 2.849% and 5.836%, respectively. Under the same flow velocity condition, the total resistance of SUBOFF-AFF-8 is larger than that of SUBOFF-AFF-1 because of the effect of the appendages, and the difference in total resistance becomes larger as the flow velocity increases.     Figure 5 shows the simulated and measured total resistance of the bare hull (SUBOFF-AFF-1) and fully appended hull (SUBOFF-AFF-8). The numerical simulation results are consistent with the experimental results, and as the flow velocity increases, the total resistance of SUBOFF gradually increases. According to the relative errors in Tables 3 and 4, the maximum relative errors of resistance for SUBOFF-AFF-1 and SUBOFF-AFF-8 are 2.849% and 5.836%, respectively. Under the same flow velocity condition, the total resistance of SUBOFF-AFF-8 is larger than that of SUBOFF-AFF-1 because of the effect of the appendages, and the difference in total resistance becomes larger as the flow velocity increases.   Figure 5 shows the simulated and measured total resistance of the bare hull (SUB-OFF-AFF-1) and fully appended hull (SUBOFF-AFF-8). The numerical simulation results are consistent with the experimental results, and as the flow velocity increases, the total resistance of SUBOFF gradually increases. According to the relative errors in Tables 3 and  4, the maximum relative errors of resistance for SUBOFF-AFF-1 and SUBOFF-AFF-8 are 2.849% and 5.836%, respectively. Under the same flow velocity condition, the total resistance of SUBOFF-AFF-8 is larger than that of SUBOFF-AFF-1 because of the effect of the appendages, and the difference in total resistance becomes larger as the flow velocity increases.     In addition to the total resistance, the nondimensional pressure coefficient C P = (P 0 − P ∞ )/ρU 2 ∞ and the nondimensional friction coefficient C f = τ w /ρU 2 ∞ (τ w is the wall shear stress) are important indicators for verifying the established model. Liu and Huang [29] measured the nondimensional pressure coefficient and the dimensionless friction coefficient of the surface of the SUBOFF submerged body under the condition that the Reynolds number is 1.4 × 10 7 and the speed V = 3.34 m/s. Figure 6 shows the simulated and measured nondimensional pressure coefficient distribution along SUBOFF-AFF-1, and the numerical simulation results of Yang and Lohner [31] and Sezen et al. [32] are also shown for comparison. Figure 6 shows that the numerical simulation results of this paper are in good agreement with the experimental results of Liu and Huang [29] and closer to the experimental results than those of Yang and Lohner [31] and Sezen et al. [32] when In addition to the total resistance, the nondimensional pressure coefficient

Comparisons between Simulated and Measured Results
wall shear stress) are important indicators for verifying the established model. Liu and Huang [29] measured the nondimensional pressure coefficient and the dimensionless friction coefficient of the surface of the SUBOFF submerged body under the condition that the Reynolds number is 1.4 × 10 7 and the speed V = 3.34 m/s. Figure 6 shows the simulated and measured nondimensional pressure coefficient distribution along SUBOFF-AFF-1, and the numerical simulation results of Yang and Lohner [31] and Sezen et al. [32] are also shown for comparison. Figure 6 shows that the numerical simulation results of this paper are in good agreement with the experimental results of Liu and Huang [29] and closer to the experimental results than those of Yang and Lohner [31] and Sezen et al. [32] when x/L ≤ 0.2. Figure 7 shows the simulated and measured nondimensional friction coefficient distribution along SUBOFF-AFF-1. For x/L ≤ 0.2, the three numerical models of the present paper, Yang and Lohner [31] and Sezen et al. [32], are basically consistent. The numerical simulation results of this paper generally agree with the experimental measurements and numerical results of Yang and Lohner [31] and Sezen et al. [32]. Since Liu and Huang [29] did not provide measured pressure data on the leading edge of the sail of SUBOFF-AFF-8, the numerical simulation results of Alin et al. [33] are introduced for comparison of the nondimensional pressure coefficient distribution along SUBOFF-AFF-8 in Figure 8. Comparing the nondimensional pressure coefficients of SUBOFF-AFF-8 and SUBOFF-AFF-1 (Figures 6 and 8) shows that the nondimensional  Figure 7 shows the simulated and measured nondimensional friction coefficient distribution along SUBOFF-AFF-1. For x/L ≤ 0.2, the three numerical models of the present paper, Yang and Lohner [31] and Sezen et al. [32], are basically consistent. The numerical simulation results of this paper generally agree with the experimental measurements and numerical results of Yang and Lohner [31] and Sezen et al. [32]. pressure coefficient has a sudden change affected by the appendages, especially on the front edge of the sail. The numerical results of this paper are consistent with the numerical simulations of Alin et al. [33] on the leading edge of the sail and they fit well with the experimental results of Liu and Huang [29] in other parts.  Overall, according to the analysis of DARPA SUBOFF's total resistance, nondimensional pressure coefficient and nondimensional friction coefficient, the numerical simulation results of the model developed in this paper are highly consistent with the experimental results. The model can be used for subsequent research on SUBOFF forces due to ISWs. Since Liu and Huang [29] did not provide measured pressure data on the leading edge of the sail of SUBOFF-AFF-8, the numerical simulation results of Alin et al. [33] are introduced for comparison of the nondimensional pressure coefficient distribution along SUBOFF-AFF-8 in Figure 8. Comparing the nondimensional pressure coefficients of SUBOFF-AFF-8 and SUBOFF-AFF-1 ( Figures 6 and 8) shows that the nondimensional pressure coefficient has a sudden change affected by the appendages, especially on the front edge of the sail. The numerical results of this paper are consistent with the numerical simulations of Alin et al. [33] on the leading edge of the sail and they fit well with the experimental results of Liu and Huang [29] in other parts. pressure coefficient has a sudden change affected by the appendages, especially on the front edge of the sail. The numerical results of this paper are consistent with the numerical simulations of Alin et al. [33] on the leading edge of the sail and they fit well with the experimental results of Liu and Huang [29] in other parts.  Overall, according to the analysis of DARPA SUBOFF's total resistance, nondimensional pressure coefficient and nondimensional friction coefficient, the numerical simulation results of the model developed in this paper are highly consistent with the experimental results. The model can be used for subsequent research on SUBOFF forces due to ISWs. Overall, according to the analysis of DARPA SUBOFF's total resistance, nondimensional pressure coefficient and nondimensional friction coefficient, the numerical simulation results of the model developed in this paper are highly consistent with the experimental results. The model can be used for subsequent research on SUBOFF forces due to ISWs.

Numerical Analysis of ISW Forces on SUBOFF
The submarine may be dragged down by the ISWs and may even crash due to exceeding the maximum submergence depth. Therefore, this section mainly focuses on the forces of SUBOFF subjected to ISWs under different submerged conditions. The nondimensional force and moment are analysed, and the specific nondimensional form is defined as follows: where C T is the nondimensional drag force, C L is the nondimensional lift force, and C M is the nondimensional moment.

Computational Domain and Simulation Conditions
The SUBOFF is magnified 20-fold as a submarine prototype of actual scale. Its total length L s is 87.12 m, and the maximum radius of gyration is 10.16 m. The actual submarine is complicated in internal structure, and it cannot be considered a homogeneous body to determine the centre of rotation. For the SUBOFF numerical simulation, the centre of rotation is taken at the centreline of the axisymmetric component of SUBOFF and is a distance of L s /3 from the bow. The direction of the moment of SUBOFF is determined according to the right-hand rule, that is, the clockwise direction perpendicular to the paper is positive, and the anticlockwise direction is negative. The choice of the computational domain size of the ISW numerical tank directly affects the simulation time and the accuracy of the calculation results. According to the observation data of the Asian Seas International Acoustics Experiment (ASIAEX) [13,14] and considering the scale of ISWs and the influence of the boundary on the SUBOFF, the length, width and water depth of the numerical tank are taken as 8 km, 1 km, and 300 m, respectively. A schematic top view of the computational domain is shown in Figure 9.

Numerical Analysis of ISW Forces on SUBOFF
The submarine may be dragged down by the ISWs and may even crash due to exceeding the maximum submergence depth. Therefore, this section mainly focuses on the forces of SUBOFF subjected to ISWs under different submerged conditions. The nondimensional force and moment are analysed, and the specific nondimensional form is defined as follows: where CT is the nondimensional drag force, CL is the nondimensional lift force, and CM is the nondimensional moment.

Computational Domain and Simulation Conditions
The SUBOFF is magnified 20-fold as a submarine prototype of actual scale. Its total length Ls is 87.12 m, and the maximum radius of gyration is 10.16 m. The actual submarine is complicated in internal structure, and it cannot be considered a homogeneous body to determine the centre of rotation. For the SUBOFF numerical simulation, the centre of rotation is taken at the centreline of the axisymmetric component of SUBOFF and is a distance of Ls/3 from the bow. The direction of the moment of SUBOFF is determined according to the right-hand rule, that is, the clockwise direction perpendicular to the paper is positive, and the anticlockwise direction is negative. The choice of the computational domain size of the ISW numerical tank directly affects the simulation time and the accuracy of the calculation results. According to the observation data of the Asian Seas International Acoustics Experiment (ASIAEX) [13,14] and considering the scale of ISWs and the influence of the boundary on the SUBOFF, the length, width and water depth of the numerical tank are taken as 8 km, 1 km, and 300 m, respectively. A schematic top view of the computational domain is shown in Figure 9. The depths of the upper (h1) and lower (h2) layers are 50 m and 250 m, respectively, the densities of the upper and lower layers are 1022 kg/m 3 and 1028 kg/m 3 , respectively, the location of the centre of the pycnocline (zpyc) is −50 m, the pycnocline thickness (dpyc) is 20 m vertically, and the initial ISW amplitude (a) is 40 m. Slip boundary conditions are applied to the bottom and both sides, while cyclic boundary conditions are assigned to the inlet and outlet boundaries. The top boundary is a rigid lid. The boundary conditions related to the density field are no-flux boundary conditions.  To analyse the influence of the submergence depths on the forces, we set up four numerical simulation conditions, as shown in Figure 10. Cases 1~4 represent four conditions of the SUBOFF centreline, which are located at the upper edge of the ISW (z = −50 m), at the middle of the ISW (z = −70 m), at the lower edge of the ISW trough (z = −90 m), and 30 m lower than the trough bottom of the ISW (z = −120 m).

Drag Force
The simulated time series of the nondimensional drag forces of SUBOFF-AFF-1 and SUBOFF-AFF-8 are shown in Figure 11. The curves of the nondimensional drag forces for SUBOFF-AFF-1 and SUBOFF-AFF-8 are very similar under the same conditions. To analyse drag forces, the ISW-induced dynamic pressure field is given in Figure 12. Figure 12 shows that the ISW-induced dynamic pressure is largest in the centre of the ISW trough and is positive in the upper layer and negative in the lower layer.

Drag Force
The simulated time series of the nondimensional drag forces of SUBOFF-AFF-1 and SUBOFF-AFF-8 are shown in Figure 11. The curves of the nondimensional drag forces for SUBOFF-AFF-1 and SUBOFF-AFF-8 are very similar under the same conditions. To analyse drag forces, the ISW-induced dynamic pressure field is given in Figure 12. Figure 12 shows that the ISW-induced dynamic pressure is largest in the centre of the ISW trough and is positive in the upper layer and negative in the lower layer.

Drag Force
The simulated time series of the nondimensional drag forces of SUBOFF-AFF-1 and SUBOFF-AFF-8 are shown in Figure 11. The curves of the nondimensional drag forces for SUBOFF-AFF-1 and SUBOFF-AFF-8 are very similar under the same conditions. To analyse drag forces, the ISW-induced dynamic pressure field is given in Figure 12. Figure 12 shows that the ISW-induced dynamic pressure is largest in the centre of the ISW trough and is positive in the upper layer and negative in the lower layer.   In Case 1, the nondimensional drag force variation undergoes three stages with ISW propagation. The first stage increases from zero to the maximum drag force, the second stage decreases from the maximum drag force to the minus maximum drag force (reverse direction of ISW propagation), and the third stage increases from the minus maximum drag force to approximately zero. The variation process of the drag force depends on the relative position of the SUBOFF in the ISW. The fluid velocities of the ISW in the upper layer and lower layer are opposite to the centre of the pycnocline. The velocity in the upper layer has the same direction as ISW propagation, while that in the lower layer is in the opposite direction (marked with a negative value). In the first stage, as the ISW propagates to the SUBOFF, the centre of the pycnocline moves downwards, causing the SUBOFF to gradually enter the upper layer with increasing fluid velocity. The corresponding ISWinduced dynamic pressure is also greater, which leads to a continuous increase in the positive horizontal force. When the SUBOFF is completely above the pycnocline, the dynamic pressure difference between the bow and stern is the largest in the horizontal direction, and the nondimensional drag force reaches the maximum, as shown in Figure 13a. In the second stage, the SUBOFF is in the upper part of the trough area of the ISW. With the propagation of the ISW, the ISW-induced dynamic pressure difference on SUBOFF in the horizontal direction gradually decreases to zero and then gradually increases in the reverse direction, as shown in Figure 12. The trend of the force variation is the same as that of the ISW-induced dynamic pressure difference in the horizontal direction. When SUB-OFF encounters the pycnocline, the horizontal force reaches the negative maximum value (as shown in Figure 13b). In the third stage, as the ISW propagates, the centreline of the pycnocline moves upwards, causing the SUBOFF to gradually return to the initial state, the dynamic pressure difference gradually decreases to zero, and the value of the negative nondimensional drag force gradually decreases and returns to the initial zero value. In Case 1, the nondimensional drag force variation undergoes three stages with ISW propagation. The first stage increases from zero to the maximum drag force, the second stage decreases from the maximum drag force to the minus maximum drag force (reverse direction of ISW propagation), and the third stage increases from the minus maximum drag force to approximately zero. The variation process of the drag force depends on the relative position of the SUBOFF in the ISW. The fluid velocities of the ISW in the upper layer and lower layer are opposite to the centre of the pycnocline. The velocity in the upper layer has the same direction as ISW propagation, while that in the lower layer is in the opposite direction (marked with a negative value). In the first stage, as the ISW propagates to the SUBOFF, the centre of the pycnocline moves downwards, causing the SUBOFF to gradually enter the upper layer with increasing fluid velocity. The corresponding ISW-induced dynamic pressure is also greater, which leads to a continuous increase in the positive horizontal force. When the SUBOFF is completely above the pycnocline, the dynamic pressure difference between the bow and stern is the largest in the horizontal direction, and the nondimensional drag force reaches the maximum, as shown in Figure 13a. In the second stage, the SUBOFF is in the upper part of the trough area of the ISW. With the propagation of the ISW, the ISW-induced dynamic pressure difference on SUBOFF in the horizontal direction gradually decreases to zero and then gradually increases in the reverse direction, as shown in Figure 12. The trend of the force variation is the same as that of the ISW-induced dynamic pressure difference in the horizontal direction. When SUBOFF encounters the pycnocline, the horizontal force reaches the negative maximum value (as shown in Figure 13b). In the third stage, as the ISW propagates, the centreline of the pycnocline moves upwards, causing the SUBOFF to gradually return to the initial state, the dynamic pressure difference gradually decreases to zero, and the value of the negative nondimensional drag force gradually decreases and returns to the initial zero value.
In Case 2, with the arrival of the ISW, the lower fluid with negative velocity causes a leftward force on the SUBOFF. As the pycnocline of the ISW propagates to the SUBOFF and gradually passes through, the upper fluid with positive velocity produces a rightward force on the SUBOFF, making the nondimensional drag force gradually decrease to zero and then gradually increase in the positive direction. When SUBOFF just passes through the pycnocline and completely enters the upper fluid, the nondimensional drag force reaches its maximum value. The SUBOFF located in the trough of the ISW is subject to ISW-induced dynamic pressure field changes, and the nondimensional drag force is continuously reduced and then gradually increases in the reverse direction. This process is similar to the second stage of Case 1. When the bow of SUBOFF is about to pass through the centreline of the pycnocline, the nondimensional drag force reaches the reverse maximum value. As the SUBOFF traverses the ISW again, the dynamic pressure difference between the bow and stern gradually decreases to zero and then gradually increases in the positive direction. The law of force is the same as that of the ISW-induced dynamic pressure difference on SUBOFF in the horizontal direction. With the propagation of the ISW, the value of the nondimensional drag force gradually decreases and returns to the initial zero value. of the ISW-induced dynamic pressure difference in the horizontal direction. When SUB-OFF encounters the pycnocline, the horizontal force reaches the negative maximum value (as shown in Figure 13b). In the third stage, as the ISW propagates, the centreline of the pycnocline moves upwards, causing the SUBOFF to gradually return to the initial state the dynamic pressure difference gradually decreases to zero, and the value of the negative nondimensional drag force gradually decreases and returns to the initial zero value. In Case 2, with the arrival of the ISW, the lower fluid with negative velocity causes a leftward force on the SUBOFF. As the pycnocline of the ISW propagates to the SUBOFF and gradually passes through, the upper fluid with positive velocity produces a rightward force on the SUBOFF, making the nondimensional drag force gradually decrease to zero and then gradually increase in the positive direction. When SUBOFF just passes through the pycnocline and completely enters the upper fluid, the nondimensional drag force reaches its maximum value. The SUBOFF located in the trough of the ISW is subject to ISW-induced dynamic pressure field changes, and the nondimensional drag force is continuously reduced and then gradually increases in the reverse direction. This process is similar to the second stage of Case 1. When the bow of SUBOFF is about to pass through the centreline of the pycnocline, the nondimensional drag force reaches the reverse maximum value. As the SUBOFF traverses the ISW again, the dynamic pressure difference between the bow and stern gradually decreases to zero and then gradually increases in the positive direction. The law of force is the same as that of the ISW-induced dynamic pressure difference on SUBOFF in the horizontal direction. With the propagation of the ISW, the value of the nondimensional drag force gradually decreases and returns to the initial zero value.
In Case 3, the SUBOFF is located in the lower fluid of the ISW, and its centreline is located at the lower waveform edge of the ISW (z = −90 m). Under this condition, the SUBOFF does not cross the ISW, but it crosses the pycnocline below the waveform. As the ISW propagates, it not only causes a dynamic pressure difference between the left and right sides of the waveform, but also causes a dynamic pressure difference in the upper and lower fluids along the horizontal direction (as shown in Figure 12). With the arrival of the ISW, the lower fluid with negative dynamic pressure exerts a leftward force on the SUBOFF. When SUBOFF encounters the centreline of the pycnocline, the horizontal force reaches the reverse maximum value. After that, because of the change in the dynamic pressure difference on the SUBOFF, the nondimensional drag force gradually decreases to zero and then continuously increases in the positive direction. When SUBOFF starts to leave the pycnocline, the nondimensional drag force reaches its positive maximum value As SUBOFF gradually leaves the pycnocline, the nondimensional drag force continues to decrease and finally returns to zero.
In Case 4, the SUBOFF traverses neither the ISW waveform nor the pycnocline below In Case 3, the SUBOFF is located in the lower fluid of the ISW, and its centreline is located at the lower waveform edge of the ISW (z = −90 m). Under this condition, the SUBOFF does not cross the ISW, but it crosses the pycnocline below the waveform. As the ISW propagates, it not only causes a dynamic pressure difference between the left and right sides of the waveform, but also causes a dynamic pressure difference in the upper and lower fluids along the horizontal direction (as shown in Figure 12). With the arrival of the ISW, the lower fluid with negative dynamic pressure exerts a leftward force on the SUBOFF. When SUBOFF encounters the centreline of the pycnocline, the horizontal force reaches the reverse maximum value. After that, because of the change in the dynamic pressure difference on the SUBOFF, the nondimensional drag force gradually decreases to zero and then continuously increases in the positive direction. When SUBOFF starts to leave the pycnocline, the nondimensional drag force reaches its positive maximum value. As SUBOFF gradually leaves the pycnocline, the nondimensional drag force continues to decrease and finally returns to zero.
In Case 4, the SUBOFF traverses neither the ISW waveform nor the pycnocline below the wave surface. The nondimensional drag force is mainly affected by the dynamic pressure of the lower fluid of the ISW, and the overall change trend is similar to that of Case 3. Since SUBOFF is not affected by the pycnocline, the peak value of the nondimensional drag in Case 4 is larger than that in Case 3, and the peak value of the positive nondimensional drag is smaller.

Lift Force
The time series of the nondimensional lift force of the SUBOFF is shown in Figure 14. The nondimensional lift forces of SUBOFF-AFF-1 and SUBOFF-AFF-8 under the same conditions are similar. The fully appended hull adds stern appendages, which makes its horizontal projection area slightly larger, resulting in its nondimensional lift force being slightly larger than that of the bare hull. The graphs of the nondimensional lift force against time for Case 1, Case 2 and Case 3 are similar, and all the lift forces first increase in the negative direction and then gradually return to near zero along the positive direction. Figure 15 shows the dynamic pressure field around SUBOFF when the nondimensional lift force reaches the minus maximum. The pressure is greater on the upper side than on the lower side, which is most obvious in the stern of the SUBOFF. When SUBOFF crosses the ISW, the amplitude of the nondimensional lift force of Case 2 is the largest (as shown in Figure 15) and is approximately twice that of Case 1 and Case 3 and approximately tenfold that of Case 4. than on the lower side, which is most obvious in the stern of the SUBOFF. When SUBOFF crosses the ISW, the amplitude of the nondimensional lift force of Case 2 is the largest (as shown in Figure 15) and is approximately twice that of Case 1 and Case 3 and approximately tenfold that of Case 4.
As the ISWs propagate, the isobaric surface in the lower fluid will bend and deform. In Case 4, SUBOFF crosses the bent isobaric surface in the lower fluid, causing the nondimensional lift force of Case 4 to differ from that of the other three cases; that is, the lift force first shows a downwards trend, then increases to the maximum upwards force, then decreases to a negative value, and finally increases back to zero. The simulation results are basically consistent with the results found by Fu et al. [9] using Fluent software.  Comparing Figures 11 and 14 shows that the nondimensional lift force is one order of magnitude larger than the nondimensional drag force; that is, the vertical force is one order of magnitude larger than the horizontal force. For the cases of the ISW acting on SUBOFF, the nondimensional drag of Case 2 is the smallest, but the magnitude of the nondimensional lift is the largest, and there is a danger of being dragged into the deep sea and causing overturning and crushing. When the submergence depth of SUBOFF is lower than the ISW waveform (Case 4), the nondimensional drag and the nondimensional lift are at the same order of magnitude.  than on the lower side, which is most obvious in the stern of the SUBOFF. When SUBOFF crosses the ISW, the amplitude of the nondimensional lift force of Case 2 is the largest (as shown in Figure 15) and is approximately twice that of Case 1 and Case 3 and approximately tenfold that of Case 4.
As the ISWs propagate, the isobaric surface in the lower fluid will bend and deform. In Case 4, SUBOFF crosses the bent isobaric surface in the lower fluid, causing the nondimensional lift force of Case 4 to differ from that of the other three cases; that is, the lift force first shows a downwards trend, then increases to the maximum upwards force, then decreases to a negative value, and finally increases back to zero. The simulation results are basically consistent with the results found by Fu et al. [9] using Fluent software.  Comparing Figures 11 and 14 shows that the nondimensional lift force is one order of magnitude larger than the nondimensional drag force; that is, the vertical force is one order of magnitude larger than the horizontal force. For the cases of the ISW acting on SUBOFF, the nondimensional drag of Case 2 is the smallest, but the magnitude of the nondimensional lift is the largest, and there is a danger of being dragged into the deep sea and causing overturning and crushing. When the submergence depth of SUBOFF is lower than the ISW waveform (Case 4), the nondimensional drag and the nondimensional lift are at the same order of magnitude.  As the ISWs propagate, the isobaric surface in the lower fluid will bend and deform. In Case 4, SUBOFF crosses the bent isobaric surface in the lower fluid, causing the nondimensional lift force of Case 4 to differ from that of the other three cases; that is, the lift force first shows a downwards trend, then increases to the maximum upwards force, then decreases to a negative value, and finally increases back to zero. The simulation results are basically consistent with the results found by Fu et al. [9] using Fluent software.
Comparing Figures 11 and 14 shows that the nondimensional lift force is one order of magnitude larger than the nondimensional drag force; that is, the vertical force is one order of magnitude larger than the horizontal force. For the cases of the ISW acting on SUBOFF, the nondimensional drag of Case 2 is the smallest, but the magnitude of the nondimensional lift is the largest, and there is a danger of being dragged into the deep sea and causing overturning and crushing. When the submergence depth of SUBOFF is lower than the ISW waveform (Case 4), the nondimensional drag and the nondimensional lift are at the same order of magnitude. Figure 16 shows the time series of the nondimensional moment of the SUBOFF under different conditions. Obviously, the nondimensional moments of SUBOFF-AFF-1 and SUBOFF-AFF-8 have very similar variation processes under the same conditions. The nondimensional moment processes of Case 1, Case 2 and Case 3 show similar trends of first increasing and then decreasing, which causes the SUBOFF to rotate clockwise. Figure 17 shows the nondimensional force vector on the surface of the SUBOFF under Case 2 when the nondimensional moment reaches the maximum. Figure 17 shows that the nondimensional force on the upper side of the SUBOFF produces a substantially larger moment than that on the lower side and makes the SUBOFF rotate clockwise. In Case 2, when the SUBOFF is about to pass through the centreline of the pycnocline for the second time, the nondimensional moment reaches the maximum, which is 1.2 × 10 −2 and twice that of Case 3, or fourfold that of Case 4 (as shown in Figure 16). The nondimensional moment of Case 4 first increases in the positive direction, then continuously decreases to zero, then increases in the negative direction to the maximum value, and finally decreases to zero. Figure 16 shows the time series of the nondimensional moment of the SUBOFF under different conditions. Obviously, the nondimensional moments of SUBOFF-AFF-1 and SUBOFF-AFF-8 have very similar variation processes under the same conditions. The nondimensional moment processes of Case 1, Case 2 and Case 3 show similar trends of first increasing and then decreasing, which causes the SUBOFF to rotate clockwise. Figure  17 shows the nondimensional force vector on the surface of the SUBOFF under Case 2 when the nondimensional moment reaches the maximum. Figure 17 shows that the nondimensional force on the upper side of the SUBOFF produces a substantially larger moment than that on the lower side and makes the SUBOFF rotate clockwise. In Case 2, when the SUBOFF is about to pass through the centreline of the pycnocline for the second time, the nondimensional moment reaches the maximum, which is 1.2 × 10 −2 and twice that of Case 3, or fourfold that of Case 4 (as shown in Figure 16). The nondimensional moment of Case 4 first increases in the positive direction, then continuously decreases to zero, then increases in the negative direction to the maximum value, and finally decreases to zero.

Effects of the Rotation Centre on Momentum
The position of the centre of rotation is very important to the magnitude and direction of the moment received by the submarine. The effects of the rotation centre position on the moment of SUBOFF under the action of the ISW are analysed in this section, and the optimal rotation centre position range is further obtained.
Since the main body of the SUBOFF is axisymmetric, the centre of rotation can be assumed to be on the revolving axis of the axisymmetric component. The position of the bow of the SUBOFF is specified as the zero point, and the position of the centre on the rotation axis is represented by S = x/Ls (0 ≤ S ≤ 1), where x is the length from the bow, and   Figure 16 shows the time series of the nondimensional moment of the SUBOFF under different conditions. Obviously, the nondimensional moments of SUBOFF-AFF-1 and SUBOFF-AFF-8 have very similar variation processes under the same conditions. The nondimensional moment processes of Case 1, Case 2 and Case 3 show similar trends of first increasing and then decreasing, which causes the SUBOFF to rotate clockwise. Figure  17 shows the nondimensional force vector on the surface of the SUBOFF under Case 2 when the nondimensional moment reaches the maximum. Figure 17 shows that the nondimensional force on the upper side of the SUBOFF produces a substantially larger moment than that on the lower side and makes the SUBOFF rotate clockwise. In Case 2, when the SUBOFF is about to pass through the centreline of the pycnocline for the second time, the nondimensional moment reaches the maximum, which is 1.2 × 10 −2 and twice that of Case 3, or fourfold that of Case 4 (as shown in Figure 16). The nondimensional moment of Case 4 first increases in the positive direction, then continuously decreases to zero, then increases in the negative direction to the maximum value, and finally decreases to zero.

Effects of the Rotation Centre on Momentum
The position of the centre of rotation is very important to the magnitude and direction of the moment received by the submarine. The effects of the rotation centre position on the moment of SUBOFF under the action of the ISW are analysed in this section, and the optimal rotation centre position range is further obtained.
Since the main body of the SUBOFF is axisymmetric, the centre of rotation can be assumed to be on the revolving axis of the axisymmetric component. The position of the bow of the SUBOFF is specified as the zero point, and the position of the centre on the rotation axis is represented by S = x/Ls (0 ≤ S ≤ 1), where x is the length from the bow, and

Effects of the Rotation Centre on Momentum
The position of the centre of rotation is very important to the magnitude and direction of the moment received by the submarine. The effects of the rotation centre position on the moment of SUBOFF under the action of the ISW are analysed in this section, and the optimal rotation centre position range is further obtained.
Since the main body of the SUBOFF is axisymmetric, the centre of rotation can be assumed to be on the revolving axis of the axisymmetric component. The position of the bow of the SUBOFF is specified as the zero point, and the position of the centre on the rotation axis is represented by S = x/L s (0 ≤ S ≤ 1), where x is the length from the bow, and L s is the total length of the SUBOFF. According to the shape characteristics of the bare hull (SUBOFF-AFF-1), the centre of rotation of the submarine should be between 1/2 of the submarine and the bow. Therefore, for the rotation centre position of SUBOFF-AFF-1, we take a series of positions within the interval of 1/3 ≤ S ≤ 1/2 for numerical simulation, as shown in Figure 18. As the fully appended hull (SUBOFF-AFF-8) adds stern appendages and a sail, the centre of rotation should be in the middle of the submarine. We take a series of positions within the interval of 2/6 ≤ S ≤ 4/6 for numerical simulation research, as shown in Figure 18. The specific values of the rotation centre positions of the two types are shown in Table 5.
submarine and the bow. Therefore, for the rotation centre position of SUBOFF-AFF-1, we take a series of positions within the interval of 1/3 ≤ S ≤ 1/2 for numerical simulation, as shown in Figure 18. As the fully appended hull (SUBOFF-AFF-8) adds stern appendages and a sail, the centre of rotation should be in the middle of the submarine. We take a series of positions within the interval of 2/6 ≤ S ≤ 4/6 for numerical simulation research, as shown in Figure 18. The specific values of the rotation centre positions of the two types are shown in Table 5.   Figure 19 shows the nondimensional moment processes of different rotation centre positions for SUBOFF-AFF-1 under four conditions (Case 1, Case 2, Case 3 and Case 4). For each condition, 16 kinds of rotation centre positions of SUBOFF-AFF-1 were simulated, while only eight curves are shown in the figure for clear vision. In Case 1, Case 2, and Case 3, the nondimensional moment obtains the positive maximum value (in clockwise) at S = 1/3 and the negative maximum value (in anticlockwise) at S =1/2. For SUBOFF-AFF-1, with the continuous increase in the distance from the rotation centre position to the bow (1/3 to 1/2), the nondimensional moment with a positive value continues to decrease, gradually appears in the form that first rotates anticlockwise and then rotates clockwise, and finally develops into a form in which the negative nondimensional moment dominates. In Case 4, the nondimensional moment curve is different from the other three cases. The appearance of the nondimensional moment is gradually delayed, and all the nondimensional moments first rotate clockwise and then anticlockwise. When the position of the rotation centre is the same among the four cases, the magnitude of the nondimensional moment is the largest in Case 2 and smallest in Case 4. Generally, for SUB-OFF-AFF-1 crossing an ISW, the nondimensional moment is the largest, which is the most dangerous condition.   Figure 19 shows the nondimensional moment processes of different rotation centre positions for SUBOFF-AFF-1 under four conditions (Case 1, Case 2, Case 3 and Case 4). For each condition, 16 kinds of rotation centre positions of SUBOFF-AFF-1 were simulated, while only eight curves are shown in the figure for clear vision. In Case 1, Case 2, and Case 3, the nondimensional moment obtains the positive maximum value (in clockwise) at S = 1/3 and the negative maximum value (in anticlockwise) at S = 1/2. For SUBOFF-AFF-1, with the continuous increase in the distance from the rotation centre position to the bow (1/3 to 1/2), the nondimensional moment with a positive value continues to decrease, gradually appears in the form that first rotates anticlockwise and then rotates clockwise, and finally develops into a form in which the negative nondimensional moment dominates. In Case 4, the nondimensional moment curve is different from the other three cases. The appearance of the nondimensional moment is gradually delayed, and all the nondimensional moments first rotate clockwise and then anticlockwise. When the position of the rotation centre is the same among the four cases, the magnitude of the nondimensional moment is the largest in Case 2 and smallest in Case 4. Generally, for SUBOFF-AFF-1 crossing an ISW, the nondimensional moment is the largest, which is the most dangerous condition.
The  Figure 20 shows that in Case 1, Case 2 and Case 3, the maximum value of the positive nondimensional moment (clockwise rotation) appears at S = 1/3, and the maximum value of the negative nondimensional moment (anticlockwise rotation) appears at S = 4/6. With the continuous increase in the position of the SUBOFF-AFF-8 rotation centre (1/3 to 2/3), the variation in the nondimensional moment is similar to that of SUBOFF-AFF-1. When the position of the rotation centre is the same, among the four cases, the magnitude of the nondimensional moment is the largest in Case 2 and smallest in Case 4, which is the same as SUBOFF-AFF-1. Comparing the nondimensional moments of SUBOFF-AFF-1 and SUBOFF-AFF-8 (Figures 19 and 20) at the same value of S, for SUBOFF within the wave height range of the ISW, the nondimensional moment of SUBOFF-AFF-1 is slightly larger than that of SUBOFF-AFF-8; for SUBOFF lower than the trough bottom of the ISW, the nondimensional moment of SUBOFF-AFF-1 is slightly smaller than that of SUBOFF-AFF-8. To compare the difference between the nondimensional moments of SUBOFF-AFF-1 and SUBOFF-AFF-8 more clearly, the nondimensional moment at S = 0.466 is given (as shown in Figure 21).  Figure 20 shows that in Case 1, Case 2 and Case 3, the maximum value of the positive nondimensional moment (clockwise rotation) appears at S = 1/3, and the maximum value of the negative nondimensional moment (anticlockwise rotation) appears at S = 4/6. With the continuous increase in the position of the SUBOFF-AFF-8 rotation centre (1/3 to 2/3), the variation in the nondimensional moment is similar to that of SUBOFF-AFF-1. When the position of the rotation centre is the same, among the four cases, the magnitude of the nondimensional moment is the largest in Case 2 and smallest in Case 4, which is the same as SUBOFF-AFF-1. Comparing the nondimensional moments of SUB-OFF-AFF-1 and SUBOFF-AFF-8 (Figures 19 and 20) at the same value of S, for SUBOFF within the wave height range of the ISW, the nondimensional moment of SUBOFF-AFF-1 is slightly larger than that of SUBOFF-AFF-8; for SUBOFF lower than the trough bottom of the ISW, the nondimensional moment of SUBOFF-AFF-1 is slightly smaller than that of SUBOFF-AFF-8. To compare the difference between the nondimensional moments of SUBOFF-AFF-1 and SUBOFF-AFF-8 more clearly, the nondimensional moment at S = 0.466 is given (as shown in Figure 21).
For SUBOFF-AFF-1 and SUBOFF-AFF-8 crossing the ISW, the nondimensional moment is the largest, which is the most dangerous condition.

Analysis of the Optimal Rotation Centre Position
According to the results of Section 4.3, the rotation centre position causes drastic variation in the nondimensional moment when the SUBOFF (AFF-1 and AFF-8) encounters the ISW. Therefore, it is meaningful to quantitatively explore the optimal range of the rotation centre position for the design of the submarine.
The peak values of the nondimensional moments corresponding to different rotation centre positions are shown in Figure 22. With the increase in the distance from the rotation centre position to the bow of SUBOFF-AFF-1 and SUBOFF-AFF-8, the nondimensional moment peak value first decreases and then increases for all four cases. The moment value For SUBOFF-AFF-1 and SUBOFF-AFF-8 crossing the ISW, the nondimensional moment is the largest, which is the most dangerous condition.

Analysis of the Optimal Rotation Centre Position
According to the results of Section 4.3, the rotation centre position causes drastic variation in the nondimensional moment when the SUBOFF (AFF-1 and AFF-8) encounters the ISW. Therefore, it is meaningful to quantitatively explore the optimal range of the rotation centre position for the design of the submarine.
The peak values of the nondimensional moments corresponding to different rotation centre positions are shown in Figure 22. With the increase in the distance from the rotation centre position to the bow of SUBOFF-AFF-1 and SUBOFF-AFF-8, the nondimensional moment peak value first decreases and then increases for all four cases. The moment value in Case 2 is the largest. For SUBOFF-AFF-1, when the rotation centre position ranges from 0.46-0.48, the magnitude of the nondimensional moment is the smallest. For SUBOFF-AFF-8, when the rotation centre position ranges from 0.45-0.49, the magnitude of the nondimensional moment is the smallest. Considering that the four cases cover all the representative submergence areas, the optimal rotation centre position ranges for the two SUBOFFs can be recommended as 0.46-0.48 for SUBOFF-AFF-1 and 0.45-0.49 for SUBOFF-AFF-8. in Case 2 is the largest. For SUBOFF-AFF-1, when the rotation centre position ranges from 0.46-0.48, the magnitude of the nondimensional moment is the smallest. For SUBOFF-AFF-8, when the rotation centre position ranges from 0.45-0.49, the magnitude of the nondimensional moment is the smallest. Considering that the four cases cover all the representative submergence areas, the optimal rotation centre position ranges for the two SUB-OFFs can be recommended as 0.46-0.48 for SUBOFF-AFF-1 and 0.45-0.49 for SUBOFF-AFF-8.

Conclusions
In this paper, a numerical model of three-dimensional Navier-Stokes equations with a modified k-ω SST turbulence model combined with a density transport equation is developed to simulate the interaction of ISWs and submarines in continuously stratified, incompressible and viscous fluid. The model adopts the fully nonlinear models of the DJL

Conclusions
In this paper, a numerical model of three-dimensional Navier-Stokes equations with a modified k-ω SST turbulence model combined with a density transport equation is developed to simulate the interaction of ISWs and submarines in continuously stratified, incompressible and viscous fluid. The model adopts the fully nonlinear models of the DJL equation to generate ISWs in continuously stratified fluid. The developed model is used to simulate the internal solitary wave forces on the DARPA SUBOFF submarine model at different submergence depths. The following conclusions were obtained.
(1) The model developed in this paper has been verified through flow resistance experiments of a submarine model. For SUBOFF's total resistance, nondimensional pressure coefficient and nondimensional friction coefficient, the numerical simulation results of the presented model are highly consistent with the experimental results.
(2) When the SUBOFF is subjected to the action of an ISW, for the SUBOFF within the wave height range of the ISW, the magnitude of the vertical force of the SUBOFF is one order larger than that of the horizontal force, and the vertical force plays a leading role for the moment of the SUBOFF. When SUBOFF traverses the ISW, the vertical force is the largest, and the horizontal force is the smallest.
(3) When the submergence depth of SUBOFF is lower than the ISW waveform, the horizontal force and the vertical force are at the same order of magnitude, and the vertical force still plays a leading role in the moment of the submerged body.
(4) In the case of the ISW acting on SUBOFF, SUBOFF crossing the ISWs is the most dangerous working condition. At this time, the vertical force (downwards) and the moment (anticlockwise rotation) are the largest, and there is a danger of being dragged into the deep sea and causing overturning and crushing.
(5) With the continuous increase in distance from the position of the SUBOFF rotation centre to the bow, the positive nondimensional moment continues to decrease, gradually appearing in the form that first rotates anticlockwise and then rotates clockwise, and finally develops into a form in which a negative nondimensional moment dominates. There exists an optimal rotation centre position for the SUBOFFs corresponding to the smallest momentum. The following range of the optimal rotation centre position can be recommended for the SUBOFFs: S = 0.46-0.48 for a bare hull and S = 0.45-0.49 for a fully appended hull.
It is worth noting that only stationary submarines encountering ISWs were investigated in this paper. The submarine certainly goes into motion and even sinks in response to encountering ISW in the actual situation and it is necessary to consider the submarine motion for future research.

Institutional Review Board Statement:
The study was conducted according to the guidelines of the Declaration of Helsinki, and approved by the Institutional Review Board.
Informed Consent Statement: Informed consent was obtained from all subjects involved in the study.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.

Conflicts of Interest:
The authors declare no conflict of interest.