Study on Dynamic Behavior of Unmanned Surface Vehicle-Linked Unmanned Underwater Vehicle System for Underwater Exploration

This paper focuses on motion analysis of a coupled unmanned surface vehicle (USV)–umbilical cable (UC)–unmanned underwater vehicle (UUV) system to investigate the interaction behavior between the vehicles and the UC in the ocean environment. For this, a new dynamic modeling method for investigating a multi-body dynamics system of this coupling system is employed. Firstly, the structure and hardware composition of the proposed system are presented. The USV and UUV are modeled as rigid-body vehicles, and the flexible UC is discretized using the catenary equation. In order to solve the nonlinear coupled dynamics of the vehicles and flexible UC, the fourth-order Runge–Kutta numerical method is implemented. In modeling the flexible UC dynamics, the shooting method is applied to solve a two-point boundary value problem of the catenary equation. The interaction between the UC and the USV–UUV system is investigated through numerical simulations in the time domain. Through the computer simulation, the behavior of the coupled USV–UC–UUV system is analyzed for three situations which can occur. In particular, variation of the UC forces and moments at the tow points and the configuration of the UC in the water are investigated.


Introduction
The guidance and control of marine vessels is an area of focus within the research community. The use of marine vehicles is increasing rapidly within several fields, such as marine biology, seafloor do so, efficient dynamic models for subsystems were integrated into the total system. The USV is modeled using 3-degree-of-freedom (DOF) rigid body dynamics, while the motion of the UUV is analyzed in 6-DOF. For UC modeling, the catenary equation to conduct motion analysis of the AUV and cable coupling system is applied, and the shooting method is then used to solve the nonlinear finite differential equations in our numerical simulation scheme to obtain more reliable solving results. Finally, the equations of motion considering the motion correlation between each subsystem are described. Several simulations were carried out using the developed equations of motion, and the movement of the USV, as well as that of the UUV, was also observed according to disturbance generated by the UC. The rest of this paper is organized as follows: Section 2 introduces the structure and hardware composition of both USV and UUV systems. Section 3 explores the dynamic model of the coupled USV and UC system. Both kinematics and kinetics are highlighted in this section. Then, Section 4 describes the UC dynamics and how to apply the UC effects to the vehicles. Next, Section 5 provides an overview of the mathematical modeling of the UUV, considering the interaction forces of the UC. Section 6 presents the motion analysis of the integrated USV-UC-UUV system moving in a series of scenarios by numerical simulation. Finally, Section 7 provides the main contributions and conclusions of the work presented throughout this paper.

Control System of the USV
The control system of the USV consisted of a USV operation control system and a winch control system for controlling depth of the UUV. The system composition is shown in Figure 2.
The architecture of the USV control system is shown in Figure 3. An AHRS (attitude heading reference system) and GPS sensor were used to control USV position and attitude, as well as for navigation using the line-of-sight method.

Control System of the USV
The control system of the USV consisted of a USV operation control system and a winch control system for controlling depth of the UUV. The system composition is shown in Figure 2.
The control signal for the UUV came through the tether cable using serial communication. The winch system was used to control the depth of the UUV. The winch system was installed on the USV as shown in Figure 1.

Structure of the UUV
The shape of the UUV is shown in Figure 4. Three horizontal thrusters were installed at the UUV for surge, sway, and yaw motion control. The configuration of thrusters in the propulsion system of this UUV is shown in Figure 5. The UUV was designed as an over-actuated system with seven thrusters. The three vertical thrusters were used for heave, pitch, and roll motion, and the four horizontal thrusters were used for surge, sway, and yaw motion. The direction of the thrusters is defined as follows: the four horizontal thrusters were defined as positive where they made a positive contribution in the x-direction. The architecture of the USV control system is shown in Figure 3. An AHRS (attitude heading reference system) and GPS sensor were used to control USV position and attitude, as well as for navigation using the line-of-sight method.
Sensors 2020, 20, 1329 4 of 27 The control signal for the UUV came through the tether cable using serial communication. The winch system was used to control the depth of the UUV. The winch system was installed on the USV as shown in Figure 1.

Structure of the UUV
The shape of the UUV is shown in Figure 4. Three horizontal thrusters were installed at the UUV for surge, sway, and yaw motion control. The configuration of thrusters in the propulsion system of this UUV is shown in Figure 5. The UUV was designed as an over-actuated system with seven thrusters. The three vertical thrusters were used for heave, pitch, and roll motion, and the four horizontal thrusters were used for surge, sway, and yaw motion. The direction of the thrusters is defined as follows: the four horizontal thrusters were defined as positive where they made a positive contribution in the x-direction. The control signal for the UUV came through the tether cable using serial communication.
The winch system was used to control the depth of the UUV. The winch system was installed on the USV as shown in Figure 1.

Structure of the UUV
The shape of the UUV is shown in Figure 4. Three horizontal thrusters were installed at the UUV for surge, sway, and yaw motion control. The control signal for the UUV came through the tether cable using serial communication. The winch system was used to control the depth of the UUV. The winch system was installed on the USV as shown in Figure 1.

Structure of the UUV
The shape of the UUV is shown in Figure 4. Three horizontal thrusters were installed at the UUV for surge, sway, and yaw motion control. The configuration of thrusters in the propulsion system of this UUV is shown in Figure 5. The UUV was designed as an over-actuated system with seven thrusters. The three vertical thrusters were used for heave, pitch, and roll motion, and the four horizontal thrusters were used for surge, sway, and yaw motion. The direction of the thrusters is defined as follows: the four horizontal  The configuration of thrusters in the propulsion system of this UUV is shown in Figure 5. The UUV was designed as an over-actuated system with seven thrusters. The three vertical thrusters were used for heave, pitch, and roll motion, and the four horizontal thrusters were used for surge, sway, and yaw motion. The direction of the thrusters is defined as follows: the four horizontal thrusters were defined as positive where they made a positive contribution in the x-direction. Meanwhile, the three vertical thrusters corresponded with a positive z-direction contribution. To avoid water being flushed through the UUV, these thrusters were tilted slightly. Figure 6 shows the control structure for the UUV. Meanwhile, the three vertical thrusters corresponded with a positive z-direction contribution. To avoid water being flushed through the UUV, these thrusters were tilted slightly. Figure 6 shows the control structure for the UUV.

USV Dynamic Modeling
A mathematical model was required for designing the controller of the USV. For this reason, the numerical modeling of the USV is presented in this section. The ocean environment (wind, waves, and currents) which affects the USV dynamics was neglected in this paper.

Assumptions
To simplify the problem, the motion of the USV is described only in the horizontal plane in this paper. The motion variables of the USV in the body-fixed coordinate are shown in Figure 7. Some simplifications were made for computer simulations of the USV motion. These simplifications were Meanwhile, the three vertical thrusters corresponded with a positive z-direction contribution. To avoid water being flushed through the UUV, these thrusters were tilted slightly. Figure 6 shows the control structure for the UUV.

USV Dynamic Modeling
A mathematical model was required for designing the controller of the USV. For this reason, the numerical modeling of the USV is presented in this section. The ocean environment (wind, waves, and currents) which affects the USV dynamics was neglected in this paper.

Assumptions
To simplify the problem, the motion of the USV is described only in the horizontal plane in this paper. The motion variables of the USV in the body-fixed coordinate are shown in Figure 7. Some simplifications were made for computer simulations of the USV motion. These simplifications were

USV Dynamic Modeling
A mathematical model was required for designing the controller of the USV. For this reason, the numerical modeling of the USV is presented in this section. The ocean environment (wind, waves, and currents) which affects the USV dynamics was neglected in this paper.

Assumptions
To simplify the problem, the motion of the USV is described only in the horizontal plane in this paper. The motion variables of the USV in the body-fixed coordinate are shown in Figure 7. Some simplifications were made for computer simulations of the USV motion. These simplifications were as follows: • The motion of the USV in roll, pitch, and heave directions was neglected.

•
The USV had neutral buoyancy and the origin of the body-fixed coordinate was located at the center of mass.

•
The USV had three planes of symmetry.

•
The dynamic equations of the USV did not include the disturbance forces (waves, wind, and ocean currents).
Sensors 2020, 20, 1329 6 of 27  The motion of the USV in roll, pitch, and heave directions was neglected.  The USV had neutral buoyancy and the origin of the body-fixed coordinate was located at the center of mass.  The USV had three planes of symmetry.  The dynamic equations of the USV did not include the disturbance forces (waves, wind, and ocean currents).

Three-Coordinate Systems
In modeling the complete USV-UC-UUV system in this study, we defined three coordinate systems composed of the earth-fixed coordinate (XE, YE, ZE), the local cable coordinates along the UC (C1, C2, C3), and the vehicle-fixed coordinate (xb, yb, zb), as shown in Figure 8. Firstly, the earth-fixed coordinate was located at the mass center of the USV, which was defined with XE pointing to the northerly direction, YE pointing to the easterly direction, and ZE pointing to the earth. Next, the UC was divided into small rigid segments instead of a continuous non-rigid system. Each segment had a coordinate system (C1, C2, C3) with C1 tangent to the UC, and C3 in the plane of {XE, YE}, whereas C1 and C2 were orthogonal vectors applied to the whole length of the UC. Lastly, the vehicle-fixed coordinate (xb, yb, zb) was defined with respect to the UUV itself. The UUV reference coordinates are described in Figure 8.

Three-Coordinate Systems
In modeling the complete USV-UC-UUV system in this study, we defined three coordinate systems composed of the earth-fixed coordinate (X E , Y E , Z E ), the local cable coordinates along the UC (C 1 , C 2 , C 3 ), and the vehicle-fixed coordinate (x b , y b , z b ), as shown in Figure 8. Firstly, the earth-fixed coordinate was located at the mass center of the USV, which was defined with X E pointing to the northerly direction, Y E pointing to the easterly direction, and Z E pointing to the earth. Next, the UC was divided into small rigid segments instead of a continuous non-rigid system. Each segment had a coordinate system (C 1 , C 2 , C 3 ) with C 1 tangent to the UC, and C 3 in the plane of {X E , Y E }, whereas C 1 and C 2 were orthogonal vectors applied to the whole length of the UC. Lastly, the vehicle-fixed coordinate (x b , y b , z b ) was defined with respect to the UUV itself. The UUV reference coordinates are described in Figure 8.
northerly direction, YE pointing to the easterly direction, and ZE pointing to the earth. Next, the UC was divided into small rigid segments instead of a continuous non-rigid system. Each segment had a coordinate system (C1, C2, C3) with C1 tangent to the UC, and C3 in the plane of {XE, YE}, whereas C1 and C2 were orthogonal vectors applied to the whole length of the UC. Lastly, the vehicle-fixed coordinate (xb, yb, zb) was defined with respect to the UUV itself. The UUV reference coordinates are described in Figure 8.  With the defined coordinates, the rotation matrix for converting from the UUV coordinate (x b , y b , z b ) to the earth-fixed coordinate (X E , Y E , Z E ), with regard to Euler angles in Reference [16], can be expressed as follows: where φ, θ, and ψ are the roll, pitch, and heading angles of the vehicles, respectively. Moreover, the transformation from the UC coordinate (C 1 , C 2 , C 3 ) to the earth-fixed coordinate (X E , Y E , Z E ) can be represented in a simple matrix form as follows [13]: Considering Equations (1)-(3), the transformation from the UC coordinate (C 1 , C 2 , C 3 ) to the UUV coordinate (x b , y b , z b ) can be determined as follows: Note that the rotation matrix, R, is orthogonal and, hence, its inverse is equal to the transpose of its matrix, which is useful for easily converting from earth-fixed coordinates to body-fixed coordinates.

Mathematical Model
A model for a system is not only useful for formulating control algorithms, but also for performing the simulation. According to Reference [9], when deriving a control design model of a USV, it can be assumed that only the motions in the horizontal plane, which are surge motion, sway motion, and yaw motion, are of interest to reduce the model complexities.
To express the kinematics of a USV, two coordinates were defined as shown in Figure 7. The first coordinate EX E Y E Z E describes the earth-fixed coordinate while the second coordinate Bx B y B z B represents the body-fixed coordinate. Furthermore, using the notation of Society of Naval Architects and Marine Engineers (SNAME) (1950) [17], a USV operating in 3-DOF can be described by six motion variables. The first mode is the Sensors 2020, 20, 1329 8 of 29 position of the USV η = x y ψ T ∈ 3 , which is referred to as surge, sway, and yaw motions, and it describes the position of the USV in the horizontal plane with reference to earth-fixed coordinate EX E Y E Z E . The second mode is the velocity of the USV v = u v r T ∈ 3 , where u and v are the surge and sway velocities and r is the heave velocity with reference to body-fixed coordinate Bx B y B z B .
The relationship between position and orientation of the USV in the earth-fixed coordinate (EX E Y E Z E ) and the linear and angular velocities in the vehicle coordinate (Bx B y B z B ) is given as where the rotation matrix R(η) is defined as The dynamic model of the USV is based on the model by Fossen [18]. This model of the USV in 3-DOF is derived from the Newton-Euler motion equation as described in Reference [9].
where M ∈ 3×3 is the symmetric positive definite inertia matrix, C(v) ∈ 3×3 is the centripetal and Coriolis matrix, and D(v) ∈ 3×3 is the damping matrix; g(η) ∈ 3 represents the gravitational forces, and we assume that represents the control input, and τ cable = τ Cx , τ Cy , τ Cn ∈ 3 represents the cable forces and moments. Additionally, the matrices M, C(v), and D(v) are expressed as where the term m describes the dry mass of the USV, x g is the coordinate between the center of gravity and vehicle origin in the x-axis expressed in the body-fixed frame, and the term I z denotes the moments of inertia about the Bz B axis.
and Y |v|r , Y |r|r , N |v|v , N |r|v are the hydrodynamic coefficients of the USV.

Configuration of Thrusters
The USV was designed to have three thrusters for 3-DOF motion control as shown in Figure 9.
the moments of inertia about the B Bz axis.

Configuration of Thrusters
The USV was designed to have three thrusters for 3-DOF motion control as shown in Figure 9. With revolution speed ni, the thrust force is expressed as shown in Equation (11). With revolution speed n i , the thrust force is expressed as shown in Equation (11).
where ρ is the density of sea water, D P is the diameter of three thrusters located on the USV, and K T is the thrust coefficient. The relationship between local thrust force and body-fixed thrust force on the USV can be described as where the individual thrust force vector is , and the generalized force vector acting on the USV is τ c = F x F y M z T ∈ 3×1 . Moreover, the thrust allocation of the USV is expressed as where t p is the thrust deduction coefficient of each thruster, D L is the distance between the forward thruster (TH3) and two stern thrusters (TH1, TH2), and D s is the distance between the two-port thruster (TH1) and starboard thruster (TH2).

Cable Dynamic Modeling
The UC was used to connect the USV with the UUV and to supply power and communication. However, water resistance to the UC interferes with or restrains the movement of the USV and UUV. In this section, the interacting forces between the USV and UUV are analyzed. For this, the mathematical modeling of the UC motion in water is presented. For modeling of the UC motion in water, the catenary equations and shooting method are proposed.

Assumptions
In this paper, to analyze interaction forces of the UC between the USV and UUV, the following assumptions are proposed: • A continuous, inextensible, and flexible UC was used in this study.

•
The UC had no bending and torsional stiffness.

•
The length of the UC was constant L = 100 m.

•
The UC acted as the axial force, UC self-weight, and hydrodynamic drag forces.

•
The stress/strain of the UC was linear.

Mathematical Model
In this paper, when the UC was suspended between the USV and UUV, there were four kinds of forces acting on the UC including gravity forces, buoyancy forces, drag forces, and residual bottom tension, which are described in Figure 10.
In this paper, to analyze interaction forces of the UC between the USV and UUV, the following assumptions are proposed:  A continuous, inextensible, and flexible UC was used in this study.  The UC had no bending and torsional stiffness.  The length of the UC was constant L = 100 m.  The UC acted as the axial force, UC self-weight, and hydrodynamic drag forces.  The stress/strain of the UC was linear.

Mathematical Model
In this paper, when the UC was suspended between the USV and UUV, there were four kinds of forces acting on the UC including gravity forces, buoyancy forces, drag forces, and residual bottom tension, which are described in Figure 10.  External forces on the UC are caused by the environmental forces such as hydrodynamic drag and gravity. In this paper, it is assumed that the attached UC is a long slender pipe, and the drag force acts on the UC. Morison's equation is used to estimate the forces acting on the UC. Thus, the tangential and normal components of the drag force can be respectively expressed as follows: where ρ is the density of seawater, and d is the diameter of the cable; V t and V n are respectively the tangential and normal components of current velocity relative to the UC. C t and C n are respectively the tangential and normal drag coefficients. The UC has weight such that it satisfies the catenary equation. The static calculation of the shape and tension of an ideal cable was given by Reference [19].
For modeling the UC, the cable with two boundary conditions was divided into many segments along its length. Then, we defined s i to be the value of s in node I; thus, the segment i was the part of the cable where s i−1 < s ≤ s i . The reaction force f 0 in s = 0 is described as where m is the number of segments of the cable, w i i ∈ R 3 is the constant distributed force acting over segment i, and f i i ∈ R 3 is nodal force acting on the nodal point i. The strain of each cable segment ε can be expressed as the function of an infinitesimal cable segment ds and the stretched length dp.
Moreover, the axial strain can be described by the relationship between the axial tension T, Young's modulus of cable E, and the cross-sectional area of the cable A.

of 29
The cable tension T(s): [0, L] → R 3 is estimated from the expression where m is the number of cable segments divided along the length of the UC, and dp is obtained from Equation (17). The position of the cable expressed in coordinate i is defined as For simplification, we express The cable tension can be solved from Using Equations (17) and (18) along with the assumption of the cable which has a linear stress-strain relationship, the cable tension can be defined as Furthermore, the relationship between the direction of an unstretched length of the cable segment ds and a stretched cable segment dp can be obtained as By using Equations (23) and (24), Equation (25) can be rearranged as In order to solve Equation (26), the ordinary numerical integration method can be applied. Sagatun [19] presented the closed-form solution of this integral, which is where the sign ⊗ describes component-wise multiplication. The tension of each cable segment T(s) is obtained from Equation (23), where Using the assumption that cable segment 1 is at the point 0, this means that r(0) = r0, and the continuities between the segments have to be fulfilled. Then, we have The integration constants can be obtained as

Boundary Conditions
In order to solve the governing equation for the UC, two boundary conditions at both ends of the UC were applied. The response of the UC ends was given by the connections to the USV and UUV. In this case, the boundary conditions were placed at both ends of the UC (the upper end attached to the USV and the lower end connected to the UUV). Due to this, the USV motion was on the water surface; thus, the first boundary condition P USV = P USV (t) as a function of time represents the position of the USV during the motion. Similarly, the second boundary condition at the point connected to the UUV P UUV = P UUV (t) describes the position of the UUV, also considered a function of time. Thus, the top and bottom boundary conditions were where x usv , y usv , z usv and x uuv , y uuv , z uuv represent the positions of the USV and the UUV during their motion, respectively.

Cable Effects
Note that the tension of the UC, T, caused the forces and moments acting on both vehicles (USV and UUV) to vary with time because the positions of vehicles changed with time when the UC moved in the water. The equations below which were mentioned by Reference [8], were used to represent the UC effects on both vehicle dynamics, where the UC forces and moments are expressed in the vehicle-fixed frame.
where → r c = (x c , y c , z c ) denotes the location of the tow points at the vehicles (USV and UUV), expressed in the vehicle frames.

UUV Dynamic Modeling
A detailed dynamic model of the UUV was used to estimate the vehicle behavior in various situations. In this section, the 6-DOF equation of the motions of the UUV with seven thrusters are described based on the physics of an over-actuated UUV and its actuators.

Assumptions
The dynamic model of the UUV was quite complex and needed many parameters; therefore, the following assumptions were used to simplify the model.

•
The UUV was fairly symmetrical about its three planes.

•
The center of buoyancy of the UUV was located on the geometric symmetry plane. • There were no environmental disturbances acting on the UUV.

•
The UUV was considered as a rigid body; thus, there were no bending and geometrical deformations.

•
The hydrodynamic coefficients of the UUV were not variable.

Mathematical Model
Following standard practice [18], a 6-DOF nonlinear model of a UUV was expressed using the earth-fixed coordinate and a body-fixed coordinate as shown in Figure 11. The body-fixed coordinate O-xyz was attached to the center of the gravity of the vehicle. The motion of the body-fixed coordinate is described relative to the earth-fixed coordinate E-XYZ.  The transformation matrix J relates the motion of the UUV in the body-fixed coordinate to the earth-fixed coordinate. Thus, the kinematic equation for the UUV expressed using generalized coordinates in 6-DOF is given as The rotation matrix between the body-fixed coordinate and the earth-fixed coordinate related to Euler angles is given by   22 1 sin tan cos tan 0 cos sin 0 sin / cos cos / cos where the notations s(.) = sin (.), c(.) = cos (.), and t(.) = tan (.) are used for notational brevity. Notice The notation used throughout this paper is based on the 6-DOF representation for the UUV, given by Reference [17].
where η 1 ∈ 3×1 is the linear position of the UUV, and η 2 ∈ 3×1 is the vector of Euler angles. Both η 1 and η 2 are defined in the earth-fixed coordinate E-XYZ. Meanwhile, v 1 ∈ 3×1 denotes the translational velocities in surge, sway, and heave motions of the UUV, and v 1 ∈ 3×1 denotes the rotational velocities in roll, pitch, and yaw motions in the body-fixed coordinate O-xyz. Finally, the vector τ ∈ 3×1 describes the generalized forces and moments acting on the UUV in the body-fixed coordinate O-xyz, with τ 1 ∈ 3×1 corresponding to the forces along x-, y-, and z-axes, while τ 2 ∈ 3×1 corresponds to the moments about the x-, y-, and z-axes.
The transformation matrix J relates the motion of the UUV in the body-fixed coordinate to the earth-fixed coordinate. Thus, the kinematic equation for the UUV expressed using generalized coordinates in 6-DOF is given as Sensors 2020, 20, 1329
where M ∈ 6×6 is an inertial matrix of UUV.C(v) ∈ 6×6 is a centripetal force and Coriolis matrix. D(v) ∈ 6×6 is a hydrodynamic damping matrix. G(η) ∈ 6×1 is a gravity and buoyancy term, τ th ∈ 6×1 represents the propulsion forces and moments acting on the UUV, and τ cable ∈ 6×1 denotes the UC forces and moments. Moreover, the aforementioned matrices are described as follows: −I yz q − I xz p + I zz r + N . r r I yz r + I xy p − I yy q − M . q q I yz q + I xz p − I zz r − N . r r 0 −I xz r − I xy q + I xx p + K . p p −I yz r − I xy p + I yy q + M . q q I xz r + I xy q − I xx p − K . p p 0 All symbols of variables used in the above equations can be explained as follows: m denotes the mass of the UUV, O G = x G y G z G T is the center of gravity of the UUV, I xx , I yy , and I zz are the moments of inertia of the UUV about the Bx B , By B , and Bz B axes, I xy = I yx , I xz = I zx , and I yz = I zy are the products of inertia, W is the weight of the UUV body expressed in the earth-fixed coordinate, and B is the submerged buoyancy force expressed in the earth-fixed coordinate; x b , y b , and z b are the center of buoyancy of the UUV expressed in the body-fixed coordinate. The partial derivative coefficients (X . u , , the components of linear drag (X u , Y v , Z w , K p , N r ), and quadratic drag coefficients (X u|u| , Y v|v| , Z w|w| , K p|p| , M q|q| , N r|r| ) are the hydrodynamic coefficients which can be directly or indirectly obtained in advance by practical experiments.
Alternatively, the dynamic model equation of the UUV in earth-fixed coordinate E-XYZ can also be obtained using the kinematic transformations, ..
to eliminate v and . v in Equation (40). Hence, the following earth-fixed vector expression of dynamic model can be expressed as where

Configuration of Thrusters
In this paper, the UUV was equipped with seven thrusters as shown in Figure 12. The UUV did not have rudders; thus, its motion was only affected by the thrusters. Four horizontal thrusters T 1 , T 2 , T 3 , and T 4 , which were installed at the bow and the stern part with inclined angle α, were responsible for the motions along the horizontal plane. Meanwhile, three vertical thrusters T 5 , T 6 , and T 7 were responsible for the motions along the vertical plane.

Configuration of Thrusters
In this paper, the UUV was equipped with seven thrusters as shown in Figure 12. The UUV did not have rudders; thus, its motion was only affected by the thrusters. Four horizontal thrusters T1, T2, T3, and T4, which were installed at the bow and the stern part with inclined angle  , were responsible for the motions along the horizontal plane. Meanwhile, three vertical thrusters T5, T6, and T7 were responsible for the motions along the vertical plane. Since the UUV was equipped with seven thrusters and controlled in 6-DOF, the relationship between the required force in each DOF and the forces was The thruster configuration matrix T can be expressed as follows: where ls , lf, lr, ds, df, dr are the lengths of the arm that create momentum in roll, pitch, and yaw, and 0 30   is the angle thruster which located in the xy-plane. Since the UUV was equipped with seven thrusters and controlled in 6-DOF, the relationship between the required force in each DOF and the forces was in which τ th is the desired force in the different DOF, T is the thruster configuration matrix, and τ u is the desired force of each thruster. τ th and τ u are defined as The thruster configuration matrix T can be expressed as follows: where l s , l f , l r , d s , d f , d r are the lengths of the arm that create momentum in roll, pitch, and yaw, and α = 30 0 is the angle thruster which located in the xy-plane.

Simulation Procedure
As the nonlinear dynamic equations are complex and difficult to solve analytically, the numerical simulation approach was adopted to simulate the motion of the USV-UC-UUV coupling system. For this, a model-based simulation process is proposed as shown in Figure 13. After the establishment of the respective dynamic equations of the coupled USV-UC-UUV system, the dynamic equations were further discretized for numerical simulation. Thus, several iteration methods were applied to solve the discretized equations.
As the nonlinear dynamic equations are complex and difficult to solve analytically, the numerical simulation approach was adopted to simulate the motion of the USV-UC-UUV coupling system. For this, a model-based simulation process is proposed as shown in Figure 13. After the establishment of the respective dynamic equations of the coupled USV-UC-UUV system, the dynamic equations were further discretized for numerical simulation. Thus, several iteration methods were applied to solve the discretized equations. In order to solve the nonlinear differential equations, different iteration methods were applied to the respective USV, UUV, and UC dynamic equations due to their different dynamic characteristics. Usually, the Runge-Kutta method is well known to solve the differential equations of USV and UUV dynamics, while the shooting method is used to solve the partial differential equations of UC dynamics. For the dynamic behavior of the UC, the partial differential Equation (26) under the boundary conditions Equations (31) and (32) could be solved numerically using the shooting method. For this, the UC was divided into n segments (or nodes) equally, such that the shooting method could be used by discretizing the UC dynamics using both the time segment ( t  ) and the cable length segment ( s  ). Note that the length of the connecting UC was fixed when the USV and UUV moved forward in the water. In this case, the length of each segment of the UC was fixed in the finite difference method. The finite difference in Equation (26) involves n cable nodes, and, according to the boundary conditions in Equations (31) and (32), there are 3n + (6 + 3) cable node variables in total. Because the UUV is a rigid body, 12 motion states were used to describe its dynamics. The USV had 6 motion states to describe the dynamics. Thus, for the combined USV-UC-UUV system, there were in total 6 + (3n + 6 + 3) + 12 dynamic equations to solve the 6 + (3n + 6 + 3) + 12 motion states, as described below.
To calculate the motion states in Equation (53), the initial condition of the system ( In order to solve the nonlinear differential equations, different iteration methods were applied to the respective USV, UUV, and UC dynamic equations due to their different dynamic characteristics. Usually, the Runge-Kutta method is well known to solve the differential equations of USV and UUV dynamics, while the shooting method is used to solve the partial differential equations of UC dynamics. For the dynamic behavior of the UC, the partial differential Equation (26) under the boundary conditions Equations (31) and (32) could be solved numerically using the shooting method. For this, the UC was divided into n segments (or nodes) equally, such that the shooting method could be used by discretizing the UC dynamics using both the time segment (δt) and the cable length segment (δs). Note that the length of the connecting UC was fixed when the USV and UUV moved forward in the water. In this case, the length of each segment of the UC was fixed in the finite difference method.
The finite difference in Equation (26) involves n cable nodes, and, according to the boundary conditions in Equations (31) and (32), there are 3n + (6 + 3) cable node variables in total. Because the UUV is a rigid body, 12 motion states were used to describe its dynamics. The USV had 6 motion states to describe the dynamics. Thus, for the combined USV-UC-UUV system, there were in total 6 + (3n + 6 + 3) + 12 dynamic equations to solve the 6 + (3n + 6 + 3) + 12 motion states, as described below.
To calculate the motion states in Equation (53), the initial condition of the system (U 0 whole ) should be provided for further calculation. Then, the previous iteration result (U k whole ) can be used as the initial guess for the next iteration (U k+1 whole ). In this paper, an integrated simulation model was made to represent the dynamics of the complete USV-UC-UUV system in the time domain. The simulation model in Matlab-Simulink consisted of the USV dynamics, the UC dynamics, and the UUV model as shown in Figure 14.
The parameters for simulation of the complete USV-UC-UUV are shown in Table 1. The detailed parameters of the USV and UUV used in this study were given in References [20][21][22], respectively. It is challenging to simulate the complete USV-UC-UUV system as it is a poorly damped system and the dynamics are highly nonlinear. Based on the simulation scheme introduced in Figure 13, the model-based motion simulations are further implemented in the later sections. The three following simulations were performed: Sensors 2020, 20, 1329 17 of 27 the initial guess for the next iteration ( 1 k whole U  ). In this paper, an integrated simulation model was made to represent the dynamics of the complete USV-UC-UUV system in the time domain. The simulation model in Matlab-Simulink consisted of the USV dynamics, the UC dynamics, and the UUV model as shown in Figure 14. The parameters for simulation of the complete USV-UC-UUV are shown in Table 1. The detailed parameters of the USV and UUV used in this study were given in References [20,21] and [22], respectively.
It is challenging to simulate the complete USV-UC-UUV system as it is a poorly damped system and the dynamics are highly nonlinear. Based on the simulation scheme introduced in Figure  13

Simulation 1
In the first simulation, the motion of the coupled UUV and UC system in the horizontal plane was studied. The USV was assumed to stay at a position by using a perfect dynamic positioning controller such that the USV-induced motion was ignored in this case. At this time, the UUV performed a turning motion as shown in Figure 15

Simulation 1
In the first simulation, the motion of the coupled UUV and UC system in the horizontal plane was studied. The USV was assumed to stay at a position by using a perfect dynamic positioning controller such that the USV-induced motion was ignored in this case. At this time, the UUV performed a turning motion as shown in Figure 15. In this simulation, the simulation duration was 35 s, with a sampling time of 0.01 s. The initial position of the UUV in the earth-fixed coordinate was η UUV = (70, 0, 50, 0, 0, 0), while the initial velocity of the UUV expressed in body-fixed coordinate was ν UUV = (0.2, 0, 0, 0, 0, 0). On the contrary, the USV was assumed at the stable position (0, 0, 0) by using a strong controller.   A shown in Figure 15, to force the UUV into a clockwise turning motion, the four horizontal thrusters T1, T2, T3, and T4 were applied using different values such as T1 = 11.5 N, T2 = 12 N, T3 = −11.5 N, and T4 = −12 N. Meanwhile, the three thrusters T5, T6, T7 were set to 0 N to achieve a pure turning motion in this simulation.

Length of USV
In general, the UUV drag increases as the length of the UC increases. Thus, the maximum affordable UC length for the UUV needs to be designed according to the power capacity of the UUV. The trajectories of the UUV runs in terms of the turning motion without the UC and with the UC are shown in Figure 16. −11.5 N, and T4 = −12 N. Meanwhile, the three thrusters T5, T6, T7 were set to 0 N to achieve a pure turning motion in this simulation.
In general, the UUV drag increases as the length of the UC increases. Thus, the maximum affordable UC length for the UUV needs to be designed according to the power capacity of the UUV. The trajectories of the UUV runs in terms of the turning motion without the UC and with the UC are shown in Figure 16. The variation of the force and moment of the UC acting upon the turning motion of the UUV is shown in Figure 17. It shows that the oscillatory heave force Fcz and the roll moment Mcx from the UC caused the oscillators to achieve the heave velocity and the roll motion of the UUV. In addition, the surge force Fxz initially increased and then decreased due to the negative value right after the UUV turning to the left side, while the sway force Fcy decreased at the beginning and then increased to a positive value when the UUV underwent the turning motion.
The effects of the UC on the position and orientation of the UUV are clearly shown in Figure 18. Obviously, when operating the UUV and the UC coupling system in the turning motion, the state variables of the UUV were significantly affected.
For the velocities of the UUV, the UUV initially moved straight for about 10 s, then turned leftward, and finally moved backward. In this motion, the surge velocity increased from the initial value 0.2 m/s and then decreased to a constant speed 0.37 m/s, while the sway speed of the UUV slightly decreased from 0 m/s to a negative small steady value of about −0.12 m/s, as shown in Figure  19. For the heave velocity and roll motion, oscillations are shown. Moreover, the pitch motion is also affected because of the UC with a negative decrease negatively initially before increasing to a positive angle. However, the results also show that the interaction of the UC and the UUV did not affect the surge motion, sway motion, and yaw motion of the UUV. The variation of the force and moment of the UC acting upon the turning motion of the UUV is shown in Figure 17. It shows that the oscillatory heave force F cz and the roll moment M cx from the UC caused the oscillators to achieve the heave velocity and the roll motion of the UUV. In addition, the surge force F xz initially increased and then decreased due to the negative value right after the UUV turning to the left side, while the sway force F cy decreased at the beginning and then increased to a positive value when the UUV underwent the turning motion. The effects of the UC on the position and orientation of the UUV are clearly shown in Figure 18. Obviously, when operating the UUV and the UC coupling system in the turning motion, the state variables of the UUV were significantly affected. For the velocities of the UUV, the UUV initially moved straight for about 10 s, then turned leftward, and finally moved backward. In this motion, the surge velocity increased from the initial value 0.2 m/s and then decreased to a constant speed 0.37 m/s, while the sway speed of the UUV slightly decreased from 0 m/s to a negative small steady value of about −0.12 m/s, as shown in Figure 19. For the heave velocity and roll motion, oscillations are shown. Moreover, the pitch motion is also affected because of the UC with a negative decrease negatively initially before increasing to a positive angle. However, the results also show that the interaction of the UC and the UUV did not affect the surge motion, sway motion, and yaw motion of the UUV.

Simulation 2
In the second simulation, the motion of the coupled USV and UC system in the horizontal plane was analyzed with time varied. Unlike the first simulation, the 3-DOF USV firstly moved straight in the forward direction for about 15 s from the origin of the earth-fixed coordinate η USV = (0, 0, 0, 0, 0, 0), then turned left afterward, while the UUV was assumed to maintain position (50, 0, 70). Furthermore, the initial velocity of the USV was ν USV = (0.1, 0, 0, 0, 0, 0) and the simulation duration was 25 s, with a sampling time of 0.01 s in this case. Figure 20 shows the thruster directions of the USV during its turning motion. In order to simulate the turning motion of the USV, two stern thrusters TH1 and TH2 were applied with different values of TH1 = 10 N and TH2 = 9 N, while the bow thruster TH3 was set to 0 N.  Figure 20 shows the thruster directions of the USV during its turning motion. In order to simulate the turning motion of the USV, two stern thrusters TH1 and TH2 were applied with different values of TH1 = 10 N and TH2 = 9 N, while the bow thruster TH3 was set to 0 N.

Figure 20.
Thruster directions of the USV in turning motion. Figure 21 shows the trajectory of the USV during the turning motion, either with or without the UC. In this study, the proposed system consisted of the USV connected to the UUV using the UC. Thus, the forces generated by the motion of the UC could affect the motion of the combined vehicles. The forces and moments of the UC affecting the USV during the turning motion are presented in Figure 22. It shows that the exerting force of the UC on the USV seemed significant when the UUV underwent the turning motion. However, all the state variables of USV (surge, sway, and yaw motions) were not affected much, as shown in Figures 23 and 24. The corresponding phenomena can be explained by the size of the USV being relatively large compared to that of the UUV, while the movement of the UC and the UUV had a small effect on the behavior of the USV.   Figure 21 shows the trajectory of the USV during the turning motion, either with or without the UC. In this study, the proposed system consisted of the USV connected to the UUV using the UC. Thus, the forces generated by the motion of the UC could affect the motion of the combined vehicles. The forces and moments of the UC affecting the USV during the turning motion are presented in Figure 22. It shows that the exerting force of the UC on the USV seemed significant when the UUV underwent the turning motion. However, all the state variables of USV (surge, sway, and yaw motions) were not affected much, as shown in Figures 23 and 24. The corresponding phenomena can be explained by the size of the USV being relatively large compared to that of the UUV, while the movement of the UC and the UUV had a small effect on the behavior of the USV.  Figure 20 shows the thruster directions of the USV during its turning motion. In order to simulate the turning motion of the USV, two stern thrusters TH1 and TH2 were applied with different values of TH1 = 10 N and TH2 = 9 N, while the bow thruster TH3 was set to 0 N.

Figure 20.
Thruster directions of the USV in turning motion. Figure 21 shows the trajectory of the USV during the turning motion, either with or without the UC. In this study, the proposed system consisted of the USV connected to the UUV using the UC. Thus, the forces generated by the motion of the UC could affect the motion of the combined vehicles. The forces and moments of the UC affecting the USV during the turning motion are presented in Figure 22. It shows that the exerting force of the UC on the USV seemed significant when the UUV underwent the turning motion. However, all the state variables of USV (surge, sway, and yaw motions) were not affected much, as shown in Figures 23 and 24. The corresponding phenomena can be explained by the size of the USV being relatively large compared to that of the UUV, while the movement of the UC and the UUV had a small effect on the behavior of the USV.

Simulation 3
In this simulation, the dynamic behaviors of the complete USV-UC-UUV system were analyzed to show the effects of the UC on USV and UUV motions. The USV went straight in the forward direction, while the UUV undertook a sideward motion. Similarly to simulation 1, a time interval of 0.01 s and a simulation time of 35 s were set. The initial positions of the connecting point to the UUV (position of the UUV in the earth-fixed coordinate) were set at (70 m, 0 m, 50 m) while the end point at the free surface near the USV was assumed set at (0 m, 0 m, 0 m). The initial velocities of both USV and UUV in the body-fixed coordinate were ν USV = (0.1, 0, 0, 0, 0, 0) and ν UUV = (0.2, 0, 0, 0, 0, 0), respectively. Figure 25 shows the thruster directions of the USV and the UUV during their motions. In order to simulate the forward motion of the USV, the thrust forces of the two stern thrusters TH1 and TH2 were set to 10 N, while the bow thruster TH3 was set to 0 N. For the case of sideward motion of the UUV, the thrust forces of T1 and T3 were set to 10 N, while the thrusters T2 and T4 were set to −10 N, as shown in Figure 25. Meanwhile, thrusters T5, T6, and T7 were set to 0 N to have sideward motion only.

Simulation 3
In this simulation, the dynamic behaviors of the complete USV-UC-UUV system were analyzed to show the effects of the UC on USV and UUV motions. The USV went straight in the forward direction, while the UUV undertook a sideward motion. Similarly to simulation 1, a time  Figure 25 shows the thruster directions of the USV and the UUV during their motions. In order to simulate the forward motion of the USV, the thrust forces of the two stern thrusters TH1 and TH2 were set to 10 N, while the bow thruster TH3 was set to 0 N. For the case of sideward motion of the UUV, the thrust forces of T1 and T3 were set to 10 N, while the thrusters T2 and T4 were set to −10 N, as shown in Figure 25. Meanwhile, thrusters T5, T6, and T7 were set to 0 N to have sideward motion only. For the above case, Figure 26 shows the trajectories of the USV and UUV without the UC and with the UC interacting force. According to the simulation results, in order to maintain the depth of the UUV during operation of the complete USV-UC-UUV system, the speed difference between the USV and the UUV should be confined. The UC tension causes additional drag on the USV and UUV motions, which affect the behavior of the vehicles. The variations of the UC force and moment at the tow points (upper and lower points) are presented for the USV and UUV, as shown in Figures 27  and 28, respectively. Comparing Figure 27 with Figure 28, it is shown that the force of the UC on the UUV is more significant than on the USV.
Unlike simulations 1 and 2, by considering the different motion operations on the horizontal plane of the USV and the UUV, the analysis on the dynamic behavior of both the USV and the UUV affected by the interaction of the UC is presented in this simulation. The effects of the UC on the USV and the UUV are presented in Figures 29-31. For USV motion, the USV was set to move in a straight line with the velocity Vusv = 0.25 m/s; after 25 s, the velocity increased to 0.8 m/s as shown in Figure  29. The results also show that the USV moved faster with the connected UC because the pushed forces of the surge force, Fx, increased gradually.
For the sideward motion of the UUV, the UUV moved leftward for 7.5 m in 35 s with a sway velocity of about 0.23 m/s, as shown in Figures 30 and 31. The results in these figures show that the UUV shifted a little in surge motion while moving sideward due to the initial surge velocity of the For the above case, Figure 26 shows the trajectories of the USV and UUV without the UC and with the UC interacting force. According to the simulation results, in order to maintain the depth of the UUV during operation of the complete USV-UC-UUV system, the speed difference between the USV and the UUV should be confined. The UC tension causes additional drag on the USV and UUV motions, which affect the behavior of the vehicles. The variations of the UC force and moment at the tow points (upper and lower points) are presented for the USV and UUV, as shown in Figures 27  and 28, respectively. Comparing Figure 27 with Figure 28, it is shown that the force of the UC on the UUV is more significant than on the USV. significant effect of the UC on the depth, roll, and pitch motion modes. In particular, both the pitch and heave motions of the UUV were significantly oscillatory because of the heave force Fcz. The UUV moved in the sway direction with up and down oscillations and changed the pitch angle of the UUV during simulation time. This occurred because the heave velocity w and the pitch angle θ of the UUV regularly oscillated. Furthermore, the results of the surge, sway, and yaw motions of the UUV showed similar motion compared to previous ones without UC.  and heave motions of the UUV were significantly oscillatory because of the heave force Fcz. The UUV moved in the sway direction with up and down oscillations and changed the pitch angle of the UUV during simulation time. This occurred because the heave velocity w and the pitch angle θ of the UUV regularly oscillated. Furthermore, the results of the surge, sway, and yaw motions of the UUV showed similar motion compared to previous ones without UC.   Unlike simulations 1 and 2, by considering the different motion operations on the horizontal plane of the USV and the UUV, the analysis on the dynamic behavior of both the USV and the UUV affected by the interaction of the UC is presented in this simulation. The effects of the UC on the USV and the UUV are presented in Figures 29-31. For USV motion, the USV was set to move in a straight line with the velocity V usv = 0.25 m/s; after 25 s, the velocity increased to 0.8 m/s as shown in Figure 29. The results also show that the USV moved faster with the connected UC because the pushed forces of the surge force, F x , increased gradually.    In summary, the simulations results show that the effects of the connected UC to the vehicles were big, especially for the UUV motions. The simulations could be very helpful when designing the capacity of the thrusters for the UUV and USV.

Conclusions
In this paper, a new mathematical modeling of a coupling system of a USV connected to a UUV using a UC was presented. To do so, analysis of the dynamics was firstly performed on each system, and the total coupled system dynamics was then studied. Because the UC connected the USV and the UUV, the dynamic equation of the UC was derived using the catenary equation. To analyze the behavior of the UC, the shooting method was applied. To demonstrate the application of the proposed equations, model-based motion simulations of the coupling system were performed. Computer simulations were conducted to analyze the interacting forces of the UC with the USV and UUV systems. In the simulation, the maneuvering behavior of the USV with the UC, the maneuvering behavior of the UUV with the UC, and the maneuvering behavior of the coupled USV-UC-UUV system were investigated and their results were discussed. Moreover, the variation of the UC forces and moments at the tow points and the configuration of the UC were analyzed.
The simulation results revealed that the UC significantly affected the motion of both the USV and the UUV in all cases (especially the UUV). The results also showed that the UC tension caused additional drag forces on the USV and UUV, and they affected the motion of both vehicles, while the variation of the configuration of the UC could result in the vehicle getting tangled, especially when the coupled system moves in currents. Based on suitable assumptions, the numerical model developed in the paper could numerically simulate the motions of the vehicles with the UC in the ocean environment. It is believed that the simulation results may provide useful guidance and reference for real USV-UC-UUV systems in design and operation. Using the results of the analysis on the interaction forces between the systems, it would be helpful to design the capacity of actuators for the USV or UUV.
Future work can be extended by taking into account the motion of the USV-UC-UUV system under the effect of various speeds of underwater currents, with verification via water tests. Furthermore, an analysis on the controller will be designed to reduce the effects of the UC.   Figure 30 shows the significant effect of the UC on the depth, roll, and pitch motion modes. In particular, both the pitch and heave motions of the UUV were significantly oscillatory because of the heave force F cz . The UUV moved in the sway direction with up and down oscillations and changed the pitch angle of the UUV during simulation time. This occurred because the heave velocity w and the pitch angle θ of the UUV regularly oscillated. Furthermore, the results of the surge, sway, and yaw motions of the UUV showed similar motion compared to previous ones without UC.
In summary, the simulations results show that the effects of the connected UC to the vehicles were big, especially for the UUV motions. The simulations could be very helpful when designing the capacity of the thrusters for the UUV and USV.

Conclusions
In this paper, a new mathematical modeling of a coupling system of a USV connected to a UUV using a UC was presented. To do so, analysis of the dynamics was firstly performed on each system, and the total coupled system dynamics was then studied. Because the UC connected the USV and the UUV, the dynamic equation of the UC was derived using the catenary equation. To analyze the behavior of the UC, the shooting method was applied. To demonstrate the application of the proposed equations, model-based motion simulations of the coupling system were performed. Computer simulations were conducted to analyze the interacting forces of the UC with the USV and UUV systems. In the simulation, the maneuvering behavior of the USV with the UC, the maneuvering behavior of the UUV with the UC, and the maneuvering behavior of the coupled USV-UC-UUV system were investigated and their results were discussed. Moreover, the variation of the UC forces and moments at the tow points and the configuration of the UC were analyzed.
The simulation results revealed that the UC significantly affected the motion of both the USV and the UUV in all cases (especially the UUV). The results also showed that the UC tension caused additional drag forces on the USV and UUV, and they affected the motion of both vehicles, while the variation of the configuration of the UC could result in the vehicle getting tangled, especially when the coupled system moves in currents. Based on suitable assumptions, the numerical model developed in the paper could numerically simulate the motions of the vehicles with the UC in the ocean environment. It is believed that the simulation results may provide useful guidance and reference for real USV-UC-UUV systems in design and operation. Using the results of the analysis on the