Numerical analysis of multi-scale pressure pulsation on the energy accumulation for submarine-based tracking and pointing systems

Abstract

One of the main goals of submarine designers and researchers is to estimate the influence of submarine fluid dynamics for submarine-based optical tracking and pointing systems. In this study, firstly, based on the basic flow governing equation and hierarchical grids, the numerical simulation method of DNS (direct numerical simulation) is adopted to simulate the seawater flow around the submarine at 6° yaw angle and 10⁷ Reynolds number. Secondly, the transformation equations from the earth coordinate system to the optical axis system have been deduced and the ultimate influence of pressure torques on the tracking system is studied. Transfer functions of the coarse channel direct current (DC) torque motor and fine channel fast steering mirror (FSM) also have been modeled and deduced. On this basis, the time domain step responses of both subsystems are analyzed by MATLAB Simulink. Finally, performance analyses have been deduced by comparing the error variation and vortices evolution. It revealed that the frequency characteristics of multi-scale pressure pulsation mainly depended on the lengths of submarine hull or its appendage, as well as the fluid dissipation and random interaction. In general, the coarse channel appears a good compensation performance at low frequency and large amplitude error that caused by the middle-scale pressure pulsation. Contrarily, the FSM fine channel exerts an excellent control effect for higher frequency and small amplitude error caused by small-scale pressure pulsations.

Introduction

Different from conventional photoelectric tracking system, laser directed energy device needs to guide and keep the laser beam continuously and stably focused on the target. ¹ The tracking accuracy must be one or two orders of magnitude higher than the conventional weapon. Originally, only equipped with periscope, submarine-based photoelectric mast has become an advanced image acquisition and automatic tracking device controlled with computer in recent years. ² In further, laser directed energy device is a promising application on submarines with the continue miniaturization of high energy laser equipment. Therefore, it is necessary to deeply study the mechanism of vortex induced vibration for the disturbance error of submarine-based laser directed energy device.

When submarine moves at high speed underwater, the boundary-layer separation and flow field perturbation around the submarine and its appendages will be generally formed. In the case of high Reynolds number, tip vortexes will be generated at the end edge of the hull or appendages of the submarine accompanied by large force and torque fluctuations. In hydrodynamics, large low pressure area will be formed in vortex cores, and cavitation pittings will be aroused once the pressure drops to a certain degree. ³ Not only will strong pulsations produce extremely unfavorable vibration and noise, but it also reduces the stability or concealment of the military submarine.

Besides tip vortices, turbulence vortices characterized by different multi-scales and shapes are induced from the submarine body and its appendages. Depending on their scales, the evolution of turbulence could mainly be divided into three stages. The first is the slowly decay period of large-scale vortices, during which stage the large-scale vortices decay slowly and mainly affected by the fluid’s convection and diffusion interactions rather than the dissipation. The second is the fast decay stage of medium scale vortices, during which period the large-scale vortices cannot keep their former shapes and then fast broken with the existing time being the least for the whole period. The third is the small-scale vortices dissipation stage; the fluid’s convection tends to weak and the convection time is the longest.

The turbulent vortex always changes randomly and rapidly over time and space. The interactions of turbulent and hull will lead the problematic noises. On the one hand, the random pressure pulsation directly produces radiated noise; on the other hand, the elastic structure of hull’s surface material will vibrate and produce radiated noise, which is collectively classify as “hydrodynamic noise”. ⁴ Unfortunately, the vibration or noise inevitably restricts the tracking precision of the laser directed energy device. ⁵

At present, a great deal of investigations have been carried out on submarines under different navigational conditions. However, numerous investigations have been widely conducted focused on submarines, under various navigation conditions. Predictive capabilities of hybrid RANS/LES models were compared for single-phase flow over a surface combatant at Re = 5.3 × 10⁶ and an appended the DARPA SUBOFF model at Re = 1.2 × 10⁷. ⁶ Alin et al. ⁷ discussed some features of submarine hydrodynamics, engaging in further studies on the DARPA SUBOFF configurations a bare hull and a fully appended model. Based on the differ DARPA Suboff models, Xing et al. ⁸ identified and analyzed the vortical and turbulent structures, instabilities, with analogy to the vortex breakdown and helical instability analysis for delta wing at large angles of attack.

The hydrodynamic characteristics of the same model are similar to the disturbance mechanism, so we only analyzed the condition of Re = 10⁷ and attitude at 6° yaw angle. In the third section of this paper, based on the basic control equations of fluid mechanics and the hierarchy grids, seawater flow around the submarine has been numerical simulated. In section 4, coordinate transformation methods have been deduced. In section 5, based on the multi-scale visual vortices, the torque disturbance in optical axis coordinate system has been investigated. The sixth section deduces the mathematical models of torque motor and FSM, according to the design mechanisms of torque motor and permanent magnet synchronous motor. And then by MATLAB Simulink toolbox, the effect of submarine disturbance on the coarse/fine channels output errors are simulated and analyzed. In section 7, according to the vortex field evolution characters, the influence factors and internal mechanism between the pressure pulsation and tracking precision of submarine-based optical system have been further revealed.

System descriptions

Except for algorithm filtering, adaptive optical frameworks are generally adopted to compensate the disturbance angle error in order to enhance the anti-disturbance capability as well as the stable tracking capability of tracking and aiming servo system. Recently, FSM, as an important optical path adjusting device, has been widely used for precise adjustment and stabilization control of optical axis, ⁹ which is primarily composed of piezoelectric ceramic stacks, a closed loop controller and a loaded mirror. Compared with the torque motors in coarse channel, the FSM in fine channel offer several advantages such as higher angle control precision (the repeatable positioning accuracy is about 5–20 urad) and high operation bandwidth (about 0.2–1 kHz), which will be discussed in the following paper.

Figure 1(a) shows the laboratory servo equipment we used, which is composed of coarse and fine two channels. The servo equipment is primarily composed of two coarse channel azimuth/pitch rotary platforms, two torque motors, a FSM, a charged coupled device (CCD) and a wedge-shaped body with 45° inclined plane. During the experiment, the laser is generated from a small laser source loaded on the right platform in front of the FSM. Laser beam is reflected by FSM, and then passes through the round hole. Finally, through 45° wedge mirror the laser is adjusted and shot toward the target. Two DC torque motors are selected to drive the altitude and the azimuth platforms, avoiding the gear clearance error.

Figure 1.

Diagram of submarine tracking and aiming equipment: (a) structure of tracking and aiming system and (b) conception of tracking and aiming servo system on the submarine.

In practical, servo equipment in Figure 1(a) will be loaded at the optical mast in Figure 2(b). Firstly, the target is detected by radars or photoelectric detection devices to guide the submarine-based tracking and aiming equipment and catch the target. Secondly, when the coarse channel intends to track stably, the fine channel FSM will be enabled and guide the laser shot on the target more precisely, as shown in Figure 1(b).

Figure 2.

Submarine model in numerical simulation.

Basic model and numerical method

Basic model

Figure 2 shows the submarine model’s scales, which is composed of a semicircular head, cylindrical body and pyramidal tail, with corresponding lengths of 0.06l, 0.69l, and 0.25l, respectively, where l is the submarine length and also the characteristic length in this simulation.

Total size of computational domain is 4l × 2l × 2l, free stream flow along the X-axis positive direction, the leading edge point of the submarine located at the origin of right-handed coordinate system and keeping 0.5l away from the inlet surface. The X-axis points to the downstream, the Y-axis points to the starboard, and the Z-axis points to the vertically upwards. In this simulation, for the computational domain, the left and right surfaces are defined as speed entry and pressure outlet boundary conditions with constant pressure and velocity in normal direction being zero. Additional, the submarine surface is defined as no-slip boundary condition.

Numerical method

Governing equations

In the computational domain the dimensionless viscous incompressible flow governing equations can be written as:

$$\nabla ·\mathit{U}=0$$

(1)

$$\frac{\partial \mathit{U}}{\partial t}+\mathit{U}·\nabla \mathit{U}=-\nabla p+\frac{1}{\mathrm{Re}}{\nabla }^2\mathit{U}$$

(2)

$\mathit{U}=(u,v,w)$ is the velocity vector and $p$ is the pressure. Reynolds number $\mathrm{Re}=u_{\infty }l/\upsilon$ , $u_{\infty }$ is the inflow velocity, $\upsilon$ is fluid kinematic viscosity coefficient, dimensionless time $t=t\prime u_{\infty }/l$ the $t\prime$ computing time. Here $u_{\infty }=l=1$ so $t=t\prime$ .

Spatial discretization

The computer domain is spatially discretized using square (cubic in 3D) finite volumes organized hierarchically as a quad-tree (octree in 3D).^10,11 This type of discretization in Gerris software has been applied to solve the governing equations for incompressible flows. ¹² Detail descriptions to spatial discretization and the corresponding tree representation are given in Figure 3 and a prior paper. ¹³

Figure 3.

Diagram of quad-tree discrete and corresponding tree structure.

Using the tree structure grid, for the $N_{\mathrm{level}}$ layer grids, the total number of grids should be $8^{{\mathrm{N}}_{\mathrm{level}}}$ in unit volume. However, the total grid number will be tremendously reduced by using mesh adaptive method than the uniform grid generation method. The mesh adaptive method can be described as

$$Vorticity\frac{h\left|\nabla \times u\right|}{\max (\left|u\right|)}<ϵ$$

(3)

This criterion will ensure that a finer resolution is used in solution areas with higher vorticity values. h the mesh size and $ϵ$ a user-defined threshold. The threshold $ϵ$ can be interpreted as the maximum angular deviation from a straight path, due to the local vorticity $\nabla \times \mathit{u}$ of a massless particle traveling at speed $max(\left|\mathit{u}\right|)$ across a cell of size h. For interfacial flows, vorticity generation is often associated with high interface curvatures. This criterion is used to ensure that small radii of curvature can be resolved properly. ¹⁴ When the maximum speed residuals less than a settled threshold value, the iterative program will stop.

Temporal discretization

The classical projection method is commonly referred to solve these government equations. Base on this Helmholtz-Hodge decomposition that a twice continuously differentiable bounded vector field can be decomposed into a curl-free component as well as another divergence-free component. ¹³ This projection method relies on Hodge decomposition of the intermediate velocity field U * as. ¹⁵

$${\mathit{U}}^{*}={\mathit{U}}^{n+1}+\nabla \phi$$

(4)

The $\phi$ is referred to the Hodge variable. The $\nabla ·\mathit{U}=0$ in the computational domain $\Omega$ and $\mathit{U}·\mathit{n}=0$ on boundary $\partial \Omega$ . Then we will apply the divergence operator on both sides of this equation.

$$\nabla ·\nabla \phi =\nabla ·{\mathit{U}}^{*}$$

(5)

By solving the Poisson equation the Hodge variable $\phi$ can be obtained. And then the intermediate velocity field U * is projected onto the divergence-free subspace to obtain ${\mathit{U}}^{n+1}$ by the following equation.

$${\mathit{U}}^{n+1}={\mathit{U}}^{*}-\nabla \phi$$

(6)

The discrete form of the projection operator depends on the discrete position of the corresponding pressure and velocity fields. The convective pressure was projected accurately at the center of the cell surface, while the velocity was selected approximately at the center of the cell.

Therefore, based on the above solution method, at time step n, the control equation can be discretized as:

$$\begin{array}{c}\multicolumn{1}{c}{\frac{{\mathit{U}}^{n+1}-{\mathit{U}}^n}{\mathrm{\Delta }t}+{\mathit{U}}^{n+\frac{1}{2}}·\nabla {\mathit{U}}^{n+\frac{1}{2}}=-\nabla {\mathbf{p}}^{n+\frac{1}{2}}}\\ \multicolumn{1}{c}{+\frac{1}{\mathrm{Re}}{\nabla }^2({\mathit{U}}^{n+1}+{\mathit{U}}^n)}\end{array}$$

(7)

$$\nabla ·{\mathit{U}}^n=0$$

(8)

$\mathrm{\Delta }t$ is time step. The equations can be further simplified by using the above time division projection method:

$$\frac{{\mathit{U}}^{*}-{\mathit{U}}^n}{\mathrm{\Delta }t}+{\mathit{U}}^{n+\frac{1}{2}}·\nabla {\mathit{U}}^{n+\frac{1}{2}}=\frac{1}{\mathrm{Re}}{\nabla }^2({\mathit{U}}^{*}+{\mathit{U}}^n)$$

(9)

$${\mathit{U}}^{n+1}={\mathit{U}}^{*}-\mathrm{\Delta }t\nabla p^{n+\frac{1}{2}}$$

(10)

$$\nabla ·{\mathit{U}}^{n+1}=0$$

(11)

Poisson equation can be written based on the equations (10) and (11).

$$\nabla ·(\mathrm{\Delta }t\nabla p^{n+\frac{1}{2}})=\nabla ·{\mathit{U}}^{*}$$

(12)

Poisson equation is solved by an Octree-based multistage solver. Discrete momentum equation (9) can be rewritten as:

$$\frac{{\mathit{U}}^{*}}{\mathrm{\Delta }t}-\frac{1}{\mathrm{Re}}{\nabla }^2{\mathit{U}}^{*}=\frac{1}{\mathrm{Re}}{\nabla }^2{\mathit{U}}^n+\frac{{\mathit{U}}^n}{\mathrm{\Delta }t}-{\mathit{U}}^{n+\frac{1}{2}}·\nabla {\mathit{U}}^{n+\frac{1}{2}}$$

(13)

$$\frac{{\mathit{U}}^{*}}{\mathrm{\Delta }t}-\frac{1}{\mathrm{Re}}{\nabla }^2{\mathit{U}}^{*}=\frac{1}{\mathrm{Re}}{\nabla }^2{\mathit{U}}^n$$

(14)

The convection term is discretized by second-order upwind scheme of Bell-Colella-Glaz. The diffusive term is discretized by Crank-Nicholson method, which has second order accuracy and unconditional stability. Above all, in this simulation, both time and space dispersions are of the second order accuracy. ¹⁶ Direct numerical simulation, or DNS, a full numerical simulation of a turbulent flow is adopted, which does not rely on any numerical calculation models of the flow.

Grid independence test and experimental verification

In order to verify that the grid layer number has no effect on the results as well as the correctness of the numerical simulation method, the calculation grids and results for four different kinds of grid layer numbers (N_level = 6, N_level = 7, N_level = 8, and N_level = 9) have been numerical investigated based on the SUBOFF bare model are shown in Figure 4. ¹⁷

Figure 4.

Pressure coefficient distributions of the submarine middle section under different grid layers (Exp. is experimental result). ¹⁸

Figure 4 shows that pressure coefficient C_p distributions on the section of the hull with four different grid layers. As shown, these C_p distributions have relatively consistent and closely change trend with the experimental results. And with the increasing of grid layers, the calculated pressure coefficients tend to more precise consistent with the experimental result ¹⁸ . In addition, based on this flow field calculation method, the flow around a cylinder has be calculated and compared with previous research results at Re = 100, which are available in our former papers.^19,20

In order to detect the fluid dynamics microscopic characteristic of the submarine more sensitive, a balance between the computer condition and the calculation complexity should be struck. Therefore, nine layer grids are adopted in this study. In this numerical simulation, the smallest length of grid is $\mathrm{l}/2^{N_{\mathrm{level}}}$ (N_level = 9), with the total number of grid being about 5,570,000 at the final time step.

Perturbation solution of the azimuth/pitch angles

It could be assumed that, for very short CCD sampling period $\mathrm{\Delta }t$ (0.01 s–0.02 s), the influence of torque disturbance would not make the target beyond the CCD view field. In other words, the target information can always be captured.

Deduction of translation movements

Figure 5 is diagram of translational motion decoupled. The influences of this motion along X-axis, Y-axis and Z-axis on the azimuth and pitch angles can be described as equations (15) to (17):

$$\{\begin{array}{c}\mathrm{\Delta }{A}_{\mathrm{d}x}=\sin A_{\mathrm{d}}·v_x\mathrm{\Delta }t/l_1\hfill \\ \multicolumn{1}{c}{\mathrm{\Delta }{E}_{\mathrm{d}x}=\cos A_{\mathrm{d}}\sin E_{\mathrm{d}}·v_x\mathrm{\Delta }t/L_1}\end{array}$$

(15)

$$\{\begin{array}{c}\mathrm{\Delta }{A}_{\mathrm{d}y}=-v_y\cos A_{\mathrm{d}}·v_y\mathrm{\Delta }t/l_1\hfill \\ \multicolumn{1}{c}{\mathrm{\Delta }{E}_{\mathrm{d}y}=\sin A_{\mathrm{d}}\sin E_{\mathrm{d}}·v_y\mathrm{\Delta }t/L_1}\end{array}$$

(16)

$$\mathrm{\Delta }{E}_{\mathrm{d}z}=-\cos E_{\mathrm{d}}·v_z\mathrm{\Delta }t/L_1$$

(17)

Figure 5.

Diagram of translational motion decoupled.

Through ignoring the coupling interaction from three different translational disturbances, the azimuth and pitch angles can be written as:

$$\{\begin{array}{c}\multicolumn{1}{c}{A_{\mathrm{d}}^{\prime }=A_{\mathrm{d}}+(\sin A_{\mathrm{d}}·v_x-\cos A_{\mathrm{d}}·v_y)·\mathrm{\Delta }t/l_1}\\ \multicolumn{1}{c}{E_{\mathrm{d}}^{\prime }=E_{\mathrm{d}}+(\cos A_{\mathrm{d}}·v_x+\sin A_{\mathrm{d}}·v_y-\cos A_{\mathrm{d}}·v_z)}\\ \multicolumn{1}{c}{\sin E_{\mathrm{d}}·\mathrm{\Delta }t/L_1}\end{array}$$

(18)

$A_{\mathrm{d}}$ , $E_{\mathrm{d}}$ and $A{\prime }_{\mathrm{d}}$ , $E{\prime }_{\mathrm{d}}$ are the initial and transformed azimuth and pitch angles, respectively. $\mathrm{\Delta }{A}_{\mathrm{di}}$ and $\mathrm{\Delta }{E}_{\mathrm{di}}$ are the changes of the azimuth and pitch angles at the i th(i = x, y, z) components. $v_x,v_y,v_z$ are the motion velocities of the pitch, swaying and heaving. $L_1$ is the straight distance from the device to the target. $l_1$ is the distance from the submarine to the target’s shadow on the sea.

In tracking and aiming systems, the necessary frame frequency for image processing $1/\mathrm{\Delta }t$ is about 50–100 Hz. Compare to the far distance target distance L₁ the $v_x\mathrm{\Delta }t$ , $v_{\mathrm{y}}\mathrm{\Delta }t$ , $v_{\mathrm{z}}\mathrm{\Delta }t$ are too small, making the angle distributions from the translational motion can be ignored. So in this study we only discuss the torque disturbance.

Deduction of swing movements

The azimuthal angle and pitch angle of target G (x, y, z) in submarine coordinate $(X_{\mathrm{d}},Y_{\mathrm{d}},Z_{\mathrm{d}})$ could be solved by ${\alpha }_{\mathrm{H}}$ , ${\alpha }_{\mathrm{P}}$ , ${\alpha }_{\mathrm{R}}$ (the yaw angler ${\alpha }_{\mathrm{H}}$ , pitch angle ${\alpha }_{\mathrm{P}}$ and roll angle ${\alpha }_{\mathrm{R}}$ ) and $\mathrm{\Delta }t$ , independently. The diagram of rotational movement decoupled presented in Figure 6.

Figure 6.

Diagram of rotational movement decoupled (yaw).

Three steps needed when convert the coordinates of the target point in the earth coordinate system $(\mathrm{X},\mathrm{Y},\mathrm{Z})$ to the submarine coordinate system $(X_{\mathrm{d}},Y_{\mathrm{d}},Z_{\mathrm{d}})$ . In the earth system the X-axis points to the east, the Y-axis points to the north and the Z-axis points to the sky, the origin is located at the meeting point of submarine axis and hull axis. Before the rotation, the submarine coordinate system $(X_{\mathrm{d}},Y_{\mathrm{d}},Z_{\mathrm{d}})$ is totally overlapped with the earth system. (1) Rotate ${\alpha }_{\mathrm{H}}$ angle around OZ (yaw) as Figure 6 shows. (2) Rotate ${\alpha }_{\mathrm{P}}$ around OY (pitch). (3) Rotate ${\alpha }_{\mathrm{R}}$ around OX axis (roll).

Three transitions could be solved in the matrix form as following equations (19) to (21):

$$\left[\begin{array}{c}{X}_{\mathrm{H}}\hfill \\ \multicolumn{1}{c}{Y_{\mathrm{H}}}\\ \multicolumn{1}{c}{Z_{\mathrm{H}}}\end{array}\right]=\left[\begin{array}{ccc}\cos {\alpha }_{\mathrm{H}}\hfill & \sin {\alpha }_{\mathrm{H}}\hfill & 0\hfill \\ \multicolumn{1}{c}{-\sin {\alpha }_{\mathrm{H}}}& \cos {\alpha }_{\mathrm{H}}\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 1\hfill \end{array}\right]\left[\begin{array}{c}X\hfill \\ \multicolumn{1}{c}{Y}\\ \multicolumn{1}{c}{Z}\end{array}\right]=T_{\mathrm{H}}\left[\begin{array}{c}X\hfill \\ \multicolumn{1}{c}{Y}\\ \multicolumn{1}{c}{Z}\end{array}\right]$$

(19)

$$\left[\begin{array}{c}{X}_{\mathrm{P}}\hfill \\ \multicolumn{1}{c}{Y_{\mathrm{P}}}\\ \multicolumn{1}{c}{Z_{\mathrm{P}}}\end{array}\right]=\left[\begin{array}{ccc}\cos {\alpha }_{\mathrm{P}}\hfill & 0\hfill & \sin {\alpha }_{\mathrm{P}}\hfill \\ \multicolumn{1}{c}{0}& 1\hfill & 0\hfill \\ \multicolumn{1}{c}{-\sin {\alpha }_{\mathrm{P}}}& 0\hfill & \cos {\alpha }_{\mathrm{P}}\hfill \end{array}\right]\left[\begin{array}{c}{X}_{\mathrm{H}}\hfill \\ \multicolumn{1}{c}{Y_{\mathrm{H}}}\\ \multicolumn{1}{c}{Z_{\mathrm{H}}}\end{array}\right]=T_{\mathrm{P}}\left[\begin{array}{c}{X}_{\mathrm{H}}\hfill \\ \multicolumn{1}{c}{Y_{\mathrm{H}}}\\ \multicolumn{1}{c}{Z_{\mathrm{H}}}\end{array}\right]$$

(20)

$$\left[\begin{array}{c}{X}_{\mathrm{d}}\hfill \\ \multicolumn{1}{c}{Y_{\mathrm{d}}}\\ \multicolumn{1}{c}{Z_{\mathrm{d}}}\end{array}\right]=\left[\begin{array}{ccc}1\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& \cos {\alpha }_{\mathrm{R}}\hfill & \sin {\alpha }_{\mathrm{R}}\hfill \\ \multicolumn{1}{c}{0}& -\sin {\alpha }_{\mathrm{R}}\hfill & \cos {\alpha }_{\mathrm{R}}\hfill \end{array}\right]\left[\begin{array}{c}{X}_{\mathrm{P}}\hfill \\ \multicolumn{1}{c}{Y_{\mathrm{P}}}\\ \multicolumn{1}{c}{Z_{\mathrm{P}}}\end{array}\right]=T_{\mathrm{R}}\left[\begin{array}{c}{X}_{\mathrm{P}}\hfill \\ \multicolumn{1}{c}{Y_{\mathrm{P}}}\\ \multicolumn{1}{c}{Z_{\mathrm{P}}}\end{array}\right]$$

(21)

Combine the above three equations it could be obtain:

$$\begin{array}{c}\multicolumn{1}{c}{\left[\begin{array}{c}{X}_{\mathrm{d}}\hfill \\ \multicolumn{1}{c}{Y_{\mathrm{d}}}\\ \multicolumn{1}{c}{Z_{\mathrm{d}}}\end{array}\right]=T_{\mathrm{H}}{T}_{\mathrm{R}}{T}_{\mathrm{P}}\left[\begin{array}{c}X\hfill \\ \multicolumn{1}{c}{Y}\\ \multicolumn{1}{c}{Z}\end{array}\right]}\\ \multicolumn{1}{c}{=[\begin{array}{c}\cos {\alpha }_{\mathrm{H}}\cos {\alpha }_{\mathrm{P}}\hfill \\ \multicolumn{1}{c}{-\sin {\alpha }_{\mathrm{H}}\cos {\alpha }_{\mathrm{R}}-\cos {\alpha }_{\mathrm{H}}\sin {\alpha }_{\mathrm{R}}\sin {\alpha }_{\mathrm{P}}}\\ \multicolumn{1}{c}{-\cos {\alpha }_{\mathrm{H}}\sin {\alpha }_{\mathrm{P}}\sin {\alpha }_{\mathrm{R}}+\sin {\alpha }_{\mathrm{R}}\sin {\alpha }_{\mathrm{H}}}\end{array}}\\ \begin{array}{c}\begin{array}{cc}\cos {\alpha }_{\mathrm{P}}\sin {\alpha }_{\mathrm{H}}\hfill & \sin {\alpha }_{\mathrm{P}}\hfill \end{array}\hfill \\ \multicolumn{1}{c}{\begin{array}{cc}\cos {\alpha }_{\mathrm{H}}\cos {\alpha }_{\mathrm{R}}-\sin {\alpha }_{\mathrm{R}}\sin {\alpha }_{\mathrm{H}}\sin {\alpha }_{\mathrm{P}}\hfill & \cos {\alpha }_{\mathrm{P}}\sin {\alpha }_{\mathrm{R}}\hfill \end{array}}\\ \multicolumn{1}{c}{\begin{array}{cc}-\cos {\alpha }_{\mathrm{H}}\sin {\alpha }_{\mathrm{R}}-\sin {\alpha }_{\mathrm{H}}\sin {\alpha }_{\mathrm{P}}\cos {\alpha }_{\mathrm{R}}\hfill & \cos {\alpha }_{\mathrm{P}}\cos {\alpha }_{\mathrm{R}}\hfill \end{array}}\end{array}]\left[\begin{array}{c}X\hfill \\ \multicolumn{1}{c}{Y}\\ \multicolumn{1}{c}{Z}\end{array}\right]\hfill \end{array}$$

(22)

The coordinate points $(X_{\mathrm{d}},Y_{\mathrm{d}},Z_{\mathrm{d}})$ in the optical axis coordinate of submarine are converted to azimuth and pitch angles $(A_{\mathrm{d}},E_{\mathrm{d}})$ , which can be expressed as:

$$\{\begin{array}{c}{A}_d=\arctan \left\{\frac{\cos E_0[\cos {\alpha }_R\sin (A_0-{\alpha }_H)-\cos (A_0-{\alpha }_H)\sin {\alpha }_p\sin {\alpha }_R]+\sin E_0\sin {\alpha }_R\cos {\alpha }_P}{\sin E_0\sin {\alpha }_P+\cos E_0\cos {\alpha }_P\cos (A_0-{\alpha }_H)}\right\}\\ E_d=\arcsin \left\{\cos {\alpha }_P\cos {\alpha }_{\mathrm{R}}\sin E_0-\cos E_0[\sin {\alpha }_P\sin (A_0-{\alpha }_H)+\sin {\alpha }_P\cos {\alpha }_{\mathrm{R}}\cos (A_0-{\alpha }_H)]\right\}\end{array}$$

(23)

Where the $(A_0,E_0)$ are the initial azimuth and pitch angles.

Submarine torque and flow field characteristics

In this paper, the evolution of torque coefficients are discussed and used to reflect the disturbance characteristics of submarines. The total pressure difference torque coefficient and viscosity torque coefficient in three directions are defined as $M_{\mathrm{total}p}=2F_{p}l/(\rho U_{\infty }^{2}S)$ , $M_{\mathrm{total}\upsilon }=2F_{\upsilon }l/(\rho U_{\infty }^{2}S)$ , respectively. Where the S = 0.3613 denotes the non-dimensional surface area of the submarine.

Figure 7 displays that the $M_{\mathrm{total}p}$ fluctuated within 0.0011–0.0015, and $M_{\mathrm{total}\upsilon }$ fluctuated within a narrow range 2.646 × 10⁻⁵ – 2.654 × 10⁻⁵. The torque coefficient produced by the pressure is about 50 times than that of viscosity force. Therefore, the viscous force disturbance has been ignored, leaving the pressure to be discussed.

Figure 7.

Time histories of the total pressure $M_{\mathrm{total}p}$ and viscous $M_{\mathrm{total}\upsilon }$ torque coefficients of submarine.

Eventually, the torque’s influence on the photoelectric mast and tracking device will be reflected in azimuthal and pitch angles, whose strength are mainly related with the submarine mass and the reference coordinate system. The transformation between the torque and disturbance angles can be regarded as a linear relationship. In numerical analysis, we define $k=\overline{M_{\mathrm{total}p}}+\overline{M_{\mathrm{total}\upsilon }}$ as the proportionality coefficient that $\overline{M_{\mathrm{total}p}}$ , $\overline{M_{\mathrm{total}\upsilon }}$ are the average value of the total pressure and the total viscous force, respectively. The disturbance coefficients $C_{M_{p}x}$ , $C_{M_{p}y}$ , $C_{M_{p}z}$ in the x, y, and z directions can be obtained by the components of differential pressure torque dividing by k, respectively.

Once considered as an unrestraint rigid body, submarine could move and roll freely with six degrees of freedom (translational motion three, rotational motion three). However, in this simulation, we suppose that pressure does not change the submarine attitude. Usually, pressure pulsation is defined as the pressure fluctuate around its mean value $p\prime =p-\overline{p}$ . $C_{M_{p}i}(5)$ , i = x, y, z at the 5 s are taken as initial values in each direction and we defined relative pressure torque pulsation coefficient as $C_{M_{p}i}(t)$ = $C_{M_{p}i}(t)$ – $C_{M_{p}i}(5)$ , i = x, y, z, which are shown in Figure 8.

Figure 8.

The pressure torque pulsation coefficient curves from 5 s.

In this discussion, these rotation tendencies or the angles around the earth coordinate system are totally caused by the fluid dynamic moment’s waves. However, during the whole simulation process, the submarine keeps a same attitude. Combined with the momentum theorem, the moments can reflect an object’s rotation motion around its rotation center. Thus, the waves of ${\alpha }_{\mathrm{H}}$ , ${\alpha }_{\mathrm{p}}$ and ${\alpha }_{\mathrm{R}}$ can be replaced by $C_{M_{p}x}$ , $C_{M_{p}y}$ , $C_{M_{p}z}$ in Figure 8 and substituted into the coordinate transformation equation (15). In this simulation, we set the non-dimensional initial azimuth angle A₀ = 0.4 and the pitch angle E₀ = 0.2, the time histories of A_d (t) and E_d (t) are solved (as shown in Figure 9).

Figure 9.

Time histories of (A_d(t), E_d(t)), (A_d0(t), E_d0(t)), (A_d1(t), E_d1(t)).

To obtain much more deeply investigate on the effects of different amplitude pressure torque pulsation introduced by different scale vortices, the average filter methods with two differ parameters have been used to convert the (A_d(t), E_d(t)) curves shown in Figure 9 into two different disturbance curves (A_d0(t), E_d0(t)) and (A_d1(t), E_d1(t)). From the perspective of turbulence analysis, mean-filtering is to make the small scale random fluid movement causing the pressure pulsation dissipate and mix together with the surrounded larger scale vortices. Figure 9 is the curves of (A_d(t), E_d(t)) without mean-filtering. (A_d1(t), E_d1(t)) have been processed by high frequency filter, and (A_d0(t), E_d0(t)) are with low-pass filter. As illustrated in Figure 9 the fluctuation amplitudes of (A_d0(t), E_d0(t)) are approximately equal to 0.1, with their periods being about 0–15 s. (A_d1(t), E_d1(t)) wave around the curves of (A_d0(t), E_d0(t)) and it also contains small period’s oscillations about 2–3 s. The original curves (A_d(t), E_d(t)) contain both the fluctuation characteristics of (A_d1(t), E_d1(t)) and also characterized by the random pulsation traces at small amplitude and higher frequency.

Mathematical connotation of vortex strength is the mean variance of the flow field velocities curl in three directions, which reflects the rotation trend of space fluid. The interaction between the rotational fluid and the submarine surface viscosity boundary causes the pulsation of pressure torques. In order to get a deeper investigate on the time histories of pressure torque pulsation coefficients shown in Figure 9. Vortices evolution characteristics from t = 30–41 s have been listed in Figure 10.

Figure 10.

Vortices’ evolution characteristics from t = 30–41 s.

Figure 10 displayed the whole development process of vortex structures from 30 s to 41 s. The structures are visualized as iso-surfaces of vorticity, where the color indicates the local pressure. It can be seen that a large-scale A-kind vortex denotes as A_vor1 begin to shed from the upside of the conning tower. The A_vor1 vortex could keep the original shape and move toward downstream during the following 3 s 31–33 s. Since 33 s, A_vor1 moves out of the submarine body and dissipated into several medium-scale vortices. With time going on, the medium-scale vortices developed from the A_vor1 vortex are gradually transformed into a variety of small scale vortices and spread into the surrounding flow field.

By comparing the flow field at 32 s and 31 s in Figure 10, multiple small B_vors1 sub-vortices appear in the purple dotted line frame within 1 s. These B_vors1 sub-vortices are followed the A_vor1 vortex closely and move to the downstream. The kind of B_vors1 sub-vortices contain lower energy, so that most of them dissipate rapidly and develop into several interlaced smaller-scale vortices when they move out of the submarine. From the perspectives of time, the series of small-scale B_vors1 sub-vortices produce about 2–3 Hz waves, which have reflected on the changes of pressure torque pulsation coefficient, as shown in Figure 9. Comparing the filtered A_d1(t) and E_d1(t) with the unfiltered curves A_d(t) and E_d(t), it is detected that the unfiltered curves contain more small amplitude waves at frequcecies about 2–3 Hz, which are aroused from the B_vors1 sub-vortices.

Moreover, the generation mechanism for the B and A-kinds vortices is different. A-kind vortices are generated from the upside of the conning tower, they are of much kinetic energy, lower-frequency and larger-amplitude. A_vor1 is followed by B_vors1, showing in the red dotted-line box from 32 s to 36 s. Similarly, A_vor2 appears at 38 s about 7 s later than A_vor1, which will then continue and repeat the same evolutionary process. A-kind vortices are deducted at the time of 9 s, 17 s, 30 s and 38 s, for the whole process, the corresponding fluctuation peaks and dips also can be detected one after the other from A_d0(t) and E_d0(t) curves in Figure 9.

Furthermore, it could be found from 39 s to 41 s, under the flow’s action, the wake vortex track swing at a large amplitude and longer period. That caused by the submarine hull body whose length is much longer than another appendage. Its fluctuation periods are about 10–20 s ranged from 5–17 s, 17–35 s, and 35–48 s that are reflected by the heaving billows of A_d0(t) and E_d0(t) curves in Figure 9. Compared A_d1(t) and E_d1(t), it is detected that the unfiltered curves A_d(t) and E_d(t) contain much more small amplitude waves, about 2–3 Hz, which are bring about by the B_vors1 sub-vortices.

In general, based on the frequency characters, the azimuth and pitch angle curves A_d0(t) and E_d0(t) can be divided into three forms: the high frequency random small perturbations accompany with B kind vortices about 2–3 Hz, the middle frequency perturbations accompany with A kind vortices about 1/4–1/6 Hz, and the low frequency perturbations about 1/15–1/20 Hz. According to the flow field analysis, the fluctuation of middle/low frequencies are mainly resulted from two parts of submarine that is the hull and appendage. High frequency random small perturbations are generated from vortices motion, which cause fluid’s small pulsation, fluid mix or energy exchange among different fluid layers.

Mathematical model of servo system

Structure and transfer function of the DC torque motor

The coarse channel is mainly composed of DC torque motor (DTM) and rotational platform as Figure 11 showing. In this paper, the math model of DC torque motor is deduced.

Figure 11.

Experimental equipment of tracking and aiming system.

In order to simplify the research, it can be assumed that the electrical parameters and transfer functions of the azimuth torque motor and its rotational inertia are strictly consistent with the electrical parameters and transfer function of the pitch torque motor. Current fluctuation can be ignored and the DC torque motor can be estimated to be a linear device. Thus, based on the principle of motor the voltage balance equation of armature circuit can be obtained:

$$U_a=R_e{i}_a+L_e\frac{d{i}_a}{\mathrm{dt}}+e_b$$

(24)

Where $R_e,i_a,L_e,U_a$ denote the armature resistance, motor-driven current, armature inductance and input control voltage in the armature circuit, respectively. The counter EMF denotes as $e_b=K_e\stackrel{·}{\theta },$ $K_e$ is the negative electromotive force coefficient of the motor, $\stackrel{·}{\theta }$ is the angle velocity of the motor. The electromagnetic torque of the motor is $M_{\mathrm{d}}=K_m{i}_a$ , $K_m$ is electromagnetic torque coefficient or torque sensitivity coefficient. Torque balance equation of the motor is. ²¹

$$M_{\mathrm{d}}-M_{\mathrm{f}}=J\stackrel{··}{\theta }$$

(25)

Where J is the load equivalent torque of inertia, $M_{\mathrm{f}}$ is the disturbance torque. $\stackrel{··}{\theta }$ is the angle acceleration of the DC motor. The transfer function of torque motor can be written as:

$$\frac{\stackrel{·}{\theta }(\mathrm{s})}{U_a(\mathrm{s})}=K_m\times \frac{1}{L_e\mathrm{s}+R_e+K_e{K}_m/J\mathrm{s}}\times \frac{1}{J\mathrm{s}}$$

(26)

$K_e{K}_m\to 0$ so the above equation (26) could be simplified as:

$$\frac{\stackrel{·}{\theta }(\mathrm{s})}{U_a(\mathrm{s})}=K_m\times \frac{1}{R_e+L_e\mathrm{s}}\times \frac{1}{J\mathrm{s}}$$

(27)

For this model, the peak torque voltage is 27V, the maximum rotational velocity $\stackrel{·}{\theta }=2000r/\min =209.4\mathrm{rad}/s$ , $K_m=0.1N·\mathrm{m}/\mathrm{A}$ , $J=4.4\times {10}^{-5}$ kg $·{\mathrm{m}}^2$ , $R_{\mathrm{e}}=7.5\Omega$ , $L_{\mathrm{e}}=2.07\mathrm{mH}$ .

$$\begin{array}{c}\multicolumn{1}{c}{M_{\mathrm{f}}=0;\frac{1}{R_e+L_e\mathrm{s}}=\frac{1}{7.5+0.00207\mathrm{s}}\approx \frac{1}{7.5};}\\ \multicolumn{1}{c}{\frac{1}{J\mathrm{s}}=\frac{1}{4.4\times {10}^{-5}\mathrm{s}};}\end{array}$$

(28)

In order to improve the response speed, the speed loop adopts PI controller and the position loop adopts PID controller, as shown in Figure 12. For the system, this steady-state error under a unit step response without distribution is about 4.15×10^-5.

$$G(\mathrm{s})=K_p+\frac{K_p}{T_i}·\frac{1}{\mathrm{s}}+K_p\tau \mathrm{s}$$

(29)

Figure 12.

Angle error transfer function of the coarse channel.

$G(\mathrm{s})$ is the common format of the PID or PI controllers $G_1(\mathrm{s})$ and $G_2(\mathrm{s})$ . Where $K_p$ the proportionality factor, $T_i$ integration time constant, $\tau$ derivative time constant. The PID controller $G_1(\mathrm{s})$ , the $K_p=275$ , $T_i=137.5$ , $\tau =0.00272$ . And in PI controller $G_2(\mathrm{s})$ the $K_p=100$ , $T_i=16.7$ , $\tau =0$ .

Structure parameters and transfer function of fine channel FSM

The physical swing magnitude for each axis of the FSM in fine channel is ±12.5 mrad and the closed loop repeat positioning accuracy is 5 μrad. The FSM’s structure is shown in Figures 1 and 11, which primarily composed of a cuboid body, with four moveable coplanar fulcrum points on its head, onto which the mirror is fixed. Through the FSM controller, the length of inside piezoelectric ceramic stacks could be adjusted with the control voltages. Meanwhile, the deflection and pitch angles of the mirror fixed onto the four moveable fulcrum points can be accurately controlled. A spectrum analyzer is used to measure the frequency response curve of FSM in the closed-loop state, and the closed-loop transfer function could be deduced as^22,23:

$$G_{\mathrm{FSM}}(\mathrm{s})=\frac{1}{(0.0036\mathrm{s}+1)({2.53}^{-6}{\mathrm{s}}^2+0.0022\mathrm{s}+1)}$$

(30)

For the fine channel FSM, PID control method is used, its steady-state error under a unit step response without distribution is about 8.1 × 10⁻¹⁵. This control precise can only be reached in theory and idea conditions.

Transfer function of compound axis tracking system

Figure 13 is the MATLAB Simulink simulation diagram of the transfer function for this compound axis tracking system. The output of coarse tracking angle errors are then sent to the fine tracking PID controller as its input signal. The fine channel PID controller has the $K_p=3.65$ , $T_i=0.01$ , $\tau =0.000312$ .

Figure 13.

Transfer function of the coarse and fine composite axis.

In order to analyze the system’s robustness, in this numerical simulation, the unit step input signal composed with a continual random disturbance is investigated. For this tracking and aiming system, the input signal contains continuous random input errors of roughly 5% total amplitude and 1 kHz wave frequency. Then the input signal is send to the whole system demonstrated on Figure 13 in MATLAB and Figure 14 is the output response curves obtained from the scope 1 and 2 on Figure 13. The corresponding rise edge of the coarse channel is about 0.15 s. Due to the limitation of the system bandwidth, there is not effective error compensation result was found for the 1 kHz system random disturbance input. Thus, the continual random disturbance has no obvious influence on the stability of the system.

Figure 14.

Unit step response added a random disturbance input: (a) coarse tracking output, Scope1 and (b) fine tracking error, Scope2.

Influence of flow field disturbance on steady error

The time histories of $(A_{\mathrm{d}}(t),E_{\mathrm{d}}(t))$ , $(A_{\mathrm{d}0}(t),E_{\mathrm{d}0}(t))$ and $(A_{\mathrm{d}1}(t),E_{\mathrm{d}1}(t))$ in Figure 9 are taken as the disturbance input of this tracking system during the numerical simulations. The final tracking errors of azimuth and pitch angle are obtained by MATLAB Simulink simulation.

Figure 15(a) to (c) show the tracking error tracks of azimuth and pitch angles through the coarse tracking, which display that under three scale’s disturbances the output errors of azimuth angle $\mathrm{Er}{r}_{A_{\mathrm{d}}\mathrm{-}\mathrm{DTM}}$ , $\mathrm{Er}{r}_{A_{\mathrm{d}1}\mathrm{-}\mathrm{DTM}}$ , $\mathrm{Er}{r}_{A_{\mathrm{d}0}\mathrm{-}\mathrm{DTM}}$ wave within $\pm 3.5\times {10}^{-4}$ , $\pm 1.8\times {10}^{-4}$ , $\pm 5\times {10}^{-5}$ , respectively. The output errors of pitch angle $\mathrm{Er}{r}_{E_{\mathrm{d}}\mathrm{-}\mathrm{DTM}}$ , $\mathrm{Er}{r}_{E_{\mathrm{d}1}\mathrm{-}\mathrm{DTM}}$ , $\mathrm{Er}{r}_{E_{\mathrm{d}0}\mathrm{-}\mathrm{DTM}}$ wave within $\pm 5\times {10}^{-4}$ , $\pm 2\times {10}^{-4}$ , $\pm 5.2\times {10}^{-5}$ , respectively. Figure 16(a) to (c) show the azimuth and pitch angle tracking error tracks for compound axis or FSM/fine channel. Table 1 describes the different output errors for the different channels and inputs, which are simulated by MATLAB Simulink. We denote ${\sigma }_0=\sqrt{Er{r}_{A_{\mathrm{d}}-\mathrm{DTM}}\times Er{r}_{E_{\mathrm{d}}-\mathrm{DTM}}}$ and ${\sigma }_1=\sqrt{Er{r}_{A_{\mathrm{d}}-\mathrm{FSM}}\times Er{r}_{E_{\mathrm{d}}-\mathrm{FSM}}}$ as the error square deviations to indicate the error area scale.

Figure 15.

Three track distributions for the coarse tracking errors (a) Without mean-filter; (b) With high frequency filter; (c) With low-pass filter

Figure 16.

Three track distributions for the fine tracking errors (a) Without mean-filter; (b) With high frequency filter; (c) With low-pass filter.

Table 1.

Error distributions of coarse and FSM/fine channel.

		Coarse channel			FSM/fine channel
	Input error	$\mathrm{Er}{r}_{A_{\mathrm{d}}\mathrm{-}\mathrm{DTM}}$ $\times {10}^{-4}$	$\mathrm{Er}{r}_{E_{\mathrm{d}}\mathrm{-}\mathrm{DTM}}$ $\times {10}^{-4}$	${\sigma }_0$ $\times {10}^{-4}$	$\mathrm{Er}{r}_{A_{\mathrm{d}}\mathrm{-}\mathrm{FSM}}$ $\times {10}^{-5}$	$\mathrm{Er}{r}_{E_{\mathrm{d}}\mathrm{-}\mathrm{FSM}}$ $\times {10}^{-5}$	${\sigma }_1$ $\times {10}^{-5}$
Case 1	$(A_{\mathrm{d}}(t),E_{\mathrm{d}}(t))$	$\pm 3.5$	$\pm 5$	≅4.2	$\pm 1.8$	$\pm 2$	≅1.9
Case 2	$(A_{\mathrm{d}1}(t),E_{\mathrm{d}1}(t))$	$\pm 1.8$	$\pm 2$	≅1.9	$\pm 0.2$	$\pm 0.2$	≅0.2
Case 3	$(A_{\mathrm{d}0}(t),E_{\mathrm{d}0}(t))$	$\pm 0.5$	$\pm 0.52$	≅0.5	$\pm 0.03$	$\pm 0.02$	≅0.024

Firstly, for the coarse channel, Case 2 is the situation that filtered out high frequency pulsation caused by small scale high-frequency random perturbation vortices. Compared with the original input Case 1, the coarse channel accuracy is improved by about 2.21 times from 4.2 to 1.9. Moreover, compared with Case 2, Case 3 filters out the middle-frequency perturbation that could further improve the coarse channel accuracy to roughly 1.9/0.5 = 3.8 times.

Secondly, for fine channel, the corresponding improvement error times for composite tracking or FSM are about 1.9/0.2 = 9.5 and 0.2/0.024 = 8.3, respectively. This difference mainly relies on the narrow bandwidth of coarse channel. The distribution errors from medium scale vortices’ pressure pulsations have been significant improved by the coarse channel, which is about 3.8 times (Case 2 to Case 3) lager than 4.2/1.9 = 2.21 times (Case 1 to Case 2). However, for composite tracking, fine channel FSM supports a wider bandwidth appearing a better error compensation effect for random small scale, higher frequency pressure pulsation. Thus, the error enhancement time for high frequency from Case 1 to Case 2 is about 1.9/0.2 = 9.5 slightly larger than 0.2/0.024 = 8.3 from Case 2 to Case 3. But this comparison values for coarse channel are 2.21 and 3.8. It can be concluded that for a FSM offers 100–200 Hz bandwidth range, the error improvement capability is stronger for the introduction of small random perturbations with a frequency of 2–3 Hz or higher. During a new system’s design, more effort should be contributed on the selection of DC torque motor to catch the error at the middle frequencies.

Figure 17(a) to (c) show the azimuth and pitch angle tracking error density distributions for compound axis or FSM/fine channel tracking. It can be seen that the maximum density values of Figure 17(a) to (c) can reach to 1.6×10¹⁷, 8×10¹⁷, and 4.5×10¹⁹, respectively. Thus, for direction energy weapon, the control method in Figure 17(c) has a better control precision and laser destroying ability.

Figure 17.

Three track density distributions for the fine tracking errors (a) Without mean-filter; (b) With high frequency filter; (c) With low-pass filter.

In general, through the simulation analysis, we can further obtain the disturbance characteristics and generation mechanism of multi-scale vortices as well as the corresponding pressure pulsation. Furthermore, the error compensation effects of coarse/fine channels for multi-scale pressure pulsation have been deeply investigated, providing more predictive guidance for the design of control system in practice.

Conclusions

Submarine hydrodynamic characteristics of 6° yaw angle underwater have been investigated in this investigation. Based on the vortices’ evolution, the relationships between the multi-scale pressure pulsation and output photoelectric servo tracking error have been numerical analyzed with its main conclusions as follows:

1. For submarine torque disturbance in high Reynolds number flow field, the influence of pressure pulsation plays a dominant role than the viscous force.

2. The fluctuations of middle/low frequencies are mainly caused by the boundary layer separation and shedding from the fluid flow around the submarine hull and appendage. High frequency random small perturbations mainly generated from the fluid’s random interaction as well as the energy exchange among different flow layers.

3. Pressure pulsations formed by different scale’s turbulence dissipative vortices are of difference impact on the error output of the tracking and targeting system. The coarse channel has more excellent performance on suppressing the middle scale pressure pulsation. While the FSM/fine channel has the best control effects for the high frequency small-scale pressure pulsation, it contributes more than the coarse channel, especially in complex flow environments.

The research on submarine’s fluid dynamics, including the inherent influence laws from the vertical/horizontal rudders and the propellers under different angles of attack needs further study. In addition, more accurate numerical methods and turbulence models will be adopted to provide a better and detailed view of the flow fluid. Moreover, the generating chaotic phenomena for the turbulence and distributions would be investigated, which are usually generated by non-linearities occurring in the examined system. The chaotic theories can help us to discuss the propeller distribution as well as detect the fine shadows from the propeller damage.