Numerical simulation on free motion response of a submarine induced by internal solitary wave



==1. Introduction==

Internal solitary waves (ISWs) are common phenomenon in the stratified ocean, which has the characteristics of medium-fine scale, large amplitude and strong velocity <sup>[1]</sup>. The shear currents and turbulence generated by ISWs contribute to global energy transport and dissipation in the ocean, but also affect biological activities, offshore engineering and submarine navigation <sup>[2]</sup>. Due to the development of ISWs below the sea surface, low propagation frequency and long distance, it is difficult to observe them which seriously threaten the safety of ocean engineering structures and underwater vehicle <sup>[3]</sup>.

In view of the influence of ISWs on ocean engineering and underwater vehicles, the current research focuses on the monitoring of ISWs <sup>[1, 4]</sup>, and the simulation of ISWs <sup>[3, 5]</sup>, as well as the influence of ISWs on submerged tunnels<sup>[6]</sup>, semi-submersible platforms <sup>[7, 8]</sup> and submarines <sup>[9-11]</sup>. Research on the interaction between ISWs and submarines has focused on the standard submarine model proposed by the Defense Advanced Research Projects Agency (DARPA), known as the DARPA SUBOFF model <sup>[12]</sup>. CHEN et al. studied the submarine load variation induced by ISWs with and without submarine speed <sup>[9]</sup>. LI et al. used the numerical model of ISWs in a continuously stratified fluid, numerically analyzed the forces and moments of SUBOFF when it encountered ISWs in different submergence depths, discussed the influence of rotation center position on the moment, and proposed the optimal rotation center position <sup>[10]</sup>. GUAN et al. numerically studied the hydrodynamic properties and load characteristics of submarine when it encountered ISWs in the ocean <sup>[11]</sup>.

It is worth noting that in the above studies, the freedom of the submarine is constrained when the submarine encounters the ISWs. The submarine position is fixed and in a static state. As matter of fact, the submarine will inevitably move or even sink when it encounters the ISWs. So the submarine movement must be taken into account in the research. DU et al. studied the kinematics characteristic and the mass center trace of the slender model in the density stratification fluid with a semi-flexible restriction technique, and obtained the properties such as heave, surge, roll, pitch, but the slender body was bound by the flexible line in the water tank <sup>[13]</sup>. CUI et al. used experimental methods to study the dynamic response of submerged body under the action of ISWs, and obtained the influence of different wave amplitudes, the thickness ratio of the upper and lower layers of fluid, and the shape of the submarine section, etc. Strong currents in the horizontal and vertical directions caused by ISWs can greatly move submerged body and cause damage to its structure. Limited by the experimental tank, the section of the submerged model is less than 10cm <sup>[14]</sup>.In general, the submerged body is affected by ISWs, which will lead to lateral, longitudinal and pitching motions.

In order to simulate the situation of large submarine encountering ISWs in the ocean, a large numerical tank and a large submarine model are constructed by numerical simulation method, and the dynamic response of submarine under the action of ISWs is studied in a numerical tank with continuous density changes.

==2 Simulation configurations==

===2.1 Submarine model===

The standard submarine model of Defense Advanced Research Projects Agency (DARPA), which called DARPA SUBOFF, has an overall length of 4.356m and a maximum diameter of 0.508m, as shown in Figure 1(a). The SUBOFF model itself is small in scale, so a 1:25.5 scale of DARPA SUBOFF is used to simulate the actual submarine. In order to simplify the model, the appendages of SUBOFF model were removed, as shown in Figure 1 (b). The enlarged bare hull has an overall length of 111m and a maximum diameter of 12.954m. The gravity center of the submarine is placed on the axis of the hull, 51 m away from the bow.

<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
 
{|
|-
| [[Image:Draft_sun_264386430-image1-c.png|276px]]
| [[Image:Draft_sun_264386430-image2-c.png|center|282px]]
|}
</div>

<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
(a) DARPA SUBOFF model            (b) SUBOFF simplified model</div>

<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
Figure 1 SUBOFF model</div>

===2.2 Numerical flume===

The numerical flume is shown in the Figure 2 and has a total length of 5000 m, width of 200 m and depth of 500 m, including the working zone which has a length of 4000 m and the wave damping zone which has a length of 1000 m. The fluid densities of the upper and lower layers are different. The upper layer fluid density is  <math display="inline">{\rho }_1=997kg/m^3</math> and has a depth of 100m (h1). The lower layer fluid density is  <math display="inline">{\rho }_2=1034kg/m^3</math> and has a depth of 400m (h2). The density varies continuously along depth in each layer of fluid, as shown in Figure 3.

<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
 [[Image:Draft_sun_264386430-image5.png|600px]] </div>

<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
Figure 2 Sketch of numerical flume</div>

<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
 [[Image:Draft_sun_264386430-image6.png|288px]] </div>

<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
Figure 3 Change of density along the water depth</div>

The distance between submarine location and the two layers interface is defined as the submergence depth H. The inlet and outlet are on the left side of the flume. Speed entrance method is adopted to make the numerical wave and pressure boundary condition is inducted for outlet. The inlet boundary is set by the user-defined function (UDF), which produces the ISWs. The bottom, top and right of the flume are all using the solid wall boundary.

===2.3 ISW theory===

The Korteweg-de Vires (KdV) equation has a wide range of applications for small amplitude and weak nonlinear ISWs <sup>[3]</sup>. For the numerical flume established in this paper, the ratio of wave amplitude to the depth of the flume  <math display="inline">{\eta }_0/(h_1+h_2)<0.1</math> ( <math display="inline">{\eta }_0</math> is the ISW amplitude), so the KdV theory can be used <sup>[15]</sup>.The propagation of a wave along the horizontal x direction is expressed as <sup>[16, 17]</sup>:

<div style="text-align: right; direction: ltr; margin-left: 1em;">
 <math display="inline">\frac{\partial \eta }{\partial t}+c_0\frac{\partial \eta }{\partial x}+</math><math>c_1\eta \frac{\partial \eta }{\partial x}+c_2\frac{{\partial }^3\eta }{\partial x^3}=</math><math>0</math> （1）</div>

where <math display="inline">c_0</math> ,  <math display="inline">c_1</math> and  <math display="inline">c_2</math> are the linear phase speed, nonlinear coefficient, measure of dispersion respectively, which can be wright as：

<div style="text-align: right; direction: ltr; margin-left: 1em;">
 <math>c_0={\left[\frac{gh_1h_2({\rho }_2-{\rho }_1)}{{\rho }_1h_2+{\rho }_2h_1}\right]}^{1/2}</math> （2）</div>

<div style="text-align: right; direction: ltr; margin-left: 1em;">
 <math>c_1=-\frac{3}{2}\frac{c_0}{h_1h_2}\frac{{\rho }_1h_2^2-{\rho }_2h_1^2}{{\rho }_1h_2+{\rho }_2h_1}</math> （3）</div>

<div style="text-align: right; direction: ltr; margin-left: 1em;">
 <math>c_2=\frac{c_0h_1h_2({\rho }_1h_1+{\rho }_2h_2)} {6({\rho }_1h_2+{\rho }_2h_1)}</math> （4）</div>

The KdV expression for the ISW is:

<div style="text-align: right; direction: ltr; margin-left: 1em;">
 <math display="inline">\eta (x,t)={\eta }_0{sech}^2(\frac{x-Vt}{\lambda })</math> （5）</div>

The soliton velocity and characteristic wavelength based on KdV theory are:

<div style="text-align: right; direction: ltr; margin-left: 1em;">
 <math>V=c_0(1+\frac{{\eta }_0c_1}{3c_0})</math> （6）</div>

<div style="text-align: right; direction: ltr; margin-left: 1em;">
 <math display="inline">\lambda ={\left(\frac{12c_2}{{\eta }_0c_1}\right)}^{1/2}</math> （7）</div>

===2.4 Meshing scheme and simulation method===

(1) Meshing scheme

The submarine will move horizontally and vertically, accompanied by pitching motion which affected by ISWs, as shown in Figure 4. The moving mesh is easy to be generated with negative volume, and the grid update takes a long time, resulting in a very long computation time. So the overset mesh technology is adopted in this paper. A separate component mesh was constructed in the SUBOFF region, the size of which was 160 m×20 m×60 m. The unstructured grids were adopted with total grids of 1.42 million. In addition, structured grid was used in the flume as the background grid with total grids of 5.31 million, as shown in Figure 5. In the overlapping regions, the overset interfaces connect cell zones by interpolating cell data. The overset mesh method can maintain mesh quality well during submarine motion.

<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
[[Image:Draft_sun_264386430-image19.png|246px]] </div>

<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
Figure 4 Sketch of submarine motion</div>


<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">

{|
|-
| [[Image:Draft_sun_264386430-image20-c.jpeg|246px]]
| [[Image:Draft_sun_264386430-image21.jpeg|center|288px]]
|}
</div>

<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
(a) Unstructered grids around the SUBOFF                      (b) Sketch of overset mesh</div>

<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
Figrue 5 The computational grid in the domain</div>

(2) Simulation method

The unsteady Reynolds- Averaged Naiver-Stokes equations with the <math display="inline">k-\epsilon </math> turbulence model are applied to simulate the flow field. The <math display="inline">k-\epsilon </math> turbulence added two variables to the governing equation which are turbulence kinetic energy  <math display="inline">k</math> and dissipation rate of turbulence kinetic energy  <math display="inline">\epsilon </math> .  <math display="inline">k</math> represents the change of velocity fluctuation and  <math display="inline">\epsilon </math> represents dissipation rate of velocity fluctuation. The mass continuity equation and momentum equation are shown as follows:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math display="inline">\frac{\partial \rho }{\partial t}+\nabla \cdot (\rho \boldsymbol{u})=</math><math>0</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (8)
|}

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math display="inline">\frac{\partial \rho \boldsymbol{u}}{\partial t}+\nabla \cdot (\rho \boldsymbol{u}\otimes \boldsymbol{u})\mbox{-}\nabla \cdot ({\mu }_{\mbox{eff}}\nabla \boldsymbol{u})\mbox{=}\nabla \cdot p^'+</math><math>\nabla \cdot {\left({\mu }_{\mbox{eff}}\nabla \boldsymbol{u}\right)}^T+</math><math>S_M</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (9)
|}


where  <math display="inline">{\mu }_{\mbox{eff}}</math> is effective viscosity; <math display="inline">p^'</math> is the modified pressure;  <math display="inline">S_M</math> is the sum of body forces.

The  <math display="inline">k</math> and  <math display="inline">\epsilon </math> can be solved by following equations:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math display="inline">\frac{\partial (\rho k)}{\partial t}+\nabla \cdot (\rho \boldsymbol{u}k)=</math><math>\nabla \cdot \left[\left(\mu +\frac{{\mu }_t}{{\sigma }_k}\right)\nabla k\right]+</math><math>P_k-\rho \epsilon </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (10)
|}

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math display="inline">\frac{\partial (\rho \epsilon )}{\partial t}+\nabla \cdot (\rho \boldsymbol{u}\epsilon )=</math><math>\nabla \cdot \left[\left(\mu +\frac{{\mu }_t}{{\sigma }_{\epsilon }}\right)\nabla \epsilon \right]+</math><math>C_{\epsilon 1}\frac{\epsilon }{k}P_k-C_{\epsilon 2}\rho \frac{{\epsilon }^2}{k}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (11)
|}


where  <math display="inline">{\sigma }_k</math> and  <math display="inline">{\sigma }_{\epsilon }</math> are Prandtl numbers corresponding to turbulence kinetic energy and dissipation rate, respectively; <math display="inline">P_k</math> is the turbulence production due to viscous forces;  <math display="inline">C_{\epsilon 1}</math> and  <math display="inline">C_{\epsilon 2}</math> are constants； <math display="inline">\mu </math> is the dynamic viscosity and  <math display="inline">{\mu }_\mbox{t}</math> is the turbulent viscosity.

The ISW wave surface was tracked by Volume of fluid (VOF) method <sup>[18]</sup>. The fluid’s volume fraction  <math display="inline">{\alpha }_i</math> represents the ratio of the  <math display="inline">i^{th}</math> fluid’s volume in a cell to the total volume of the cell.  For the  <math display="inline">{\alpha }_i</math> , the equation has the following form：

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math display="inline">\frac{\partial {\alpha }_\mbox{i}}{\partial t}+\frac{\partial (u{\alpha }_\mbox{i})}{\partial x}+</math><math>\frac{\partial (v{\alpha }_\mbox{i})}{\partial y}=</math><math>0</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (12)
|}


For every control volume, the volume fractions of all phases sum to unit. When  <math display="inline">{\alpha }_i\mbox{=1}</math> , the cell is full of the  <math display="inline">i^{th}</math> fluid；when  <math display="inline">{\alpha }_i\mbox{=}0</math> , the cell is empty of the <math display="inline">i^{th}</math>  fluid；when  <math display="inline">0<{\alpha }_i<1</math> ,the cell contains the interface <sup>[19]</sup>.For the two-layer fluid in this paper,  <math display="inline">i=1</math> or 2.

==3．Numerical results and discussions==

The submarine is located in the numerical flume and the distance between the gravity center and the flume inlet is 1500 m. The submarine is in a state of equilibrium with its own buoyancy. The power of the submarine is not taken into account, so it has no initial speed and its freedom of motion is not constrained.

In order to study the interactions between different ISWs and submarine, three ISW amplitudes and two submergence depths are used in this part, which include  <math display="inline">{\eta }_0\mbox{=10m}</math> , <math display="inline">{\eta }_0\mbox{=20m}</math> ,  <math display="inline">{\eta }_0\mbox{=30m}</math> and H=45m, H=15m,as shown in Table 1.

<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
Table 1 Calculation cases</div>

{| style="width: 70%;margin: 1em auto 0.1em auto;border-collapse: collapse;" 
|-
|  style="border-top: 1pt solid black;border-bottom: 1pt solid black;text-align: center;vertical-align: top;"|H（m）

<math display="inline">{\eta }_0</math> （m）
|  style="border-top: 1pt solid black;border-bottom: 1pt solid black;text-align: center;"|45
|  style="border-top: 1pt solid black;border-bottom: 1pt solid black;text-align: center;"|15
|-
|  style="border-top: 1pt solid black;text-align: center;vertical-align: top;"|10
|  style="border-top: 1pt solid black;text-align: center;vertical-align: top;"|C11
|  style="border-top: 1pt solid black;text-align: center;vertical-align: top;"|C12
|-
|  style="text-align: center;vertical-align: top;"|20
|  style="text-align: center;vertical-align: top;"|C21
|  style="text-align: center;vertical-align: top;"|C22
|-
|  style="border-bottom: 1pt solid black;text-align: center;vertical-align: top;"|30
|  style="border-bottom: 1pt solid black;text-align: center;vertical-align: top;"|C31
|  style="border-bottom: 1pt solid black;text-align: center;vertical-align: top;"|C32
|}


The case C31 was selected to verify the numerical flume established in this paper, and the formation of ISWs was monitored at x=1500m. As can be seen from Figure 6, the wave generation of the numerical flume was highly consistent with that of KdV theory, which met the requirements of numerical simulation.

<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
 [[Image:Draft_sun_264386430-image50.png|330px]] </div>

<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
Figure 6 Numerical result verification</div>

===3.1 Trajectory of submarine motion===

The initial submergence depth of SUBOFF was 45m or 15m below the interface. When the ISWs flow field was formed by the numerical flume, the motion of the submarine was observed. Figure 7 shows the motion trajectory of SUBOFF's center of gravity under different submergence depths. As can be seen from Figure 7, no matter the depth is 45m or 15m, SUBOFF's center of gravity moved in a similar trajectory, moving to the lower left under the influence of the ISWs. And the longitudinal displacement was significantly greater than the lateral displacement. Restricted by the size of the numerical flume, SUBOFF rebounded after impacting the bottom of the flume, and then continued to move to the upper left. At the same depth, the ISWs amplitude had no effect on the trajectory of submarine motion. Different wave amplitudes had no effect on the position of the impact point of SUBOFF, and the trajectory of the submarine basically coincided after rebound.

<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
 
{|
|-
| [[Image:Draft_sun_264386430-image51.png|288px]]
| [[Image:Draft_sun_264386430-image52.png|center|288px]]
|}
</div>

<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
:(1) H=45m (b) H=15m</div>

<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
Figure 7 Trajectory of submarine motion</div>

The change of submarine longitudinal position with time is shown in Figure 8. For the first 100s, SUBOFF kept its position unchanged. Because the flow field of ISWs had not reached the zone where the submarine was located, and the submarine was in a state of equilibrium. Then it began to dive under the influence of the ISWs flow field, and finally hit bottom and rebounded. As can be seen from Figure 8, with the increase of ISWs amplitude, the time required for the submarine to reach the bottom was longer, and the time required for C11, C21 and C31 to reach the bottom were 269s, 278s and 282s, respectively. According to KdV theory, the ISWs speed decreased with the increase of the ISWs amplitudes, which was the reason for the change of submarine bottoming time. The position of the submarine was 1500m away from the flume inlet, and it took less than 300s for the submarine sink more than 300m. It can be seen that the ISWs can cause serious damage to the submarine.

<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
 
{|
|-
| [[Image:Draft_sun_264386430-image53.png|288px]]
| [[Image:Draft_sun_264386430-image54.png|center|288px]]
|}
</div>

<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
:(1) H=45m (b) H=15m</div>

<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
Figure 8 Change of submarine longitudinal position with time </div>

In order to compare the influence of different submergence depth on the submarine motion trajectory, C31 and C32 with a wave amplitude of 30m were selected, as shown in Figure 9. The distance between the submarine and the flume bottom became greater with the decrease of the submergence depth, which lead to the lateral displacement was greater when the submarine hit the bottom, and the longitudinal displacement was smaller after the rebound which was affected by the ISWs. As can be seen from Figure 9 (b), as the depth H decreased, the time required for the submarine to reach the bottom increased.

<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
 
{|
|-
| [[Image:Draft_sun_264386430-image55.png|288px]]
| [[Image:Draft_sun_264386430-image56.png|center|288px]]
|}
</div>

<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
:(1) Motion trajectory (b) longitudinal position change with time</div>

<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
Figure 9 Trajectory of submarine motion with different H</div>

It is well known that the submarine will be dragged into the seabed by ISWs, and the impact on the bottom will produce strong load and lead to serious damage to its structure. The data after the rebound on the bottom has little significance. Therefore, this paper focuses on the analysis before the submarine hits the bottom.

===3.2 Characteristics of the motion===

The motion characteristics of the marine in the flume were studied when the wave amplitude was 20m and the submergence depth was 15m (C22). Figure 10 shows the movement of SUBOFF at six moments before it hits the bottom. The color mark the upper and lower layers of fluid with different densities respectively, and the submarine is below the interface. During the first 100s, the position of the submarine remained basically unchanged. After that, the submarine began to sink, and Figure 10(d) showed obvious sinking of the submarine. At t=200s, it can be seen that the submarine was starting to pitch, the submarine's bow was up. At t=250s, the change of pitching and sinking was obvious. At this time, the ISWs peak has not reached the location of the submarine. The flow field formed by ISWs had a strong influence on the submarine before the peak of the ISW reached the position of the submarine, resulting in not only submarine transverse and longitudinal displacement but also the pitching.

<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
 
{|
|-
| [[Image:Draft_sun_264386430-image57.jpeg|246px]]
| [[Image:Draft_sun_264386430-image58.jpeg|center|246px]]
|}
</div>

<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
(a) 0s (b)50s</div>

<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
 
{|
|-
| [[Image:Draft_sun_264386430-image59.jpeg|246px]]
| [[Image:Draft_sun_264386430-image60.jpeg|center|246px]]
|}
</div>

<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
(c)100s (d)150s</div>

<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
 
{|
|-
| [[Image:Draft_sun_264386430-image61.jpeg|246px]]
| [[Image:Draft_sun_264386430-image62.jpeg|center|246px]]
|}
</div>

<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
(e)200s (f)250s</div>

<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
Figure 10 Submarine motion characteristics at different times</div>

===3.3 Motion parameters===

(1)The effect of different ISWs amplitudes

The lateral velocity of the submarine under different ISW amplitudes is shown in Figure 11. After t=100 s, the absolute lateral value of velocity increased along the negative x-direction. Under the same submergence depth, lateral acceleration decreased with the increase of wave amplitude. When H=45 m and t=260 s, the lateral velocities of C11, C21 and C31 were -0.51 m/s, -0.41 m/s and -0.39 m/s, respectively. By comparison with Figure11 (a) and Figure 11 (b), it is found that the submergence depth has little influence on the lateral velocity variation trend.

<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
 
{|
|-
| [[Image:Draft_sun_264386430-image63.png|288px]]
| [[Image:Draft_sun_264386430-image64.png|center|288px]]
|}
</div>

<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
:(1) H=45m (b) H=15m</div>

<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
Figure 11 Lateral velocity at different ISW amplitudes</div>

The longitudinal velocity of the submarine under different depths is shown in Figure 12. Similar to the lateral velocity, the longitudinal velocity increased along the negative y-direction, and the change of longitudinal velocity was one order of magnitude larger than that of the lateral velocity. Before t=180s, the submarine was accelerating longitudinally. After t=180s, the longitudinal velocity began to keep basically stable at about -2.5m/s.

<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
 
{|
|-
| [[Image:Draft_sun_264386430-image65.png|288px]]
| [[Image:Draft_sun_264386430-image66.png|center|288px]]
|}
</div>

<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
:(1) H=45m (b) H=15m</div>

<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
Figure 12 Longitudinal velocity at different ISW amplitudes</div>

The pitch angle of the submarine under different ISW amplitudes is shown in Figure 13. When the submergence depth was 45m, as shown in Figure 13(a), the submarine basically keeps level without pitching before t=180s; after t=180s, the submarine began to pitch. The pitch angle gradually increased with time, and the pitch angle variation decreased with the increase of ISW amplitudes. When t=260s, the pitch angles of C11, C21 and C31 were -21.10°, -17.16° and -14.50°, respectively. When the submergence depth H=15m, the changes of pitch angle were similar to H=45m.

<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
 
{|
|-
| [[Image:Draft_sun_264386430-image67.png|288px]]
| [[Image:Draft_sun_264386430-image68.png|center|288px]]
|}
</div>

<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
:(1) H=45m (b) H=15m</div>

<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
Figure 13 Pitch angle at different ISW amplitudes</div>

（2）The effect of submergence depths

Figure 14 shows the comparison of the submarine velocity under different submergence depths. When the wave amplitude was constant, the influence of the submergence depth on the change of the lateral velocity and longitudinal velocity can be ignored.

<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
 
{|
|-
| [[Image:Draft_sun_264386430-image69.png|288px]]
| [[Image:Draft_sun_264386430-image70.png|center|288px]]
|}
</div>

<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
(a) Comparison of lateral velocity  (b) Comparison of longitudinal velocity </div>

<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
Figure 14 Comparison of velocity at different depths</div>

The variations of pitch angle with time at the same wave amplitude but different submergence depths are shown in Figure 15. As shown in the figure, submergence depths had little influence on the variation of pitch angle.

<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
 [[Image:Draft_sun_264386430-image71.png|288px]] </div>

<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
Figure15 Comparison of pitch angle at different depths</div>

==4 Conclusions==

The dynamic response of large submarine induced by ISWs was simulated by establishing the ISWs numerical flume. The flume was a stratified fluid of two layers, and the fluid density varied with the depth which was more in line with the actual ocean. The track of gravity center, variations in lateral, longitudinal and pitch motion are studied. The conclusions of this study are listed as follows:

(1) When the submarine is located under the two-layer interface, the submarine moves to the lower left, then impacts on bottom and rebounds which is affected by ISWs. The submarine motion trajectory and the impact location are not affected by the initial submergence depth.

(2) As the submergence depth decreases, the distance between the initial position of the submarine and the bottom increases. Under the same ISW, the time for the submarine impacting the flume bottom is longer. The ISW amplitude also affects the time for the submarine impacting the flume bottom, and as the ISW amplitude increases, the longer the time is required.

(3) The submarine will move which is affected by ISWs. The longitudinal velocity is bigger than lateral velocity. The motion of the submarine is accompanied by pitching. The influence of submergence depth on the submarine lateral velocity and longitudinal velocity is negligible. Before the submarine reaches the bottom, the horizontal velocity is increasing. The higher the ISW amplitude, the smaller the acceleration is. The longitudinal velocity also accelerates first, and then the value fluctuation is maintained. In the process of the sinking, the pitch angle gradually becomes larger, and the pitch angle decreases as the ISW amplitude increases. The changes in the pitch angle are not affected by the submergence depth.

==Conflict of interest statement==

The authors declare that there are no conflict of interests, we do not have any possible conflicts of interest.

==Funding==

This work was supported by Natural Science Foundation of Sichuan Province of China No. 2022NSFSC0478, CAFUC Foundation J2022-035 and  Foundation of Key Laboratory of Flight Techniques and Flight Safety，CAAC, FZ2020ZZ05.

==References==

<span id='_nebA1B3473E_F351_4D32_B948_3DA0003B2C7B'></span>[1] Huang X, Chen Z, Zhao W, et al. An extreme internal solitary wave event observed in the northern South China Sea. Scientific Reports, 2016,6(1):30041.

<span id='_nebDA44E60A_2C14_425A_BBE8_EB5D1EB74680'></span> [2] Guo C, Chen X. A review of internal solitary wave dynamics in the northern South China Sea. Progress in oceanography, 2014,121:7-23.

<span id='_neb9E7607E2_8CC5_4704_8EDC_65232DC3188C'></span> [3] Cui J, Dong S, Wang Z. Study on applicability of internal solitary wave theories by theoretical and numerical method. Applied Ocean Research, 2021,111:102629.

<span id='_nebB0B99C6E_F214_49DB_8071_899C36AE0AD7'></span> [4] Alford MH, Lien R, Simmons H, et al. Speed and evolution of nonlinear internal waves transiting the South China Sea. Journal of Physical Oceanography, 2010,40(6):1338-1355.

<span id='_nebEA94831D_DAD0_43BF_ADD5_AE940F12B5B4'></span> [5] Chen CY, Hsu RC, Chen HH, et al. Laboratory observations on internal solitary wave evolution on steep and inverse uniform slopes. Ocean Engineering, 2007,34(1):157-170.

<span id='_neb70BE7145_9EDC_4C2E_9BA7_36B1DEEAC8FB'></span> [6] Zou PX, Bricker JD, Uijttewaal WSJ. The impacts of internal solitary waves on a submerged floating tunnel. Ocean Engineering, 2021,238:109762.

<span id='_neb6406592C_E197_46DF_A377_D720346DEC31'></span> [7] Wang X, Zhou J, Wang Z, et al. A numerical and experimental study of internal solitary wave loads on semi-submersible platforms. Ocean engineering, 2018,150:298-308.

<span id='_neb7E1675DA_6A67_46D0_BD31_11035A3811E1'></span> [8] Ding W, Ai C, Jin S, et al. 3D Numerical investigation of forces and flow field around the semi-submersible platform in an internal solitary wave. Water, 2020,12(1):208.

<span id='_nebC1AB4008_0311_447F_B3B9_585AFCDB0050'></span>[9] Chen J, You YX, Liu X D, et al. Numerical simulation of interaction of internal solitary waves with a moving submarine. Chinese Journal of Hydrodynamics, 2010.

<span id='_nebF9F37563_C5E0_47BA_A3E6_AB43CB4808D0'></span>[10] Li J, Zhang Q, Chen T. Numerical investigation of internal solitary wave forces on submarines in continuously stratified fluids. Journal of marine science and engineering, 2021,9(12):1374.

[11] Guan H, Wei G, Hui DU. Hydrodynamic properties of interactions of three-dimensional internal solitary waves with submarine. Journal of PLA University of Science and Technology (Natural Science Edition), 2012,13(05):577-582.

[12] Groves NC, Huang TT, Chang MS. Geometric characteristics of DARPA SUBOFF models (DTRC Model Nos. 5470 and 5471). 1989, Technical Report (No.DTRC/SHD-1298-01).

[13] Du H , Wei G , Zeng WH , et al. Experimental investigation on the kinematics characteristic of submerged slender body under internal solitary wave of depression. Journal of Ship Mechanics, 2017, 21(10):1210-1217.

<span id='_neb0F8BB81C_4ACB_4A46_A5D8_18B1A98389B0'></span>[14] Cui J, Dong S, Wang Z, et al. Kinematic response of submerged structures under the action of internal solitary waves. Ocean Engineering, 2020,196:106814.

<span id='_nebEC3BF0DB_6ED8_4811_B654_126EBF70E930'></span>[15] Koop CG, Butler G. An investigation of internal solitary waves in a two-fluid system. Journal of Fluid Mechanics, 1981,112(-1):225.

<span id='_neb5D2BFFF9_2F16_43A6_AB2F_57B09F6E6048'></span>[16] Korteweg DJ, de Vries G. XLI. On the change of form of long waves advancing in a rectangular canal, and on a new type of long stationary waves. The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science, 1895,39(240):422-443.

<span id='_neb7CE8CB30_1770_4BDA_8C4F_A16DCD76F9E1'></span>[17] Ostrovsky LA, Stepanyants YA. Do internal solitons exist in the ocean?. Reviews of Geophysics, 1989,27(3):293, 310.

<span id='_neb84F0F059_30BA_4914_B0DA_14413D389C9F'></span>[18] Hirt CW, Nichols BD. Volume of fluid (VOF) method for the dynamics of free boundaries. Journal of Computational Physics, 1981,39(1):201-225.

<span id='_nebD785F48B_4786_4E6C_9DEE_7ABA5D7FCD5E'></span>[19] Ma W, Li Y, Ding Y, et al. Numerical investigation of internal wave and free surface wave induced by the DARPA SUBOFF moving in a strongly stratified fluid. Ships and offshore structures, 2020,15(6):587-604.