Modeling Analysis and Simulation of Viscous Hydrodynamic Model of Single-DOF Manipulator

Hydrodynamic modeling is the basis of the precise control research of underwater manipulators. Viscous hydrodynamics, an important part of the hydrodynamic model, directly affects the accuracy of the dynamic model and the control model of the manipulator. Considering the limited research on viscous hydrodynamics of underwater manipulators and the difficulty in measuring viscous hydrodynamic coefficients, the viscous hydrodynamic model in the form of Taylor expansion is analyzed and established. Through carrying out simulation calculations, curve fitting and regression analysis, positional derivatives, rotational derivatives, and coupling derivatives in the viscous hydrodynamic model, are determined. This model provides a crucial theoretical foundation and reference data for subsequent related research.


Introduction
With the rapid development of the marine industry, the application of robots mounted on underwater vehicles has become increasingly widespread. Underwater manipulators are mainly employed in the acquisition and exploration, etc. Their underwater performance can be determined by their speed and accuracy. Underwater manipulators are usually subjected to large forces and moments during their operation. These forces and moments mainly include gravity, inertial hydrodynamics, and viscous hydrodynamics. There are more solutions to the determination of gravity and inertial hydrodynamics. Therefore, the viscous hydrodynamic model has become a research focus that needs to be broken through because of its large number of coefficients and the difficulty of measurement.
References [1][2][3] verified the accuracy of fluid simulation calculations using software, such as FLUENT, by comparing the simulated calculation values of hydrodynamics with the measured values of engineering methods. References [4][5][6] studied the drag and additional mass force of the manipulator underwater and obtained the inertial hydrodynamic model by the CFD (Computational Fluid Dynamics) simulation, but the viscous hydrodynamics were not considered as key targets in the study. References [7][8][9][10] numerically calculated the viscous hydrodynamic coefficients of submarine and ship models with complex shapes under the specified navigation conditions and further improved the handling performance during navigation. Reference [11] determined the viscous kinetic coefficients of the manipulator through aerodynamic experiments and carried out a sea trial.
At present, research on viscous hydrodynamics mostly focus on specific ship models and underwater vehicles, and less on manipulators loaded on them. Therefore, a numerical calculation method of the viscous hydrodynamic coefficient based on single-DOF (Degree of Freedom) manipulators is put forward. The method is based on ANSYS Fluent. UDF (user-defined function) and dynamic mesh which are used to simulate the rotational motion of the manipulator. The inlet and outlet conditions, as well as boundary conditions of the fluid domain, are changed. Also, the viscous hydrodynamic model could be obtained by curve fitting and regression analysis. and outlet conditions, as well as boundary conditions of the fluid domain, are changed. Also, the viscous hydrodynamic model could be obtained by curve fitting and regression analysis.

Mathematical Modeling of Viscous Hydrodynamics
Since its surrounding flow field is changed during the operation of the underwater manipulator; the arm is subjected to the reaction of the water body caused by force generated by the water body. The hydrodynamics of the moving manipulator in water include two parts: Inertial hydrodynamics caused by acceleration and viscous hydrodynamics caused by friction. Inertial hydrodynamics is generally determined by empirical formulas, while coefficients of viscous hydrodynamics are complex and difficult to measure.
This paper is intended to use a single-DOF manipulator model mounted on a fixed base. The arm is a homogeneous lightweight rod with a circular cross-section, 50 mm in diameter (d), 500 mm in arm length (l), at length to diameter ratio of 10, belonging to the slender pole. The establishment of its coordinate system is shown in Figure 1. When the manipulator performs the constant motion, the influence of acceleration and angular acceleration on the motion variable could be ignored based on the "slow-motion" assumption. So the mere considerations for viscous hydrodynamics are the effects of the velocity and angular velocity in the motion variable. It can be represented as: [ ] T U = u v w p q r , where u stands for the velocity in the OX direction, v for the velocity in the OY direction, w for the velocity in the OZ direction, p for the angular velocity of the rotation around the OX axis, q for the angular velocity of the rotation around the OY axis, and r for the angular velocity of the rotation around the OZ axis. The viscous hydrodynamics can be expressed as a multivariate function  When the manipulator performs the constant motion, the influence of acceleration and angular acceleration on the motion variable could be ignored based on the "slow-motion" assumption. So the mere considerations for viscous hydrodynamics are the effects of the velocity and angular velocity in the motion variable. It can be represented as: U = [u v w p q r] T , where u stands for the velocity in the OX direction, v for the velocity in the OY direction, w for the velocity in the OZ direction, p for the angular velocity of the rotation around the OX axis, q for the angular velocity of the rotation around the OY axis, and r for the angular velocity of the rotation around the OZ axis.
The viscous hydrodynamics can be expressed as a multivariate function F = f(u, v, w, p, q, r), and the six-dimensional component of viscous hydrodynamics F = [X Y Z K M N] T could also be displayed in the above functional form, the direction of which is shown in Figure 2. and outlet conditions, as well as boundary conditions of the fluid domain, are changed. Also, the viscous hydrodynamic model could be obtained by curve fitting and regression analysis.

Mathematical Modeling of Viscous Hydrodynamics
Since its surrounding flow field is changed during the operation of the underwater manipulator; the arm is subjected to the reaction of the water body caused by force generated by the water body. The hydrodynamics of the moving manipulator in water include two parts: Inertial hydrodynamics caused by acceleration and viscous hydrodynamics caused by friction. Inertial hydrodynamics is generally determined by empirical formulas, while coefficients of viscous hydrodynamics are complex and difficult to measure.
This paper is intended to use a single-DOF manipulator model mounted on a fixed base. The arm is a homogeneous lightweight rod with a circular cross-section, 50 mm in diameter (d), 500 mm in arm length (l), at length to diameter ratio of 10, belonging to the slender pole. The establishment of its coordinate system is shown in Figure 1.
where U k = u 0 v 0 w 0 p 0 q 0 r 0 T , U k is the initial constant. .
Considering the left-right and front-back symmetry and the motion limit of the manipulator, several coefficients in the second-order expansions are zero, and the remaining others are non-negligible hydrodynamic coefficients, as in Equation (2): Viscous hydrodynamics of the manipulator could be calculated in FLUENT. The regression analysis of the calculated viscous hydrodynamics and corresponding velocities could be performed to obtain unknown coefficients in Equation (2) in MATLAB. The first derivative coefficients associated only with the linear velocity (u, v, w) are the position derivatives, and the first derivative coefficients related to the angular velocity (p, q, r) are the rotational derivatives. The coefficients caused by joint changes of two or more parameters are the coupling derivatives.

Control Equation
To analyze hydrodynamics of underwater manipulators, it is generally assumed that the fluid is isothermal and incompressible, also as a constant flow that magnitude and direction do not change with time. The continuity equation (Equation (3)) and the Navier-Stokes equation (Equation (4)) serve as the two basic equations necessary to solve the flow problem of viscous fluid. These equations are generally described in the form of partial differentials: The form of expression of the time-averaged continuity equation does not change. Instead, the tensor of the Reynolds stress is added to the formula after the N-S equation time-averaged, which leads to the closure problem of the original equation. The Reynolds stress is about 10 −2 Pa, which is indicative of the turbulent flow and cannot be directly ignored. Therefore, it is necessary to introduce a proper turbulence model to make a modeling description of the increased Reynolds stress.
For the incompressible isothermal turbulent water environment, the basic equations of turbulence include the DNS equations (direct numerical simulation), LES equations (large eddy simulation), and RANS equations (Renault time average). The two equations (DNS equations and LES equations) are of limited use due to the requirements of a large number of computational grids. At present, most engineering calculations adopt the RANS equation to solve closed equations formed by introducing the turbulence model and, thus, obtain the time-average value of the turbulent elements. To solve the viscous hydrodynamics under steady-state, the appropriate turbulence model is the key to the numerical simulations in this paper.

Calculation Domain Establishment
Generally speaking, the larger the size of the calculation domain is, the closer it is to the real working condition. The downside is that it will increase the amount of calculation and prolong the calculation period. If the calculation domain is too small, the boundary conditions and calculation results are difficult to match the real working conditions. Therefore, it is very important to choose the size of the calculation domain reasonably.
Based on previous experience and multiple numerical practices [6,12,13], this paper establishes the domain as the computational domain, shown in Figure 3b. The specific sizes are as follows: For the incompressible isothermal turbulent water environment, the basic equations of turbulence include the DNS equations (direct numerical simulation), LES equations (large eddy simulation), and RANS equations (Renault time average). The two equations (DNS equations and LES equations) are of limited use due to the requirements of a large number of computational grids. At present, most engineering calculations adopt the RANS equation to solve closed equations formed by introducing the turbulence model and, thus, obtain the time-average value of the turbulent elements. To solve the viscous hydrodynamics under steady-state, the appropriate turbulence model is the key to the numerical simulations in this paper.

Calculation Domain Establishment
Generally speaking, the larger the size of the calculation domain is, the closer it is to the real working condition. The downside is that it will increase the amount of calculation and prolong the calculation period. If the calculation domain is too small, the boundary conditions and calculation results are difficult to match the real working conditions. Therefore, it is very important to choose the size of the calculation domain reasonably.
Based on previous experience and multiple numerical practices [6,12,13], this paper establishes the domain as the computational domain, shown in Figure 3b. The specific sizes are as follows:  The motion parameters of the simulation are shown in Figure 3, α stands for the rotation angle around the OZ axis, β for the angle of the arm to the flow direction, V for the inlet flow velocity of the domain, and the linear velocity components of the manipulator are as shown in Equation (5):

Meshing
The calculation of CFD requires well-distributed grids. Usually, grids are classified into structured and unstructured grids. The unstructured grid means that internal points in the mesh area do not have the same adjacent unit, with no regular topology, not arranged in layers, and the distribution of mesh nodes is arbitrary. Therefore, it is more flexible than structured grid. The motion parameters of the simulation are shown in Figure 3, α stands for the rotation angle around the OZ axis, β for the angle of the arm to the flow direction, V for the inlet flow velocity of the domain, and the linear velocity components of the manipulator are as shown in Equation (5):

Meshing
The calculation of CFD requires well-distributed grids. Usually, grids are classified into structured and unstructured grids. The unstructured grid means that internal points in the mesh area do not have the same adjacent unit, with no regular topology, not arranged in layers, and the distribution of mesh nodes is arbitrary. Therefore, it is more flexible than structured grid. Unstructured grids can be optimized by using certain criteria in the process of their generation. Ultimately, they can be displayed as high-quality meshes, which are suitable for complex geometry, easy to control grid size, and node density. Moreover, the adoption of random data structures is conducive to mesh adaptation. It is difficult to model the structured mesh due to the shape of the manipulator, instead, the unstructured mesh is easy to combine with the dynamic mesh technology. Thus, the unstructured mesh is adopted in the research.
We compare basic mesh to dense mesh for the sake of grid independence verification. Basic mesh in the vicinity of the manipulator is shown in Figure 4a and dense mesh is shown in Figure 4b. The number of elements and nodes are shown in Table 1. Unstructured grids can be optimized by using certain criteria in the process of their generation. Ultimately, they can be displayed as high-quality meshes, which are suitable for complex geometry, easy to control grid size, and node density. Moreover, the adoption of random data structures is conducive to mesh adaptation. It is difficult to model the structured mesh due to the shape of the manipulator, instead, the unstructured mesh is easy to combine with the dynamic mesh technology. Thus, the unstructured mesh is adopted in the research. We compare basic mesh to dense mesh for the sake of grid independence verification. Basic mesh in the vicinity of the manipulator is shown in Figure 4a and dense mesh is shown in Figure 4b. The number of elements and nodes are shown in Table 1.   The rotation process of the manipulator is realized by UDF and dynamic mesh technology in Fluent. Dynamic meshing is performed by using the spring smoothing model, which approximates connecting lines between the grid nodes as springs, and the position of nodes after smoothing is obtained by calculating force equilibrium equations between them. In the calculation process, meshes with a large aberration rate or huge change in size are grouped together to re-divide the partial meshes.

The Solution of Position Derivatives, Rotation Derivatives, and Coupling Derivatives
Solving the unknown coefficients in Equation (2) is the key to the research task, where Yv, Nv, Zw, and Mw are positional derivatives, Yr and Nr are rotational derivatives, and the remaining unknown coefficients are coupled derivatives.
The positional derivatives could be obtained by simulating the low-speed wind tunnel test, and the rotational derivatives are obtained by measuring the viscous hydrodynamic of the model at different rotational angular velocities. The number of coupled derivatives is large, and the viscous hydrodynamic subjected to the rotational manipulator is measured when β is not 0. After extensive experiments, the viscous hydrodynamic coefficients were obtained by least-square regression analysis.
The calculation domain inlet is set as the velocity boundary, the outlet as the free-flow condition, and the calculation domain wall as the non-slip fixed wall.
According to the research of the turbulence model in the reference [13], the standard k-ω model has many advantages, such as good numerical stability, accurate solution of pressure gradient, low Reynolds number influence, compressibility effect, and shear flow diffusion. It has better adaptability in calculating the flow problem of the boundary layer separation. It is one of the most widely used turbulence models for the viscous hydrodynamic solution. Therefore, the standard k-ω model is used as the turbulence model in this paper. The equation of the kinetic energy k and the turbulent  The rotation process of the manipulator is realized by UDF and dynamic mesh technology in Fluent. Dynamic meshing is performed by using the spring smoothing model, which approximates connecting lines between the grid nodes as springs, and the position of nodes after smoothing is obtained by calculating force equilibrium equations between them. In the calculation process, meshes with a large aberration rate or huge change in size are grouped together to re-divide the partial meshes.

The Solution of Position Derivatives, Rotation Derivatives, and Coupling Derivatives
Solving the unknown coefficients in Equation (2) is the key to the research task, where Y v , N v , Z w , and M w are positional derivatives, Y r and N r are rotational derivatives, and the remaining unknown coefficients are coupled derivatives.
The positional derivatives could be obtained by simulating the low-speed wind tunnel test, and the rotational derivatives are obtained by measuring the viscous hydrodynamic of the model at different rotational angular velocities. The number of coupled derivatives is large, and the viscous hydrodynamic subjected to the rotational manipulator is measured when β is not 0. After extensive experiments, the viscous hydrodynamic coefficients were obtained by least-square regression analysis.
The calculation domain inlet is set as the velocity boundary, the outlet as the free-flow condition, and the calculation domain wall as the non-slip fixed wall.
According to the research of the turbulence model in the reference [13], the standard k-ω model has many advantages, such as good numerical stability, accurate solution of pressure gradient, low Reynolds number influence, compressibility effect, and shear flow diffusion. It has better adaptability in calculating the flow problem of the boundary layer separation. It is one of the most widely used turbulence models for the viscous hydrodynamic solution. Therefore, the standard k-ω model is used as the turbulence model in this paper. The equation of the kinetic energy k and the turbulent frequency ω are as shown in Equation (6). The determined model parameters are shown in Table 2, according to the reference [13]:

Result Analysis
Comparing the viscous hydrodynamics of the manipulator measured by adopting the basic and dense mesh, respectively, we conclude that the differences are in the range from 1.7% to 4.3%. The differences are sufficiently low and negligible.
Computations are executed in a 64-bit processor consist of CPU (Intel Core i5-8400 @ 2.80 GHz) and 8 GB accessible memory, and take 37 h with at least 40 iterations per time step.
Using the size of calculation domain set in Section 3.1, the flow velocity of the calculation domain is kept constant, the angle β between the manipulator and the incoming flow direction is adjusted, the viscous hydrodynamics at different angles are calculated, and then the viscous hydrodynamic position derivative is obtained by linear analysis. The position derivative calculation contents are shown in Table 3. Table 3. Solving cases of position derivatives.
Since the different degrees of α and β cause the manipulator to have different linear velocities in the OX, OY, and OZ directions, viscous hydrodynamics of the arm is measured in the horizontal plane XOZ and the vertical plane XOY, respectively. The values of the vertical force Y and the pitching moment N measured at the different linear velocity v in the OY direction are shown in Table 4. The position derivatives Yv and Nv are the first derivative coefficients of the linear velocity v, so the least squares curve fitting of the vertical force Y and the pitching moment N in viscous hydrodynamics to the linear velocity v is performed, as shown in Figures 5 and 6. The derivative value at the median 0 points of the line velocity v is taken as the position derivative.  Same as above, the values of the lateral force Z and the yaw moment M measured at the different linear velocity w in the OZ direction are shown in Table 5. The position derivatives Zw and Mw are the first derivative coefficients of the linear velocity w, so the least squares curve fitting of the viscous hydrodynamic lateral force Z and the yawing moment M to the linear velocity w is performed, as shown in Figures 7 and 8. The derivative at 0 point of the median line velocity w is taken as the position derivative. The position derivatives Yv and Nv are the first derivative coefficients of the linear velocity v, so the least squares curve fitting of the vertical force Y and the pitching moment N in viscous hydrodynamics to the linear velocity v is performed, as shown in Figures 5 and 6. The derivative value at the median 0 points of the line velocity v is taken as the position derivative.  Same as above, the values of the lateral force Z and the yaw moment M measured at the different linear velocity w in the OZ direction are shown in Table 5. The position derivatives Zw and Mw are the first derivative coefficients of the linear velocity w, so the least squares curve fitting of the viscous hydrodynamic lateral force Z and the yawing moment M to the linear velocity w is performed, as shown in Figures 7 and 8. The derivative at 0 point of the median line velocity w is taken as the position derivative. Same as above, the values of the lateral force Z and the yaw moment M measured at the different linear velocity w in the OZ direction are shown in Table 5. The position derivatives Z w and M w are the first derivative coefficients of the linear velocity w, so the least squares curve fitting of the viscous hydrodynamic lateral force Z and the yawing moment M to the linear velocity w is performed, as shown in Figures 7 and 8. The derivative at 0 point of the median line velocity w is taken as the position derivative.  The positional derivatives could be obtained as shown in Table 6 below: Table 6. Position derivatives of the manipulator in viscous hydrodynamics. Since the manipulator studied in this paper is a single DOF, only the viscous hydrodynamics are measured when it rotates around the OZ axis. UDF is applied to adjust the rotational velocity of the arm. The measurement conditions of the rotation derivatives are shown in Table 7.  This paper calculates the hydrodynamics of the rotating underwater manipulator when the flow is stationary. The arm rotates in the XOY plane. The angular velocity r is adjusted according to Table  7, and the instantaneous viscous hydrodynamics of the arm during the rotation from −10° to 10° around the OZ axis are measured. The measured vertical force Y and pitching moment N are shown in Table 8.   The positional derivatives could be obtained as shown in Table 6 below: Since the manipulator studied in this paper is a single DOF, only the viscous hydrodynamics are measured when it rotates around the OZ axis. UDF is applied to adjust the rotational velocity of the arm. The measurement conditions of the rotation derivatives are shown in Table 7. This paper calculates the hydrodynamics of the rotating underwater manipulator when the flow is stationary. The arm rotates in the XOY plane. The angular velocity r is adjusted according to Table  7, and the instantaneous viscous hydrodynamics of the arm during the rotation from −10° to 10° around the OZ axis are measured. The measured vertical force Y and pitching moment N are shown in Table 8.  The positional derivatives could be obtained as shown in Table 6 below: Since the manipulator studied in this paper is a single DOF, only the viscous hydrodynamics are measured when it rotates around the OZ axis. UDF is applied to adjust the rotational velocity of the arm. The measurement conditions of the rotation derivatives are shown in Table 7. This paper calculates the hydrodynamics of the rotating underwater manipulator when the flow is stationary. The arm rotates in the XOY plane. The angular velocity r is adjusted according to Table 7, and the instantaneous viscous hydrodynamics of the arm during the rotation from −10 • to 10 • around the OZ axis are measured. The measured vertical force Y and pitching moment N are shown in Table 8. The rotational derivatives Y r and N r are the first derivative coefficients of the angular velocity r, similar to the solution method of the position derivatives. The least squares curve fitting of the vertical force Y and the pitching moment N in viscous hydrodynamics to the angular velocity r is performed, as shown in Figures 9 and 10. The derivative at the median angular velocity r of 1.25 rad/s is taken as the rotation derivative. The rotational derivatives Yr and Nr are the first derivative coefficients of the angular velocity r, similar to the solution method of the position derivatives. The least squares curve fitting of the vertical force Y and the pitching moment N in viscous hydrodynamics to the angular velocity r is performed, as shown in Figures 9 and 10. The derivative at the median angular velocity r of 1.25 rad/s is taken as the rotation derivative.  The rotational derivatives could be obtained as shown in Table 9 below: The coupled derivative is the second-order viscous hydrodynamic coefficient of the arm under complex motion conditions. Through the transient iterative calculation of the motion of the manipulator at different α and β angles and various angular velocities, the viscous hydrodynamics of each motion state are measured, and many coupling derivatives are obtained by least-squares regression analysis. The coupling derivative solution conditions are shown in Table 10 below.  The rotational derivatives Yr and Nr are the first derivative coefficients of the angular velocity r, similar to the solution method of the position derivatives. The least squares curve fitting of the vertical force Y and the pitching moment N in viscous hydrodynamics to the angular velocity r is performed, as shown in Figures 9 and 10. The derivative at the median angular velocity r of 1.25 rad/s is taken as the rotation derivative.  The rotational derivatives could be obtained as shown in Table 9 below: The coupled derivative is the second-order viscous hydrodynamic coefficient of the arm under complex motion conditions. Through the transient iterative calculation of the motion of the manipulator at different α and β angles and various angular velocities, the viscous hydrodynamics of each motion state are measured, and many coupling derivatives are obtained by least-squares regression analysis. The coupling derivative solution conditions are shown in Table 10 below.  The rotational derivatives could be obtained as shown in Table 9 below: Table 9. Rotation derivatives of the manipulator in viscous hydrodynamics. The coupled derivative is the second-order viscous hydrodynamic coefficient of the arm under complex motion conditions. Through the transient iterative calculation of the motion of the manipulator at different α and β angles and various angular velocities, the viscous hydrodynamics of each motion state are measured, and many coupling derivatives are obtained by least-squares regression analysis. The coupling derivative solution conditions are shown in Table 10 below. The coupled derivatives are obtained by regression analysis, and the obtained coupled derivatives are processed without dimensioning. The values of the parameters are as shown in Table 11.

Conclusions
The hydrodynamics of underwater manipulators during operation are complex and difficult to predict. It is analyzed that components of the hydrodynamics include inertial hydrodynamics caused by acceleration and viscous hydrodynamics caused by friction. This paper takes viscous hydrodynamics as the research target.
According to the research on dynamics of AUV (Autonomous Underwater Vehicle), ROV(Remote Operated Vehicle) and ships in other references, the viscous hydrodynamic model of the manipulator is analyzed in the form of Taylor expansion, and the viscous hydrodynamics are measured by using a single-DOF manipulator to simulate the underwater motion, and the viscous hydrodynamic coefficients among the model are calculated via regression analysis. An accurate viscous hydrodynamic model is obtained to predict viscous hydrodynamics of manipulators during operation at any attitude.
This model is the basis for the feedforward control and is helpful to study control stability of underwater manipulators. We believe that although the simulations in this paper were performed for single-DOF manipulators, the modeling method may be extended for the manipulator with multi-DOF and more complicated shapes as is in the case of real underwater manipulators.

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