Hydrodynamic Model of Diver–DPV Coupled Multi-Body and Its Underwater Cruising Numerical Simulation

Diver propulsion vehicles (hereinafter referred to as DPV) are a kind of small vehicle with underwater high-speed used by divers, who are able to grasp or ride on, and operate the volume switch to change the speed. Different from unmanned underwater vehicles (UUVs), the interference caused by diver’s posture changing is a unique problem. In this paper, a Diver–DPV multi-body coupling hydrodynamic model considering rigid body dynamics and fluid disturbance is established by analyzing the existing DPV related equipment. The numerical simulation of multi-body articulated motion is realized by using Star-CCM+ overlapping grid and DFBI 6-DOF body motion method. Five cases of DPVs underwater cruising in a straight-line when restraining diver movement is simulated, and five cases with free diver movement are simulated too. Finally, the influence of the diver’s posture changing on the cruising speed resistance is analyzed, and the motion equation including the disturbance is solved. The final conclusion is that, the disturbance is favorable at high speed, which can reduce the cruising resistance, and unfavorable at low speed, which increases the cruising resistance. The friction resistance Ff always accounts for the main part in all speed cases.


Introduction
Diver propulsion vehicle (hereinafter referred to as DPV) is a kind of small vehicle with underwater high-speed used for divers [1], who are able to grasp or ride on, and operate the volume switch to change the speed. The numerical simulation of DPV cruising in water is different from that of HOVs (human-occupied vehicles) and UUVs (unmanned underwater vehicles). That is, when DPV cruising, the diver's body will swing with the flow field, and its posture change will have a nonlinear disturbance to the surrounding flow field, which has a huge impact on the speed and stability of the DPV. Meanwhile, the DPV's hydrodynamic shape and propeller arrangement will also have an impact on the diver body's posture. There is also the problem of the propeller wake impacting the diver's body. In the traditional numerical simulation method, diver and DPV are regarded as a relatively fixed rigid body. While in the multi-body coupling model, the joints of diver body and DPV are regarded as a complex hinge connecting rod structure, which has active and driven motion, and the two have the influence of fluid disturbance and rigid hinge force. Using CFD (compulation fluid dynamics) method to simulate the underwater cruising of a multi-body coupling model, and finally predicting the performance and analyzing the posture of the diver are the keys to solving the design optimization problem of high-performance DPVs.
Since the last century, scholars from all over the world have begun to study the design theory of DPVs and other diving equipment. In 1987, Presswood, Clark G, Godshall, David and others [2] first proposed a set of propeller matching theory applicable to DPV, and applied to the optimization of Tekna DV-3X and MIL-UNIT MOD S-5100 DPV.
CFD method is widely used in the numerical simulation of DPV linear cruising. In 2017, M.R. Sadeghizadeh and B.Daranjam [3] studied the numerical simulation of the diver-DPV coupling model by CFD in detail, and put forward the DPV optimization theory of improving the rapidity and reducing the impact on the human face is proposed. The problem that part of diver body is exposed to the external flow field is considered in the modeling. Finally, the model experiment of towing tank is carried out and the accuracy of the numerical simulation is verified, and its CFD velocity cloud chart is presented in Figure 1. In particular, the k-ω SST Turbulence model is used to simulate the complex vortex region at the connection between the diver's body and DPV, and the mesh of the local prism layer is encrypted on the human face. ble to DPV, and applied to the optimization of Tekna DV-3X and MIL-UNIT MOD S-5100 DPV.
CFD method is widely used in the numerical simulation of DPV linear cruising. In 2017, M.R. Sadeghizadeh and B.Daranjam [3] studied the numerical simulation of the diver-DPV coupling model by CFD in detail, and put forward the DPV optimization theory of improving the rapidity and reducing the impact on the human face is proposed. The problem that part of diver body is exposed to the external flow field is considered in the modeling. Finally, the model experiment of towing tank is carried out and the accuracy of the numerical simulation is verified, and its CFD velocity cloud chart is presented in Figure 1. In particular, the k-ω SST Turbulence model is used to simulate the complex vortex region at the connection between the diver's body and DPV, and the mesh of the local prism layer is encrypted on the human face.
In the CFD method of DPV, the movable human body modeling is also the research focus. At the end of the 20th century, the body hydrodynamic modeling and CFD numerical simulation of swimmers and divers were began to take. In 1996, Bixler, B. S. and M. Schloder et al. [4] from University of Calgary firstly put forward CFD method to simulate the dynamic and posture model of swimmers, so as to predict, analyze, and improve the performance of swimmers. In terms of 3D action model of swimmers, Motomu Nakashima，Ken Datou, and Yasufumi Miura [5] proposed a method for the establishment of a numerical simulation model of swimming human body considering rigid body dynamics and nonlinear hydrodynamic interference in 2007. In this model, the human limb model rotates in the water, and its rotation angle and root bending moment are measured, and the splitting steps are shown in Figure 2. Although the error of the model is 10%, it can better reflect the unsteady fluctuation of experimental data. In addition, it also simulates the position of the human body that can generate the gliding lift to determine the tangent resistance coefficient. Finally, a complete six-step freestyle simulation example is given. The simulation speed is 7.5% lower than the actual swimming speed.
CFD method of swimming human is divided into grid method and grid free method. COHEN, R.C.Z. et al. [6] proposed a computational method for human swimming using smoothed particle hydrodynamics (SPH method) in the year of 2009. It points out that the SPH method can better solve the influence of turbulence, complex free surface motion (including gas splash and entrainment) and the rapid deformation of swimmers' geometric structure. The preliminary simulation of SPH traction for male and female swimmers at different speeds with fixed gliding posture is carried out. In the experiment, laser scanning is used to establish a more accurate three-dimensional model of human body, and In the CFD method of DPV, the movable human body modeling is also the research focus. At the end of the 20th century, the body hydrodynamic modeling and CFD numerical simulation of swimmers and divers were began to take. In 1996, Bixler, B. S. and M. Schloder et al. [4] from University of Calgary firstly put forward CFD method to simulate the dynamic and posture model of swimmers, so as to predict, analyze, and improve the performance of swimmers.
In terms of 3D action model of swimmers, Motomu Nakashima, Ken Datou, and Yasufumi Miura [5] proposed a method for the establishment of a numerical simulation model of swimming human body considering rigid body dynamics and nonlinear hydrodynamic interference in 2007. In this model, the human limb model rotates in the water, and its rotation angle and root bending moment are measured, and the splitting steps are shown in Figure 2. Although the error of the model is 10%, it can better reflect the unsteady fluctuation of experimental data. In addition, it also simulates the position of the human body that can generate the gliding lift to determine the tangent resistance coefficient. Finally, a complete six-step freestyle simulation example is given. The simulation speed is 7.5% lower than the actual swimming speed.
CFD method of swimming human is divided into grid method and grid free method. COHEN, R.C.Z. et al. [6] proposed a computational method for human swimming using smoothed particle hydrodynamics (SPH method) in the year of 2009. It points out that the SPH method can better solve the influence of turbulence, complex free surface motion (including gas splash and entrainment) and the rapid deformation of swimmers' geometric structure. The preliminary simulation of SPH traction for male and female swimmers at different speeds with fixed gliding posture is carried out. In the experiment, laser scanning is used to establish a more accurate three-dimensional model of human body, and the dynamic model of articulation is generated based on the surface mesh deformation of human bones. Finally, the underwater video shot of swimmers is captured to verify the model. the dynamic model of articulation is generated based on the surface mesh deformation of human bones. Finally, the underwater video shot of swimmers is captured to verify the model.

Figure 2.
Swimming dynamics modeling and dynamic process simulation proposed by Nakashima et al. [5].
In the aspect of grid method, M.R. Sadeghizadeh and B. Daranjam et al. [7] conducted numerical simulation of three-dimensional two-phase turbulent flow on the basis of fluid volume method (VOF) of FLUENT software and tetrahedral mesh in 2012, so as to evaluate the effective power required by high-speed underwater swimmers. The speed of the swimmer model is increased to 8 m/s by using the towing tank experiment, and the results are consistent with CFD.
In 2010, Zaidi H. et al. [8] analyzed different types of turbulence models in CFD of swimmers in order to better predict resistance in fluid environment. The test of two turbulence models shows that the standard k-ω model can accurately predict the resistance and capture the vortex structure formed by the back and buttocks of underwater swimmers, while the standard K-ω model underestimates the value of resistance. Finally, the flow visualization experiment in insep (National Institute of Sports) verifies the rationality of vortex structure under K-ω model.
In 2012, Ramos et al. [9] used CFD method to analyze the influence of water depth on the resistance of swimmers during underwater gliding, and calculated the resistance coefficient (increment is 0.5 m/s) in the speed range of 1.5-2.5 m/s. It is found that the resistance decreases with the increase of depth when the center line of swimmer's modeling is in different depths of 0-0.75 m (increment is 0.25 m), while when the depth is greater than 0.75 m, the resistance value is almost the same stay the same.
Compared with swimmers, the underwater motion of diver wearing scuba diving equipment is more complex. In 2016, Darko Valenko et al. [10] firstly established the dynamic model and simulation method of scuba diver and buoyancy control device (BCD), simulated the complete heave motion process, and compared the simulation or results with experimental results.
The problem of divers driving DPV in the water is essentially the mutual disturbance of multiple underwater objects, such as the problem of the submarine running near the wall, the problem of two vehicles meeting or exceeding in narrow waters, the mutual interference between the legs of the jacket platform, and so on. Previous studies on the two- In the aspect of grid method, M.R. Sadeghizadeh and B. Daranjam et al. [7] conducted numerical simulation of three-dimensional two-phase turbulent flow on the basis of fluid volume method (VOF) of FLUENT software and tetrahedral mesh in 2012, so as to evaluate the effective power required by high-speed underwater swimmers. The speed of the swimmer model is increased to 8 m/s by using the towing tank experiment, and the results are consistent with CFD.
In 2010, Zaidi H. et al. [8] analyzed different types of turbulence models in CFD of swimmers in order to better predict resistance in fluid environment. The test of two turbulence models shows that the standard k-ω model can accurately predict the resistance and capture the vortex structure formed by the back and buttocks of underwater swimmers, while the standard K-ω model underestimates the value of resistance. Finally, the flow visualization experiment in insep (National Institute of Sports) verifies the rationality of vortex structure under K-ω model.
In 2012, Ramos et al. [9] used CFD method to analyze the influence of water depth on the resistance of swimmers during underwater gliding, and calculated the resistance coefficient (increment is 0.5 m/s) in the speed range of 1.5-2.5 m/s. It is found that the resistance decreases with the increase of depth when the center line of swimmer's modeling is in different depths of 0-0.75 m (increment is 0.25 m), while when the depth is greater than 0.75 m, the resistance value is almost the same stay the same.
Compared with swimmers, the underwater motion of diver wearing scuba diving equipment is more complex. In 2016, Darko Valenko et al. [10] firstly established the dynamic model and simulation method of scuba diver and buoyancy control device (BCD), simulated the complete heave motion process, and compared the simulation or results with experimental results.
The problem of divers driving DPV in the water is essentially the mutual disturbance of multiple underwater objects, such as the problem of the submarine running near the wall, the problem of two vehicles meeting or exceeding in narrow waters, the mutual interference between the legs of the jacket platform, and so on. Previous studies on the twodimensional analytical solutions of multi-body perturbation problems under the condition of potential flow have been carried out in detail.
In 2010, A.A. Tchieu, D. Crowdy, et al. [11] studied the interaction between two arbitrary bodies in two-dimensional non viscous fluid, and gave the linear velocity and angular velocity of the object. For the more complex problem with divers, it can be regarded as the problem of underwater multi-body articulated motion, which has a wider application in the field of numerical simulation and experimental research of bionic fish robots. In 2018, Lijun Li, Gen Li, et al. [12] carried out numerical simulation and motion control simulation on the multi fin propulsion of puffer fish. Ida Louise G. Borlaug and Kristin Y. Pettersen proposed a generalized over twist algorithm in 2019 [13,14], which is used to track the angle, position, and direction of the base joint of the 6-DOF of AIAUV (a biotic snake AUV, shown in Figure 3).  [11] studied the interaction between two arbitrary bodies in two-dimensional non viscous fluid, and gave the linear velocity and angular velocity of the object. For the more complex problem with divers, it can be regarded as the problem of underwater multi-body articulated motion, which has a wider application in the field of numerical simulation and experimental research of bionic fish robots. In 2018, Lijun Li, Gen Li, et al. [12] carried out numerical simulation and motion control simulation on the multi fin propulsion of puffer fish. Ida Louise G. Borlaug and Kristin Y. Pettersen proposed a generalized over twist algorithm in 2019 [13,14], which is used to track the angle, position, and direction of the base joint of the 6-DOF of AIAUV (a biotic snake AUV, shown in Figure 3). In the CFD method of multi joint biotic fish robots, the main difficulty lies in using more accurate moving grid technology and grid deformation technology to discretize the region, and using efficient solver to make the convergence speed catch up with the movement and deformation speed. Qing Xiao, Ke sun, Hao Liu, et al. [15] have studied the hydrodynamic performance of the wave NACA0012 plate near wake cylinder, which can simulate the multi-body disturbance of the rigid AUV shell to the BCF tail fin. The nephogram of vorticity is shown in Figure 4  Ruoxin Li, Qing Xiao, Yuan Liu, Jianxin, et al. [16] proposed a multi-body dynamics algorithm in 2016, which combines the incompressible fluid flow with the bionic multibody dynamics system. The tool has shown its powerful ability to solve a variety of biomechanical fish swimming problems, including self-propelled multi DOF with rigid undulator.
Yang Luo, Qing Xiao, Guangyu Shi, et al. proposed the multi body dynamic algorithm advanced code in 2019 [17] to carry out numerical simulation on the robot fish equipped with soft bionic fins. It is a fluid structure interaction (FSI) solver, which carries out finite element calculation based on CalculiX software. In the CFD method of multi joint biotic fish robots, the main difficulty lies in using more accurate moving grid technology and grid deformation technology to discretize the region, and using efficient solver to make the convergence speed catch up with the movement and deformation speed. Qing Xiao, Ke sun, Hao Liu, et al. [15] have studied the hydrodynamic performance of the wave NACA0012 plate near wake cylinder, which can simulate the multi-body disturbance of the rigid AUV shell to the BCF tail fin. The nephogram of vorticity is shown in Figure 4. dimensional analytical solutions of multi-body perturbation problems under the condi tion of potential flow have been carried out in detail.
In 2010, A.A. Tchieu, D. Crowdy, et al. [11] studied the interaction between two arbi trary bodies in two-dimensional non viscous fluid, and gave the linear velocity and angu lar velocity of the object. For the more complex problem with divers, it can be regarded a the problem of underwater multi-body articulated motion, which has a wider application in the field of numerical simulation and experimental research of bionic fish robots. I 2018, Lijun Li, Gen Li, et al. [12] carried out numerical simulation and motion control sim ulation on the multi fin propulsion of puffer fish. Ida Louise G. Borlaug and Kristin Y Pettersen proposed a generalized over twist algorithm in 2019 [13,14], which is used t track the angle, position, and direction of the base joint of the 6-DOF of AIAUV (a bioti snake AUV, shown in Figure 3). In the CFD method of multi joint biotic fish robots, the main difficulty lies in usin more accurate moving grid technology and grid deformation technology to discretize th region, and using efficient solver to make the convergence speed catch up with the move ment and deformation speed. Qing Xiao, Ke sun, Hao Liu, et al. [15] have studied th hydrodynamic performance of the wave NACA0012 plate near wake cylinder, which can simulate the multi-body disturbance of the rigid AUV shell to the BCF tail fin. The neph ogram of vorticity is shown in Figure 4  Ruoxin Li, Qing Xiao, Yuan Liu, Jianxin, et al. [16] proposed a multi-body dynamic algorithm in 2016, which combines the incompressible fluid flow with the bionic multi body dynamics system. The tool has shown its powerful ability to solve a variety of bio mechanical fish swimming problems, including self-propelled multi DOF with rigid un dulator.
Yang Luo, Qing Xiao, Guangyu Shi, et al. proposed the multi body dynamic algo rithm advanced code in 2019 [17] to carry out numerical simulation on the robot fish equipped with soft bionic fins. It is a fluid structure interaction (FSI) solver, which carrie out finite element calculation based on CalculiX software. Ruoxin Li, Qing Xiao, Yuan Liu, Jianxin, et al. [16] proposed a multi-body dynamics algorithm in 2016, which combines the incompressible fluid flow with the bionic multibody dynamics system. The tool has shown its powerful ability to solve a variety of biomechanical fish swimming problems, including self-propelled multi DOF with rigid undulator.
Yang Luo, Qing Xiao, Guangyu Shi, et al. proposed the multi body dynamic algorithm advanced code in 2019 [17] to carry out numerical simulation on the robot fish equipped with soft bionic fins. It is a fluid structure interaction (FSI) solver, which carries out finite element calculation based on CalculiX software.
In this paper, the flexible deformation of divers' fins is ignored in numerical simulation, so the function of Star-CCM+ software can be used to realize the numerical simulation of underwater multi-body articulated motion. There has been accurate CFD methods for the underwater motion of human body, which splits the swimming posture, simulates and integrates each step separately, and proposes the grid model and turbulence model with the most accurate calculation results. However, for the DPV underwater cruising, predecessors have fixed the diver body with DPV rigidly in the model, without active and passive motion of human joints. That is to say, the influence of nonlinear disturbance of multi rigid body is ignored. On this basis, the main research contents are as follows: (1) A diver-DPV coupled model with changeable posture is established.
(2) The cruising resistance of rigid connection immobility model is solved by two CFD methods of overlapping grids and DFBI multibody motion. (3) The influence of DPV resistance and movable human driving posture disturbance on the underwater cruising is studied.

Numerical Methodology
In this paper, Star-CCM+ software is used to simulate the nonlinear motion of underwater multibody. The following laws are satisfied for the Solver: conservation of mass, conservation of momentum, and conservation of energy. If the flow state is turbulent, the system also needs to follow the turbulent transport equation. It is defined as: Re = ρUL/µ = UL/γ, where U is the velocity of simulated free flow, and γ is the kinematic viscosity of fluid. According to the minimum length of the smallest part in this study, the minimum Reynolds number Re min occurs on the smallest appendage, and Re min is about 44005, which is shows that the flow field is turbulent, and the Reynolds number of the whole model's cruising speed Re is 6816272, which belongs to high Reynolds number.
In this paper, K-ω SST turbulence model is adopted, which is a two-way eddy viscosity model. The formula used in the interior of the boundary layer makes the model can be directly used to the wall through the viscous sublayer. Therefore, K-ω SST model can be used as turbulence model with low Re number and without any additional damping function. The SST formula is also converted to K-ω method in free flow, thus avoiding the common problem that the model is too sensitive to the turbulent characteristics of the inlet free flow. The K-ω SST model has good performance in adverse pressure gradient and separation flow [18].
In this paper, K-ω SST turbulence model is used. It meets the following requirements: Continuity equation: where P is the preesure, t ij is the stress tensor, τ t ij is Reynolds stress term. u i , u j is the fluctuating velocity component in the i, j direction fluid density, µ is the fluid dynamic viscosity, and ρ is the fluid density.
Because in turbulence, the actual velocity of each particle is the sum of instantaneous and fluctuating components. The last term ρu i u j , the Reynolds stress, should be obtained using a suitable turbulence model.
Predecessors have proved which model is more effective than other models with the same conditions. For example, Zaidi et al. (2010) proposed that k-ω SST model is more accurate for predicting human underwater swimming resistance. According to their model, k and ω are where the parameter "ω" is the energy dissipation rate (or specific dissipation rate) per unit volume, and "k" is the turbulent kinetic energy per unit mass. Γ k and Γ w are the effective diffusivities of k and ω. G k and G w are the turbulence results of k and ω. Y k and Y ω are the turbulent dissipations of k and ω. S k and S ω are source terms of k and ω. Therefore, the resistance of DPV can be obtained by this method, and the propulsion power P p related to the required cruising speed V can be calculated as where F t is real time resistance of hull.
On the other hand, if the propulsion power is electric, other parameters will come into question, such as propulsion efficiency η − p (composed of propeller and motor efficiency) and the total weight of diver and DPV. Under this condition, the input power P i is defined as the formula

6-DOF Body Splitting and Model Building
Firstly, the hydrodynamic shape of DPV is modeled. Whether DPV or AUV, its hydrodynamic shape directly affects the resistance and posture of underwater cruising. In this paper, a streamlined DPV with pump-jet propulsion is designed as the research object, which is referring to the Rotinor®Blackshadow 730 DPV [19]. Its parameters are in Table 1. The model abandons the exposed propeller layout of the normal DPV and adopts pump jet propulsion [20]. Which is, the inflow enters from the bottom-front openings, accelerates through the propeller in the duct, and ejects at a high velocity. In this model, when calculating the resistance of DPV, the model of propeller only retains the hub part.
In this model, the posture of diver driving DPV is considered as drag-and-drop. Which is, the diver navigates by holding the operating lever of the DPV with both hands, and uses torque to change the head direction of the DPV to carry out oblique navigation and steering movement like typical SUEX DPV [21].
This model assumes that the X-direction motion of DPV is locked, and a fixed velocity X-direction flow is used to simulate cruising. The local coordinate system was established with the centroid point P 0 of the hip part of diver. At this time, the typical cruising posture and force situation of diver and DPV are shown in Figure 5. The equation of motion is where a x and a y is the acceleration in X, Y directions.
with the centroid point P0 of the hip part of diver. At this time, the typical cruising posture and force situation of diver and DPV are shown in Figure 5. The equation of motion is where ax and ay is the acceleration in X, Y directions. For the convenience of calculation, the diver body posture at design speed is given as T0 posture. According to the Professional Association of Diving Instructors (PADI) DPV driving course, the typical driving posture should be: DPV, diver's trunk and fins are parallel to the cruising direction, DPV is parallel to x-axis; thighs are slightly upward bending; knees and elbows are bending at a certain angle; forearms are parallel to DPV. After quantification, the model in Figure 6 is obtained, in which the angles are assumed to be the positive direction of acute angle with x-axis.
That is, θp1 is the acute angle between the body's central axis and the x-axis; θp1 is the acute angle between the body's central axis and the x-axis; θp2 and θp3 is the acute angle between the upper arm's central axis and the x-axis, and θp3=θp2; θp4 is the acute angle between the thigh's central axis and the x-axis; θp5 is the acute angle between the crus's central axis and the x-axis; For the convenience of calculation, the diver body posture at design speed is given as T 0 posture. According to the Professional Association of Diving Instructors (PADI) DPV driving course, the typical driving posture should be: DPV, diver's trunk and fins are parallel to the cruising direction, DPV is parallel to x-axis; thighs are slightly upward bending; knees and elbows are bending at a certain angle; forearms are parallel to DPV. After quantification, the model in Figure 6 is obtained, in which the angles are assumed to be the positive direction of acute angle with x-axis. In order to make the hydrodynamic modeling of this paper quantitative and more universal, the dimensions and driving posture of four typical DPV as SEADOO 、 U.FLIGHT, SUEX, ROTINOR are analyzed through image data Figure 7, and the following Table 2 is obtained. That is, θ p1 is the acute angle between the body's central axis and the x-axis; θ p1 is the acute angle between the body's central axis and the x-axis; θ p2 and θ p3 is the acute angle between the upper arm's central axis and the x-axis, and θ p3 = θ p2 ; θ p4 is the acute angle between the thigh's central axis and the x-axis; θ p5 is the acute angle between the crus's central axis and the x-axis; In order to make the hydrodynamic modeling of this paper quantitative and more universal, the dimensions and driving posture of four typical DPV as SEADOO, U.FLIGHT, SUEX, ROTINOR are analyzed through image data Figure 7, and the following Table 2 is obtained.
In this paper, a typical diver's body is selected, equipped with back-flying BCD (Buoyancy Control Device) and open breathing apparatus [23]. According to Table 1, the posture angles of DPV with different models and speeds are slightly different. Since the shape and max cruising speed of DPV in this paper are similar to ROTINOR, the T 0 posture angle is selected as Table 3. In order to make the hydrodynamic modeling of this paper quantitative and more universal, the dimensions and driving posture of four typical DPV as SEADOO 、 U.FLIGHT, SUEX, ROTINOR are analyzed through image data Figure 7, and the following Table 2 is obtained.  In this paper, a typical diver's body is selected, equipped with back-flying BCD (Buoyancy Control Device) and open breathing apparatus [23]. According to Table 1, the posture angles of DPV with different models and speeds are slightly different. Since the shape and max cruising speed of DPV in this paper are similar to ROTINOR, the T0 posture angle is selected as Table 3.  In addition, the elbow joint angle θ p3 of all the four DPV is within 12%, and the knee joint angle θ p5 -θ p4 of the other three DPV is within 7.4% except SEADOO. Therefore, in order to further facilitate the calculation, the model is simplified as follows: (1) The angle of knee joint and elbow joint remained unchanged at T 0 ; (2) All joints only rotate around y-axis; (3) The DPV is rigidly fixed to the forearm and parallel to the x-axis when cruising; (4) All solids have the same density as water According to the above conditions, the diver-DPV coupling model are divided into three articulated 6-DOF rigid bodies. As shown in Figure 8, they are named DPV and arms, body and hip, and legs respectively. The hinge points of them are P2 and P4, and the quality attributes are shown in Table 4. Three reference coordinate systems of them named O 1 -X 1 Y 1 Z 1 , O 2 -X 2 Y 2 Z 2 , O 3 -X 3 Y 3 Z 3 are established. The symbols for the force components are defined in Figure 8.

Force Components and Disturbance Analysis
The causes of the external force components are explained as follows: According t the driving principle of DPV, when the propeller is started, the DPV generates axial thrus Ft, which is transmitted to the diver's arms. At the shoulder hinge point P2, FP2Y, and FP2 are applied to the body and hips. At the same time, the DPV and arms are subjected t

Force Components and Disturbance Analysis
The causes of the external force components are explained as follows: According to the driving principle of DPV, when the propeller is started, the DPV generates axial thrust F t , which is transmitted to the diver's arms. At the shoulder hinge point P 2 , F P2Y , and F P2X are applied to the body and hips. At the same time, the DPV and arms are subjected to hydrodynamic forces, which are decomposed into F l,da and F f,da . Body and hips, and legs are also subject to the tension at the hinge point and hydrodynamic forces, and the fluid-solid coupling and solid-solid coupling between multi rigid bodies exist in it [24].
Next, the coupling disturbance components are analyzed: for DPV and arms, in addition to the theoretical value calculated from hydrodynamic coefficient, the disturbance generated in the flow field due to the motion and posture of body and hip, and legs need to add a disturbance term F dis1 . According to the external force analysis of the 6-DOF body in the Figure, the disturbance terms of DPV and arms in the x-axis and z-axis directions of geodetic coordinate system 0-XYZ at the initial T0 time are obtained as where where F f,da = 1/2C t,da ρ w v(t) 2 S da C L,da is the lift coefficient, ρ w is the density of water, F P is the hinge point force between the 6-DOF bodies, F b is the buoyancy of the 6-DOF bodies, G is the gravity of the 6-DOF bodies, F t is the propeller trust, F dis is the disturbing force of the other two 6-DOF bodies. Similarly, for body and hip for legs, we get m legdaZ,legs (t) = G legs + F P4y − F b,legs − F dis3,y To sum up, when the DPV curise at a uniform speed, the acceleration a is zero. Combined with the above formula, along the x-axis of the geodetic coordinate system, there are where F f,da + F f,bh + F f,legs = F f0 . That is, the resistance in the direction of x-axis is restrained. F dis1,x + F dis2,x + F dis3,x = F dis,x . That is, the posture disturbance term in the x-axis direction.
In the fourth section, by solving and discussing the size and direction of F dis , the disturbance factors caused by diver motion are determined.

Computational Region and Meshing
Firstly, the overlapping mesh is generated. In the common CFD calculation of underwater vehicle, the background region is usually set as the water tank involved in. While in order to save computational resources, the small area around the model and the area of its wake flow always set finer grid, which is called the encryption region, and the other flow area uses thicker grid.
In order to fully develop the incoming flow and avoid the interference from the boundary of the external field, the size of the background region (water tank) is adjusted to be 4-times the total length of the diver-DPV coupled model, 2.5-times of the full width, and 5-times of the full height. Furthermore, the encryption grid area is set to increase the calculation accuracy. After integer, the size is as follows in Table 5. The side length of the grid in the water tank region is 10cm, and the side length of the grid in the encryption area is 1cm (0.4% of the model length). The grid number of water tank is 6704321, and that of the encryption is 1052876. Finally, the calculated areas are shown in Figure 9.
J. Mar. Sci. Eng. 2021, 9, x FOR PEER REVIEW The side length of the grid in the water tank region is 10cm, and the side l the grid in the encryption area is 1cm (0.4% of the model length). The grid nu water tank is 6704321, and that of the encryption is 1052876. Finally, the calculat are shown in Figure 9. In order to simulate the disturbance of diver joint motion on the hydrodyna formance simulation calculation, the following provisions are made on the basis o division: (1) The rigid body of DPV and arms is bound on the water tank region, that is, th moves together with the motion of DPV, which can ensure that the DPV wil out of the water tank during cruising; (2) The dynamic region is set based on Body and Hip and Legs. In order to imp calculation accuracy, its contour envelops the whole diver-DPV model, shown in Figure 10. (3) In particular, the overlapping grid relationship is established between the su the three overlapping mesh continuum to ensure that the fluid converges wh ing through the interface. (4) The encrypted region is set large enough to ensure that the overlapped grid i included in when in motion, so as not to diverge due to different scales after ing with the external flow field Finally, the side length of the overlapped grid region is also 1cm (the sam encryption region), and the overlapped grid number of body and hip, and legs ar and 462466. The number of grids in all regions is with a total of 7659620. In order to simulate the disturbance of diver joint motion on the hydrodynamic performance simulation calculation, the following provisions are made on the basis of region division: (1) The rigid body of DPV and arms is bound on the water tank region, that is, the region moves together with the motion of DPV, which can ensure that the DPV will not run out of the water tank during cruising; (2) The dynamic region is set based on Body and Hip and Legs. In order to improve the calculation accuracy, its contour envelops the whole diver-DPV model, which is shown in Figure 10. (3) In particular, the overlapping grid relationship is established between the surfaces of the three overlapping mesh continuum to ensure that the fluid converges when passing through the interface. (4) The encrypted region is set large enough to ensure that the overlapped grid is always included in when in motion, so as not to diverge due to different scales after contacting with the external flow field.
Finally, the side length of the overlapped grid region is also 1cm (the same as the encryption region), and the overlapped grid number of body and hip, and legs are 492833 and 462466. The number of grids in all regions is with a total of 7659620. After all the regions are established, the DFBI (Dynamic Fluid Body Interaction) motion function of Star-CCM+ software is used to make the 6-DOF body produce constrained motion. The DFBI body model can simulate the passive motion of rigid body under the action of fluid. That is, the DOF of the rigid body can be determined according to the needs, and the specific DOF can be constrained/released. The mass, moment of inertia, and centroid coordinates of the three rigid bodies of diver-DPV model are described in Table 2, while the specific degrees of freedom constraints of each DFBI body are described in Section 3.2. The joint and shoulder coordinate are needed to be supplement: The upper, lower and the inflow surface of the water tank region are set as the velocity-inlet, the back surface is the pressure-outlet, and the left-and right-side surface are set as symmetrical-planes. An overlapping grid boundary is established between the outer surface of the leg region, the body region, and the water tank region, while the model surfaces in the three overlapping regions are set as the wall-surface. The k-ω SST turbulence model is used to close the N-S equation. The full y+ wall function method suitable for high Reynolds number is used for wall treatment. First boundary layer mesh is obtained by y+=30, and the case y+ distribution is shown in Figure 11. Finally, 10 layers of prismatic mesh are adopted, and the thickness of the first prismatic mesh is 0.0073 mm according to y+. After all the regions are established, the DFBI (Dynamic Fluid Body Interaction) motion function of Star-CCM+ software is used to make the 6-DOF body produce constrained motion. The DFBI body model can simulate the passive motion of rigid body under the action of fluid. That is, the DOF of the rigid body can be determined according to the needs, and the specific DOF can be constrained/released. The mass, moment of inertia, and centroid coordinates of the three rigid bodies of diver-DPV model are described in Table 2, while the specific degrees of freedom constraints of each DFBI body are described in Section 3.2. The joint and shoulder coordinate are needed to be supplement: The boundary conditions of regions are set as follows: The upper, lower and the inflow surface of the water tank region are set as the velocity-inlet, the back surface is the pressure-outlet, and the left-and right-side surface are set as symmetrical-planes. An overlapping grid boundary is established between the outer surface of the leg region, the body region, and the water tank region, while the model surfaces in the three overlapping regions are set as the wall-surface. The k-ω SST turbulence model is used to close the N-S equation. The full y+ wall function method suitable for high Reynolds number is used for wall treatment. First boundary layer mesh is obtained by y+=30, and the case y+ distribution is shown in Figure 11. Finally, 10 layers of prismatic mesh are adopted, and the thickness of the first prismatic mesh is 0.0073 mm according to y+.
In the aspect of solver, the time step is selected first. When the time step is too small, the calculation time is too long, and when the time step is too large, the calculation accuracy is not enough. The resistance time step is generally 0.005 L/V~0.01 L/V, where V is the velocity, and the corresponding Courant number of this model is 0.2-0.4. Therefore, 0.005 s is taken as the time step, which satisfies the CFL condition. Set the maximum number of internal iterations to 5. The SIMPLE (Semi Implicit Method for Pressure Linked Equation) algorithm is used to solve the Navier-Stokes equations. The viscous term in the discrete equation adopts the second-order central difference scheme, and the convective term adopts the second-order upwind scheme.
In terms of initial conditions, assuming that the experimental condition is infinite water depth, the fresh water temperature is 24°C, g = 9.81 m/s 2 , and the kinematic viscosity coefficient v = 0.9167 × 10 6 m 2 /s. The necessary DOF of the rigid bodies are constrained, and the x-axis cruising is simulated by uniform inflow, that is, the absolute motion is replaced by relative motion to save calculation time. In the aspect of solver, the time step is selected first. When the time step is too small, the calculation time is too long, and when the time step is too large, the calculation accuracy is not enough. The resistance time step is generally 0.005L/V~0.01L/V, where V is the velocity, and the corresponding Courant number of this model is 0.2-0.4. Therefore, 0.005s is taken as the time step, which satisfies the CFL condition. Set the maximum number of internal iterations to 5. The SIMPLE (Semi Implicit Method for Pressure Linked Equation) algorithm is used to solve the Navier-Stokes equations. The viscous term in the discrete equation adopts the second-order central difference scheme, and the convective term adopts the second-order upwind scheme.
In terms of initial conditions, assuming that the experimental condition is infinite water depth, the fresh water temperature is 24 ℃, g = 9.81 m/s 2 , and the kinematic viscosity coefficient v = 0.9167*10 6 m 2 /s. The necessary DOF of the rigid bodies are constrained, and the x-axis cruising is simulated by uniform inflow, that is, the absolute motion is replaced by relative motion to save calculation time.

Results
This section shows the CFD results the straight-line cruising motion of the diver-DPV coupling model in the water. That is, the motion in the x-axis direction is simulated through the uniform incoming flow, and the translation movement in the z-axis direction is constrained. Firstly, the diver and DPV are regarded as a rigid body by restraining the human motion; then the human motion is released as a multi-rigid-articulated model; finally, the diver-DPV disturbance components are analyzed by comparison.

Case of Restricting Human Motion
This section only simulates constrained straight-line cruising. All DOF of the three 6-DOF models are fixed and restrained.
Firstly, mesh independence and convergence are analyzed by taking the incoming velocity with the DPV's maximum speed V5=2.5 m/s as an example. The sub relaxation factor of velocity and pressure is 0.5. After 10000 time steps of simulation, the convergence result of resistance Ff is shown in Figure 12. The velocity nephogram is shown in Figure  13.

Results
This section shows the CFD results the straight-line cruising motion of the diver-DPV coupling model in the water. That is, the motion in the x-axis direction is simulated through the uniform incoming flow, and the translation movement in the z-axis direction is constrained. Firstly, the diver and DPV are regarded as a rigid body by restraining the human motion; then the human motion is released as a multi-rigid-articulated model; finally, the diver-DPV disturbance components are analyzed by comparison.

Case of Restricting Human Motion
This section only simulates constrained straight-line cruising. All DOF of the three 6-DOF models are fixed and restrained.
Firstly, mesh independence and convergence are analyzed by taking the incoming velocity with the DPV's maximum speed V 5 = 2.5 m/s as an example. The sub relaxation factor of velocity and pressure is 0.5. After 10000 time steps of simulation, the convergence result of resistance F f is shown in Figure 12. The velocity nephogram is shown in Figure 13.   Then the straight-line thrust Ft required by V5=2.5 m/s is obtained. Because of the resistance tends to converge after 1s, the resistance value in the period of 1-3 s is taken as the arithmetic average. That is, Ff=Ft=300.4N，z-axis lift Fl =-14.2N.
Next, the resistance prediction of constrained straight-line cruising under V5=2.5 m/s is selected for grid independent. According to ITTC (International Towing Tank Conference) [25], the grid side length is set as 1/2000-5/2000 of the model length L, which is 1 cm in this paper. There set three different length as 1 cm,√2 cm, and 2 cm to generate coarse mesh. The total number of meshes is 7659620, 2708494, and 957453. When the other parameters and settings are the same, the change trend of resistance value with the foundation size is observed as Figure 14. Then the straight-line thrust F t required by V 5 = 2.5 m/s is obtained. Because of the resistance tends to converge after 1s, the resistance value in the period of 1-3 s is taken as the arithmetic average. That is, Next, the resistance prediction of constrained straight-line cruising under V 5 = 2.5 m/s is selected for grid independent. According to ITTC (International Towing Tank Conference) [25], the grid side length is set as 1/2000-5/2000 of the model length L, which is 1 cm in this paper. There set three different length as 1 cm, √ 2 cm, and 2 cm to generate coarse mesh. The total number of meshes is 7659620, 2708494, and 957453. When the other parameters and settings are the same, the change trend of resistance value with the foundation size is observed as Figure 14. Finally, the mean resistance Ff,mesh1=300.4N, Ff,mesh2=321.2N, Ff,mesh3=378.1N. So there are: Therefore, with the gradual refinement of the grid, the numerical simulation value of Finally, the mean resistance F f,mesh1 = 300.4 N, F f,mesh2 = 321.2 N, F f,mesh3 = 378.1 N. So there are: Therefore, with the gradual refinement of the grid, the numerical simulation value of resistance gradually decreases, and the calculation accuracy gradually converges as 6.9%, and the number of grids increases sharply as 180%. Considering the calculation accuracy and cost, it is reasonable to take the grid side length of 1 cm.
Then, the numerical simulation of constrained straight-line cruising under the other four speed conditions is carried out. The four inflow velocities are set as V 1 = 0.5 m/s, V 2 = 1.0 m/s, V 3 = 1.5 m/s, and V 4 = 2.0 m/s. So the linear relationship between V n 2 and F f is obtained. The resistance coefficient C t is introduced here. Because the rigid body under water straight-line cruising has F f = 1/2C t ρ w v 2 S, the ideal value C t0 can be obtained from the slope of the fitting curve. After numerical simulation, the relationship between resistance F f and V 2 is shown in Figure 15. Finally, the resistance coefficient of diver-DPV model is obtained as Ct0= 27.2984/10 −3 .

Case of Releasing Human Motion
This section only carries out the simulation of straight-line cruising with released human motion. That is, releasing the x-and z-axis translation DOF and the rotation DOF around the y-axis of two 6-DOF bodies of body and hips, and legs to simulate the articulated motion of diver body in water.

Maximum Cruising Speed
Firstly, the numerical simulation of the maximum cruising speed V5=2.5 m/s is carried out. The velocity clouds at the moment of 2 s is shown in Figure 16  Finally, the resistance coefficient of diver-DPV model is obtained as C t0 = 27.2984/10 −3 .

Case of Releasing Human Motion
This section only carries out the simulation of straight-line cruising with released human motion. That is, releasing the xand z-axis translation DOF and the rotation DOF around the y-axis of two 6-DOF bodies of body and hips, and legs to simulate the articulated motion of diver body in water.

Maximum Cruising Speed
Firstly, the numerical simulation of the maximum cruising speed V 5 = 2.5 m/s is carried out. The velocity clouds at the moment of 2 s is shown in Figure 16. And the variation curve of resistance F f and lift force F l are shown in Figures 17 and 18. The posture angle changes of 6-DOF model legs' θ 1 , and body and hip's θ 4 are obtained as shown in Figures 19 and 20.

Maximum Cruising Speed
Firstly, the numerical simulation of the maximum cruising speed V5=2.5 m/s is carried out. The velocity clouds at the moment of 2 s is shown in Figure 16. And the variation curve of resistance Ff and lift force Fl are shown in Figures 17 and 18. The posture angle changes of 6-DOF model legs' θ1, and body and hip's θ4 are obtained as shown in Figures  19 and 20.   The analysis shows that the angle θP1 of Body and Hip gradually converges to the stable range, and θP4 of legs first fluctuates in a certain range, and finally converges to the stable range. That is, the diver body will swing up and down during cruising. At the same time, the resistance also fluctuates in a certain range. Because the resistance converge after 1s and fluctuates in a small range,and the human body posture reached the inflection point angle at 2 s, the resistance value in the period of 1-2 s is taken as the arithmetic average. That is, Ff,V5 = 272.72N, and the lift Fl,V5 finally converges near 0.

Other Cruising Speeds
Next, the numerical simulation of other speeds is carried out. The x-axis and z-axis translation and rotation around y-axis of body and hip, and legs are released, and the zaxis translation of DPV and arms is restrained. Set five inflow velocities as V1=0.5 m/s, V2=1.0 m/s, V3=1.5 m/s, V4=2.0 m/s, and V5=2.5 m/s, then the posture angle of diver body and its influence on the speed of straight-line cruising are obtained.
In fluid mechanics, a dimensionless variable characterizing the relative magnitude of inertial force and gravity of fluid is called Froude number Fr, which represents the ratio of the inertial force to the magnitude of gravity. That is, Fr=v 2 /gL, where v is the straightline cruising speed of rigid body, g is the acceleration of gravity, L is the characteristic length of the object.
In this paper, the resistance of the other four cases of different speeds are calculated according to Fr number from low to high. After CFD convergence, the mean resistance, mean lift and posture are shown in Table 6, and the velocity scalar nephogram is shown The analysis shows that the angle θ P1 of Body and Hip gradually converges to the stable range, and θ P4 of legs first fluctuates in a certain range, and finally converges to the stable range. That is, the diver body will swing up and down during cruising. At the same time, the resistance also fluctuates in a certain range. Because the resistance converge after 1 s and fluctuates in a small range,and the human body posture reached the inflection point angle at 2 s, the resistance value in the period of 1-2 s is taken as the arithmetic average. That is, F f,V5 = 272.72 N, and the lift F l,V5 finally converges near 0.

Other Cruising Speeds
Next, the numerical simulation of other speeds is carried out. The x-axis and z-axis translation and rotation around y-axis of body and hip, and legs are released, and the z-axis translation of DPV and arms is restrained. Set five inflow velocities as V 1 = 0.5 m/s, V 2 = 1.0 m/s, V 3 = 1.5 m/s, V 4 = 2.0 m/s, and V 5 = 2.5 m/s, then the posture angle of diver body and its influence on the speed of straight-line cruising are obtained.
In fluid mechanics, a dimensionless variable characterizing the relative magnitude of inertial force and gravity of fluid is called Froude number Fr, which represents the ratio of the inertial force to the magnitude of gravity. That is, Fr = v 2 /gL, where v is the straight-line cruising speed of rigid body, g is the acceleration of gravity, L is the characteristic length of the object.
In this paper, the resistance of the other four cases of different speeds are calculated according to Fr number from low to high. After CFD convergence, the mean resistance, mean lift and posture are shown in Table 6, and the velocity scalar nephogram is shown in Figure 21. The instantaneous value of posture angle is obtained at 2 s moment, and the resistance is the average value of 1~2 s.  The analysis shows that with the increase of speed V, the Body and Hip posture angle θP1 decreases with the increase of speed, the Legs posture angle θP4 is also the same, and the resistance Ff gradually decreases. The fitting curve of the relationship between Ff and V 2 is shown in Figure 22, and the fitting equation is in Table 7. The analysis shows that with the increase of speed V, the Body and Hip posture angle θ P1 decreases with the increase of speed, the Legs posture angle θ P4 is also the same, and the resistance F f gradually decreases. The fitting curve of the relationship between F f and V 2 is shown in Figure 22, and the fitting equation is in Table 7.
It can be seen that the resistance velocity relationship of traditional submersible is not consistent with F f = 1/2C t ρ w v 2 S, This is due to the change of human body posture angle under the action of flow field, which leads to the real-time change of resistance coefficient C t , and finally tends to the minimum value. The specific process is as follows: According to the Figure 20, the absolute value of lift |F l | under various working conditions has experienced a state from large to small, and finally tends to 0N. Therefore, in straight-line navigation, the lift force brought by the real-time posture angle θ p (t) will drive the human body to approach the posture with the minimum resistance F l (t). When the final lift tends to 0, the resistance is the minimum, and the posture angle reaches the target value.
When the resistance F f converges, there is still a difference between the resistance of the constrained diver-DPV model, which is defined as the posture disturbance term F dis . When F dis > 0, it is called favorable disturbance, and when F dis < 0, it is called adverse disturbance. Its relationship with other resistance components and speed is discussed in the next section. The analysis shows that with the increase of speed V, the Body and Hip posture angle θP1 decreases with the increase of speed, the Legs posture angle θP4 is also the same, and the resistance Ff gradually decreases. The fitting curve of the relationship between Ff and V 2 is shown in Figure 22, and the fitting equation is in Table 7.

Discussion
This section analyzes the relationship between the X-axis resistance disturbance term F dis , speed and posture. Firstly, the total resistance components of the diver-DPV cruising in the water are analyzed. The results show that the total resistance of the diver-DPV coupled model in the water is F f = F f0 + F dis , where F f0 is the resistance value of restricting the direct navigation movement. Ignoring the influence of near surface and bottom effect, the component of the resistance of F f0 is usually divided into the following points: friction resistance R f , residual resistance R rs , Therefore, the total resistance can be expressed as According to Hughes, the ratio of viscous resistance coefficient to frictional resistance coefficient is constant. By introducing the shape factor to deal with the three-dimensional flow of the diver-DPV model, the resistance conversion between the model and the real vehicle is realized. The resistance coefficient relationship is where C tm is the resistance coefficient of the model, (1 + k) is the shape coefficient, C fm is the friction resistance coefficient, C fs is the actual friction resistance coefficient, C f is the roughness conversion subsidy coefficient.
According to the formula of ITTC (International Towing Tank Conference) [25], the friction and resistance coefficient of a real vehicle can be calculated by Reynolds number (Re) where v is the water viscosity coefficient affected by temperature.
The corresponding subsidy coefficient of roughness conversion can be expressed as K s is the height of rough surface. This CFD model is a full-scale model. The fresh water temperature is 25 • C and the kinematic viscosity coefficient of water is v = 0.9167 × 10 6 m 2 /s. Without considering the cruising wear, this paper takes K s = 150 × 10 −6 . In conclusion, according to the main scale of diver-DPV model, the roughness conversion subsidy coefficient can be obtained, and then the resistance coefficient conversion between the CFD model and the actual DPV can be realized.
In order to effectively analyze the change of resistance components, the model resistance is dimensionless. Assuming that the wet surface area of diver-DPV model is constant, there is still friction resistance f = 1/2C f ρ w v 2 S. The empirical formula of friction resistance coefficient is The total resistance coefficient is still calculated as F f = 1/2C t ρ w v 2 S, while That is: C t = C t0 + C dis = C f + C rs + C f + C dis (22) According to that the wet surface area S w = 3.521381 m 2 , we get C t0 = 27.2984/10 −3 . The results are shown in Table 8. According to Table 8, the resistance coefficient C t of the diver-DPV coupled model is plotted as Figure 23 according to the variation curve.
It can be seen from Figure 23 that the total resistance coefficient of model tends to decrease with the increase of Froude number F r . The friction resistance R fs in those speed cases is the main part of the total resistance, and the correlation between the value and the speed change trend gradually decreases. When F r > 0.06, C dis < 0, which means it is favorable disturbance; when F r < 0.06, C dis > 0, which means it is adverse disturbance.  Table 8, the resistance coefficient Ct of the diver-DPV coupled model is plotted as Figure 23 according to the variation curve.  According to various resistance coefficients, the actual resistance F f , resistance of restraint movement F f0 , disturbance term F dis and speed V are shown in Table 9. The fitting curves of F dis and V n is as follows in Figure 24: It can be seen from Figure 23 that the total resistance coefficient of model tends to decrease with the increase of Froude number Fr. The friction resistance Rfs in those speed cases is the main part of the total resistance, and the correlation between the value and the speed change trend gradually decreases. When Fr > 0.06, Cdis < 0, which means it is favorable disturbance; when Fr < 0.06, Cdis > 0, which means it is adverse disturbance.
According to various resistance coefficients, the actual resistance Ff, resistance of restraint movement Ff0, disturbance term Fdis and speed V are shown in Table 9. The fitting curves of Fdis and Vn is as follows in Figure 24：   Finally, the relationship between the total resistance Ff and Vn is obtained as Finally, the relationship between the total resistance F f and V n is obtained as The diver-DPV disturbance in the x-axis is m . V n (t) = 1 2 C t0 ρ w V n (t) 2 S + A + B * V n (t) + C * V n (t) 2 + D * V n (t) 3 The values of the constants are shown in the Table 10. From the conclusion, it can be concluded that in the design of DPV, if the maximum cruising speed index is too low, the diver's attitude will always be in the range of adverse disturbance, and the navigation efficiency will be reduced. However, if the maximum cruising speed is too high, the disturbance coefficient will not change, and the cruising efficiency will not be further improved. It can be understood that the high-speed flow field has already impacted the diver's body to the horizontal attitude, and will not change further. Through this CFD model, we can get the most suitable speed range index in the DPV design, which is of great significance to the design of various shape of DPV in the future.

Conclusions
In summary, by analyzing the existing DPV (diver propulsion vehicle) equipment, a set of diver-DPV multi-body coupling model considering rigid body dynamics and fluid disturbance is established in this paper. The split and motion of human joints is realized by overlapping grid of Star-CCM+ software and DFBI 6-DOF model motion method. The nonlinear disturbance caused by diver posture transformation under 10 different cruising cases and its influence on the rapidity and diver posture characteristics of DPV are concluded as follows: When a diver driving DPV and cruising straight line, the angle of diver body's posture is related to the speed of the vehicle. The movement of this posture change produces a disturbing potential to the surrounding flow field, and produces a resistance disturbance term to the diver-DPV coupled model, which makes the overall cruising resistance different from the resistance at the same speed under constraint human body conditions. The relationship between the sign and volume of the disturbance term F dis and the speed is as follows: F dis = A + B*V n + C*V n 2 + D*V n 3 , when Fr > 0.06, F dis < 0, which means it is favorable disturbance and can reduce cruising resistance. When Fr < 0.06, F dis > 0, which means it is adverse disturbance. In addition, within the design maximum speed of DPV, the friction resistance F f is always the main part.
Through this multi-body hydrodynamic model, we can get the most suitable speed range index in the DPV design. It is of great significance for the future design of DPV with various shapes to improve the underwater cruising efficiency. Acknowledgments: First of all, I would like to appreciate for Harbin DepTech marine technology Co., Ltd. and other corporations for the DPV product data. Second, thanks to every scuba diving instructors from PADI for their detailed introduction on the professional DPV driving posture. Finally, I would like to thank all the editors and reviewers for their careful guidance and help. This is the first time that I have written an international top-level journal paper as the first author. Thanks a lot for your tolerance of minor errors in the process of paper revision. Wish there will be more cooperation in the future.

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

Re min
Reynolds number of minimum rigid body in multi joint model Acceleration of x-axis and y-axis cruising P n Hinge points of multi rigid body joint model, n = 1, 2, 3, 4, 5 θ Pn Angle between diver joint and x-axis at hinge point n C t Total resistance coefficient C dis Disturbance resistance coefficient C f Friction resistance coefficient C rs Residual resistance coefficient C f Rough subsidy coefficient C dis Disturbance coefficient V n Underwater cruising speed of DPV, n = 1, 2, 3, 4, 5 S w Wet surface area of the model I x , I y , I z Moment of inertia for X, Y, Z axis F t Propeller thrust F f Cruising resistance of the diver-DPV model F f,Vn Cruising resistance at different speeds F f,meshn Resistance at maximum speed V 5 with different computational grid, n = 1, 2, 3 F b Buoyancy F l Lift force F l,da Z-axis lift force acting on DPV and arms F f,da Resistance in x-axis direction acting on DPV and arms F b,da Buoyancy of DPV and arms G da Gravity of DPV and arms F P2X , F P2Y Tensile force at hinge point P 2 F l,bh Lift force acting on body and hip in z-axis direction F f,bh Resistance in x-axis direction acting on body and hip F b,bh Buoyancy of body and hip G bh Gravity of body and hip F P4Y ,F P4X Tensile force at hinge point P 4 F l,leg Z-axis lift acting on legs F f,legs Resistance in x-axis direction acting on legs F b,legs Buoyancy of legs G legs Gravity of legs F dis Disturbance force F disn,x Disturbance force on different 6-DOF body in x-axis, n = 1, 2, 3