Hydrodynamic Modelling for a Transportation System of Two Unmanned Underwater Vehicles: Semi-Empirical, Numerical and Experimental Analyses

: Underwater transportation is an essential approach for scientiﬁc exploration, maritime construction and military operations. Determining the hydrodynamic coefﬁcients for a complex underwater transportation system comprising multiple vehicles is challenging. Here, the suitability of a quick and less costly semi-empirical approach to obtain the hydrodynamic coefﬁcients for a complex transportation system comprising two Unmanned Underwater Vehicles (UUVs) is investigated, where the interaction effects between UUVs are assumed to be negligible. The drag results were veriﬁed by Computational Fluid Dynamics (CFD) analysis at the steady state. The semi-empirical results agree with CFD in heave and sway; however, they were overpredicted in surge due to ignoring the wake effects. Furthermore, experiments were performed for the validation of the time-domain motion simulations with semi-empirical and CFD results. The simulations which were performed with the CFD drags were close to the experiments. The semi-empirical approach could be relied on once a correction parameter is included to account for the interactive effect between multiple UUVs. Overall, this work makes a contribution by deriving a semi-empirical approach for the dynamic and controlling system of dual UUVs, with CFD and experiments applied to ascertain its accuracy and potential improvement.


Introduction
Underwater transportation is an important requirement for the scientific, commercial and military sectors as it can place an underwater structure at a precise location [1][2][3]. Unmanned Underwater Vehicles (UUVs) can be used in different systems to undertake underwater transportation based on the connection between the vehicles and the payload, i.e., Rigid Connection Transportation System (RCTS) and Flexible Connection Transportation System (FCTS). In RCTS, UUVs are connected to the payload via solid links as used in [4][5][6][7]. The whole system is considered to be a single body. Alternatively, in FCTS the vehicles are connected to the payload via flexible links, as adopted in [8][9][10][11][12], where there is a relative motion between the vehicles and the payload. In this paper, hydrodynamic modelling is studied for an RCTS system operating underwater.
The development of the dynamic model for an underwater transportation system is vital to design stable, robust, high-accuracy and low-power-consumption control systems. The dynamic model consists of inertia, hydrostatic and hydrodynamic terms. The inertia and hydrostatic parameters can be directly measured with accuracy as they mainly depend on the physical system. The hydrodynamic parameters describe the opposition to the motion of the system due to the medium in which it is moving. Obtaining an accurately tuned dynamic model for an ROV and designing an appropriate control system are challenging due to the uncertainties in the hydrodynamic parameters [13,14]. The difficulty increases Hydrodynamic models have been developed for a single UUV based on a semiempirical approach. Prestero [41] developed a semi-empirical hydrodynamic model for a torpedo-shaped Autonomous Underwater Vehicle (AUV). The empirical formula was used to work out the drag term in the axial direction, whereas strip theory was used for the drag terms in other direction separately for the cylindrical hull and fins. The coefficients were used from the empirical data. The roll drag was estimated by considering only the fins, as the major portion of rolling resistance comes from fins. For the axial added mass term, the vehicle hull was approximated by an ellipsoid for which the major axis is half the vehicle length and the minor axis is half the vehicle diameter. The added mass terms in other directions were worked out using the strip theory separately for the cylindrical and cruciform hull cross-sections. The rolling added mass was calculated only for the fins as the smooth hull sections only generate a small added mass in roll which was neglected. The empirical formulas were used to work out the body and fin lift coefficients in all directions. Humphreys [44] worked out the added mass terms in detail for each part of a torpedo-shaped AUV using the semi-empirical approach. In [43], the lift, drag and pitching moment coefficients were calculated for a bare hull AUV using the semi-empirical equations which were derived in the literature, and CFD analysis was carried out to verify them. In [42], analytical and semi-empirical (ASE) methods were used to approximate the hydrodynamic parameters for AUV. Eidsvik [40] calculated the hydrodynamic parameters for a box-shaped Remotely Operated Vehicle (ROV) using the DNV-standard [48], which is based on the Applied Fluid Dynamics book by Blevins [49]. In his approach, the ROV was assumed to be a solid prism where two to three sides are equal. This assumption could be true for some specific ROVs such as SF-30k, AC-ROV 100, Seabotics LBV600-6 and Videoray PRO-4, where the height and width are approximately the same. Appropriate scaling coefficient was included to account for the penetrating flow through the ROV. The translational added mass terms were calculated using the potential flow data. However, the 3D data for the rotational added mass terms could not be found. Therefore, they were approximated by the 2D data and strip theory. The translational drag terms were calculated using the empirical. Again, the 3D data for the rotational drag parameters were not available. Therefore, they were approximated for small angles where rotational motion can be transformed to translational motion.
The calculation of hydrodynamic parameters for complex-shaped underwater vehicles or the underwater systems which consist of multiple UUVs is lacking in the literature. Therefore, hydrodynamic models are required to be developed for such systems. Moreover, the semi-empirical approach needs to be improved to reduce hydrodynamic uncertainties. The contribution of this work is twofold: (1) an effective semi-empirical approach is introduced to estimate the hydrodynamic coefficients for a system of multiple UUVs; (2) the method is verified against CFD results and validated with experimental measurements.
In this paper, the hydrodynamic parameters were not calculated in the experiments. Instead, the system was run in the towing tank and its motion response was recorded in the time domain. On the other hand, the dynamic model was developed for the system with the calculated hydrodynamic parameters from, respectively, the semi-empirical approach and CFD. The time-domain motion simulations were performed on the dynamic model. The motion response of the simulations was compared to the experiments for validation. A similar approach was used by Singh [50] for an underwater glider in which the hydrodynamic parameters were calculated using CFD analysis. These coefficients were used in the motion equations proposed by Leonard and Graver [51] in the vertical plane to obtain a simulation model. Subsequently, the motion response in simulations was compared to the experiments for validation.
The rest of the paper is arranged as follows. Section 2 describes the fabrication of the rigid connection transportation system (RCTS), which was tested in the towing tank and velocities were measured. This system, comprising two UUVs, represents the case study analysed in this work. In Section 3, the hydrodynamic parameters are calculated using the semi-empirical approach and the drag forces are verified with CFD analysis. A dynamic model is also developed to run time-domain motion simulations. In Section 4, the drag forces, which are obtained by semi-empirical and CFD methods, are compared and the difference is discussed to analyse the potential deficiency of the semi-empirical approach. The time-domain motion simulations are carried out with the coefficients computed with both methods, which are compared to the experiments for validation. Finally, concluding remarks are given in Section 5.

Experiment
The hydrodynamic modelling is performed on an available system that can be tested in the towing tank for validation.

Fabrication
For the dual UUVs transportation system, which is rigidly connected to the payload, two identical modified Seaperch UUVs [52,53] were fabricated using PVC pipes which were connected through elbows and T joints, as shown in Figure 1. The sizes of pipe, elbow and T-joint are shown in Table 1.
J. Mar. Sci. Eng. 2021, 9, 500 4 of 23 velocities were measured. This system, comprising two UUVs, represents the case study analysed in this work. In Section 3, the hydrodynamic parameters are calculated using the semi-empirical approach and the drag forces are verified with CFD analysis. A dynamic model is also developed to run time-domain motion simulations. In Section 4, the drag forces, which are obtained by semi-empirical and CFD methods, are compared and the difference is discussed to analyse the potential deficiency of the semi-empirical approach. The time-domain motion simulations are carried out with the coefficients computed with both methods, which are compared to the experiments for validation. Finally, concluding remarks are given in Section 5.

Experiment
The hydrodynamic modelling is performed on an available system that can be tested in the towing tank for validation.

Fabrication
For the dual UUVs transportation system, which is rigidly connected to the payload, two identical modified Seaperch UUVs [52,53] were fabricated using PVC pipes which were connected through elbows and T joints, as shown in Figure 1. The sizes of pipe, elbow and T-joint are shown in Table 1.  This design of Seaperch was selected as it is highly stable with a COG well below its COB. Four thrusters were attached to each Seaperch: two to move the vehicle in the vertical direction and two to move the UUV in the axial direction. Each vehicle is underactuated as no thrusters are applied in the transverse direction, so that only surge, heave, pitch and yaw can be controlled.
After the fabrication, the mass and volume of each UUV were measured to determine the weight and buoyancy. The COG of each vehicle was also measured, whilst the COB of the vehicle was calculated using the COB and volume of each part of the vehicle.
In practice, the underwater vehicles are made slightly positively buoyant so that they rise to the sea surface in case of an emergency [54]. However, for research studies, especially for designing the control systems, the test models are made to be neutrally buoyant to ease analysis [40,[55][56][57][58]. Therefore, each Seaperch UUV was made neutrally buoyant, and after attaching the manipulators and payload to the UUVs, the whole system was made neutrally buoyant by attaching small patches of buoyancy sheet. The power cables were made neutrally buoyant to nullify their impact on the hydrostatic parameters. After fabrication and connections, the whole system is shown in Figure 2.
This design of Seaperch was selected as it is highly stable with a COG well below its COB. Four thrusters were attached to each Seaperch: two to move the vehicle in the vertical direction and two to move the UUV in the axial direction. Each vehicle is underactuated as no thrusters are applied in the transverse direction, so that only surge, heave, pitch and yaw can be controlled.
After the fabrication, the mass and volume of each UUV were measured to determine the weight and buoyancy. The COG of each vehicle was also measured, whilst the COB of the vehicle was calculated using the COB and volume of each part of the vehicle.
In practice, the underwater vehicles are made slightly positively buoyant so that they rise to the sea surface in case of an emergency [54]. However, for research studies, especially for designing the control systems, the test models are made to be neutrally buoyant to ease analysis [40,[55][56][57][58]. Therefore, each Seaperch UUV was made neutrally buoyant, and after attaching the manipulators and payload to the UUVs, the whole system was made neutrally buoyant by attaching small patches of buoyancy sheet. The power cables were made neutrally buoyant to nullify their impact on the hydrostatic parameters. After fabrication and connections, the whole system is shown in Figure 2. A MATLAB program was developed to provide input to the thrusters of the system using the Arduino hardware. An initial estimate of the dynamic stability of the system was obtained from the Routh Hurwitz criteria [59]. However, after the fabrication of the system, several initial tests were performed to ensure the dynamic stability of the system. The parameters of the fabricated RCTS are shown in Table 2, including instrument precision.

Velocity Measurement
In tests conducted in the UCL towing tank, the system was given a constant thrust input via the Arduino hardware and the time taken to achieve the marked distances on the towing tank glass boundary was measured. The acceleration stage was assumed to last a short time from simulations and was thus neglected. The velocities, which were measured in surge and heave, are shown in Table 3. An uncertainty analysis was undertaken to quantify the systematic and random errors, with the uncertainty being driven low by repeating the tests multiple times.

Hydrodynamic Modelling
In this section, the hydrodynamic forces are calculated for the fabricated RCTS using the semi-empirical approach, with each piece of the Seaperch UUV, manipulator and payload being analysed separately and then summed together. Additionally, in this section, a CFD model was built and used to obtain the drag forces.

Geometry
A 3D model of the system was produced in the SolidWorks software [60], as shown in Figure 3. The plan, profile and front views are shown in Figures 4-6, respectively. A MATLAB program was developed to provide input to the thrusters of the system using the Arduino hardware. An initial estimate of the dynamic stability of the system was obtained from the Routh Hurwitz criteria [59]. However, after the fabrication of the system, several initial tests were performed to ensure the dynamic stability of the system. The parameters of the fabricated RCTS are shown in Table 2, including instrument precision.

Velocity Measurement
In tests conducted in the UCL towing tank, the system was given a constant thrust input via the Arduino hardware and the time taken to achieve the marked distances on the towing tank glass boundary was measured. The acceleration stage was assumed to last a short time from simulations and was thus neglected. The velocities, which were measured in surge and heave, are shown in Table 3. An uncertainty analysis was undertaken to quantify the systematic and random errors, with the uncertainty being driven low by repeating the tests multiple times.

Hydrodynamic Modelling
In this section, the hydrodynamic forces are calculated for the fabricated RCTS using the semi-empirical approach, with each piece of the Seaperch UUV, manipulator and payload being analysed separately and then summed together. Additionally, in this section, a CFD model was built and used to obtain the drag forces.

Geometry
A 3D model of the system was produced in the SolidWorks software [60], as shown in

Assumptions
The following assumptions were made when calculating the hydrodynamic parameters for the studied system using the semi-empirical and CFD methods: The system is a rigid body, i.e., there is no relative motion between the mass particles of the system. This is generally true as the parts on the system were not found to bend during testing. 2.
The system's mass and its distribution do not change during the motion. This is a valid assumption and ensured during testing. 3.
The system is operating below base wave depth, i.e., the depth at which the sea wave effects are negligible. The baseline depth is equal to λ/2 where λ is equal to the wavelength [61]. This was ensured in the experiments. 4.
The system only has translational motion, i.e., in surge, sway and heave. This was ensured during experiments as only the respective thrusters were used for motion in each direction and no disturbances were included. 5.
The vehicle does not experience a memory effect, i.e., the vehicle does not pass through its own wake. 6.
The vehicle does not experience the ocean current. The calm water in the towing tank would have negligible water currents. 7.
Interaction with the seabed or other underwater bodies is neglected. This was ensured during the experiments. 8.
The flow interaction amongst the parts of the system is ignored in the semi-empirical approach. This is something of concern and could be the main reason to cause an error. This effect is accounted for in CFD. 9.
The blades of the thrusters are ignored. They are small in size; therefore, their effect on the overall hydrodynamic results is small. 10. The elbows and T-joints are considered to have pipes of the same inner and outer diameters as of the other pipes, and the length is the same as of the respective elbow or T-joint. This would have a small impact as the extended portion of elbows and T-joint is only 6.5 mm. 11. The impact of buoyancy sheet patches on the hydrodynamic parameters is also ignored.
Although the system is underactuated and cannot move in pure sway, the added mass and drag forces were calculated in all three translational directions, i.e., surge, sway and heave. The velocity in sway is assumed to be the same as in surge.

Semi-Empirical Approach
The semi-empirical methods which are mostly applied in the literature to calculate the hydrodynamic parameters for a single UUV consider the whole vehicle, either a cylinder with fins or a solid prism. However, they are only true for specific types and shapes of UUVs. For the complex shaped UUV or a system of multiple UUVs, there is a lack of research in the literature to work out the hydrodynamic parameters. Therefore, in this study, a different approach of estimating the hydrodynamic parameters for an underwater vehicle is proposed, in which the system is cut down to simple individual parts. The parameters are calculated for each part and combined to obtain the net effect.
The hydrodynamic parameters were calculated for one Seaperch UUV and one manipulator plate, which were then doubled to obtain the effect of two Seaperch UUVs and two manipulators. The payload plate was evaluated separately. For the Seaperch UUV, the hydrodynamic parameters were evaluated separately for each pipe and thruster motor. They were then summed to obtain the net effect in each direction. The system parts are shown in Table 4.

Added Mass
The added mass is the inertia added to the vehicle as it displaces the surrounding water during acceleration or deceleration. The calculation of added mass does not take into account the viscous effects. The added mass part of Morison's equation for a solid body is given as [47]: where A i is the added mass term due to acceleration i or, in other words, the derivative of the added mass w.r.t acceleration i. For instance, the added mass terms in surge, sway and heave are X . u , Y . v and Z . w , respectively. ρ is the density of water (in this case, freshwater), C a is the added mass coefficient and V R is the reference volume. Axial thruster motors 04 6.
Payload plate 01 The added mass terms were calculated separately for each pipe and thruster motor of the Seaperch UUV and for the manipulator and payload plates using the data which are based on the DNV standards [48]. The detailed calculations are shown in Appendix A. After calculations, the added mass terms for each part of the system are shown in Table 5. They were summed together to obtain the total added mass terms in each surge, sway and heave, as shown in Table 6. The impact of added mass on the overall motion response would be less significant due to the short acceleration duration. The added mass terms are added to the total mass of the system to obtain the net inertia terms, which are then multiplied by the accelerations to obtain the inertial forces.

Drag
The resistance to the motion of the vehicle due to the viscosity of water is called drag. The drag part of Morison's equation for a solid body is given as [47]: where B i is the drag term due to velocity i. For instance, the drag terms in surge, sway and heave are X u|u| , Y v|v| and Z w|w| , respectively. C D/ f is the drag coefficient where D or f depends on the form or frictional drag, which further depends on the shape of the part in front of the flow and the flow velocity. The drag coefficients were also calculated from the data based on DNV standards given in [48]. A c/s is the cross-sectional or surface area depending on the form or frictional drag. Initially, it is important to know whether the flow over a body predominantly exerts form or frictional drag.
Similar to the added mass, the drag terms were calculated separately for each part of the Seaperch UUVs, manipulators and payload. Initially, it was worked by Reynold's number (Re) whether the flow over a part body predominantly exerts form or frictional drag.
The detailed calculations are shown in Appendix B. The total drag terms for the system in surge, sway and heave are shown in Table 7. The drag forces in surge, sway and heave were obtained by multiplying the drag terms by the square of the velocities in each direction.

CFD Analysis
The computational domain was built for the studied transportation system introduced in Section 2, using the STAR-CCM+ CFD software [62]. A three-dimensional cuboid computational domain was established, and the transportation system is fixed in the middle of the domain, as shown in Figure 7. The domain size is sufficiently large to avoid the system feeling any boundaries, in line with the experiments. The boundary conditions consist of a velocity inlet and a pressure outlet to propagate a constant flow against the device. The other four boundaries were defined as zero-gradient to model far fields. The domain is filled with water, and the water was defined as flowing with a constant velocity (U water ) against the system. Thus, a relative velocity exists between the system and water, where U water indicates the advancing speed of the system in calm water (U water = U), which is known in surge, sway and heave. This approach is similar to previous work studying the hydrodynamic performance of objects in water [63,64]. The direction of the water flow is set according to different motions of the system, as shown in Figure 8, in which the corresponding velocity fields at steady state are presented.
The solution of the fluid domain was obtained by solving the Reynolds-averaged Navier-Stokes (RANS) equations for an incompressible Newtonian fluid: where v is the time-averaged velocity, v is the velocity fluctuation, ρ is the fluid density, p denotes the time-averaged pressure, τ = µ [∇v+ (∇v) T ] is the viscous stress term, µ is the dynamic viscosity and g is gravitational acceleration set at 9.81 m/s 2 . The density of the water was set at 1000 kg/m 3 and the dynamic viscosity was 8.90 × 10 −4 N·s/m 2 . Since the RANS equations have been adopted to account for the turbulent effects, a turbulence model needs to be applied; here, the Shear Stress Transport (SST) k − ω model was adopted to close the equations.
J. Mar. Sci. Eng. 2021, 9, 500 10 of 23 is set according to different motions of the system, as shown in Figure 8, in which the corresponding velocity fields at steady state are presented.    The drag force is the steady-state force from fluid on the transportation system against the motion. This can be calculated as the integration of pressure and viscous forces on the system's surface depending on the direction of motion.
The mesh was refined around the device and for the nearby flow field. Following the ITTC guideline [65], the mesh was globally scaled by a factor of √ 2, which resulted in three mesh sets of, respectively, 1.56, 2.2 and 3.1 million cells. The corresponding drag forces in heave are 0.75, 0.80 and 0.80 N, which means the drag achieved convergence at 2.2 million cells, and further increasing the cells to 3.1 million could not alter the results. Therefore, the mesh set with 2.2 million cells was applied for further studies to save computational costs. The time step size was set based on the Courant number equalling to one. The applied CFD is a standard method and the numerical setup follows mature guidelines of the field [65], which has been extensively validated to be highly accurate, while this specific CFD model will be further validated by an experimental model introduced later.

Dynamic Modelling
The dynamic model is developed for the transportation system to run the time-domain motion simulations. Fossen's approach [66] is used in the development of the dynamic model in which the position and orientation of the system are defined in the earth-fixed frame (EFF), whereas velocities, forces and moments are defined in the body-fixed frame (BFF). Though the power cables were made neutrally buoyant, their inertia would have an impact on the overall inertia of the system. However, the inertia of the power cables was ignored in the system's dynamic model due to the difficulty of taking their effect on the centre of origin of the system (CO). Moreover, assumptions 1 to 7 in Section 3.2 were also made in the development of the dynamic model for the studied transportation system.
The kinematics provides the change in the position and orientation of the system, which is written in the vectorial form as: where J is the transformation matrix to convert velocity vector (v) in the earth-fixed frame (EFF), which is equated to the change in position and orientation vector . η . The kinetics provides the change in the velocity terms, which can be written in the vectorial form as [2]: where: Only the translational motions are investigated for the system in this study. The same is also ensured during the experiments. Therefore, the rigid body mass matrix (M RB ) consists of only the mass terms, given as: The mass of the system (m) was measured, as shown in Table 2. The added mass matrix is shown in Equation (8). The added mass terms which were used in the timedomain motion simulations for both the semi-empirical and drag methods were the ones calculated by the semi-empirical method. They were ignored in the calculations using the CFD method and would have less impact as the acceleration lasts a short time.
The damping matrix contains the terms which were calculated either by the semiempirical or CFD method. The separate motion simulations would take place with the damping results of semi-empirical and CFD methods. The damping matrix is written as: The Coriolis and centripetal terms are due to the rotation of the body-fixed frame (BFF) about the earth-fixed frame (EFF). Due to the consideration of no rotational motion, C RB (v) and C A (v) become zero. g(η) is the vector of hydrostatic forces and moments, which is obtained by the difference between weight and buoyancy. This is given as: τ is the vector of thrust forces and moments of the system. The thrust vector of each thruster (τ i ) can be calculated from [66]: where f x , f y and f z are the thrust forces applied in surge, sway and heave by each thruster, l x , l y , l z is the position of each thruster w.r.t the centre of the combined body (CO), which is taken at the centre of the payload.
Due to the consideration of only the translational force effects, τ i for each thruster is written as: The combined thrust vector (τ) can be written as the product of the thrust allocation matrix (T a ) and the vector of thrust forces of the thrusters (f), given as: The thrust allocation matrices for the two Seaperch UUVs in the system are shown as: All these matrices were incorporated in MATLAB [67] and the 4th order Runge-Kutta method was applied to simulate the motion response of the system over time.

Results and Discussion
For the experimental velocities, the corresponding CFD drag forces, as well as those calculated by the semi-empirical method, are shown in Table 8.

Comparison between the Semi-Empirical and CFD Predictions
It can be seen that the CFD and the semi-empirical values of heave and sway drag forces agree very well. For surge, CFD predicted a considerably lower value of drag than the semi-empirical method. This is because the semi-empirical method considers the parts of the system as encountering the flow directly and separately. Hence, the assumption that the semi-empirical method ignores the interaction between the individual components of the two UUVs may not be appropriate. In Figure 8b, it can be seen that the aft of the system stays in the wake of the forward side in surge, which significantly reduces the drag of the aft parts, though this was ignored in the semi-empirical equation. For heave and sway, the interaction effect of the front part on the back part is not notable, thus the semi-empirical equation reasonably predicted the drag forces.
Overall, the CFD analysis proved the general accuracy of the semi-empirical method. For surge, a correction coefficient needs to be added to the semi-empirical equation to calibrate the drag prediction. The value of the correction parameter should be a function of velocity and could be derived through further CFD tests. Once the calibration is done, the semi-empirical method is expected to be reliable for the prediction of the UUV's resistance. Despite its higher accuracy, CFD cannot provide immediate assessment when a new set of input conditions is given, limiting its application for real-time purposes. Therefore, using CFD to derive the parameters for the semi-empirical method is a superior solution to enable real-time UUV control.

Validation Against Experiments
The time-domain motion simulations were carried out separately for the semi-empirical and CFD drag coefficients and they were then compared to the experimental velocity results in surge and heave. Figure 9 compares the heave results of the system which were obtained during simulations and experiment. The experimental results assume zero initial acceleration and the immediate achievement of steady-state conditions. It can be seen that the simulation results which were obtained using the semi-empirical and CFD drag parameters are close to the experimental result. These findings show that the simulations using both the semi-empirical and CFD drag terms have well predicted the motion response in heave. Nonetheless, it can be seen that the heave response shows a slight shift for either the CFD or semi-empirical method. There could be a small impact of ignoring the blades of thruster and buoyancy sheet patches in the hydrodynamic calculations and considering the elbows and T-joints having the same diameters as pipes (as can be seen between Figures 2 and 3). Therefore, the deviation may result from the fact that the dynamic model ignores the inertia of the power cables.
semi-empirical method. There could be a small impact of ignoring the blades of thruster and buoyancy sheet patches in the hydrodynamic calculations and considering the elbows and T-joints having the same diameters as pipes (as can be seen between Figures 2 and 3). Therefore, the deviation may result from the fact that the dynamic model ignores the inertia of the power cables.  Figure 10 shows the surge response of the system. A substantial difference can be seen between the results of experiments and motion simulations using the semi-empirical drag terms. Conversely, when the CFD drag terms were used in the simulations, the obtained surge response is much closer to the experimental results. This justifies the assumption of the limited impact of the neutrally buoyant tether on the system's dynamics.   Figure 10 shows the surge response of the system. A substantial difference can be seen between the results of experiments and motion simulations using the semi-empirical drag terms. Conversely, when the CFD drag terms were used in the simulations, the obtained surge response is much closer to the experimental results. This justifies the assumption of the limited impact of the neutrally buoyant tether on the system's dynamics.
semi-empirical method. There could be a small impact of ignoring the blades of thruster and buoyancy sheet patches in the hydrodynamic calculations and considering the elbows and T-joints having the same diameters as pipes (as can be seen between Figures 2 and 3). Therefore, the deviation may result from the fact that the dynamic model ignores the inertia of the power cables.  Figure 10 shows the surge response of the system. A substantial difference can be seen between the results of experiments and motion simulations using the semi-empirical drag terms. Conversely, when the CFD drag terms were used in the simulations, the obtained surge response is much closer to the experimental results. This justifies the assumption of the limited impact of the neutrally buoyant tether on the system's dynamics.  An estimate of the time and cost involved in the three methods is shown in Table 9. The whole experimental procedure had the highest time involvement due to several steps involved in the fabrication, programming and tank tests, as shown in Section 2. On the other hand, CFD analysis was carried out in less time compared to experiments, which involved the development of a 3D CAD model in SolidWorks, building a 3D computational model in STAR-CCM+ [62], CFD mesh convergence study and obtaining the drag forces.
For the semi-empirical method, the practical system or the 3D model of the system was not required. Therefore, this method was completed in comparatively less time. For the cost analysis, the material costs were only involved in the fabrication and attachments of the physical system in the experiments. Other costs were estimated based on the cost of the man-hours. The cost of experiments is the highest as this involves the cost of fabrication and attachments, the manpower cost of fabrication, attachments, programming and tank tests. For the CFD and semi-empirical methods, only the manpower cost was involved, which is higher for CFD compared to the semi-empirical approach.
The time and cost requirements would increase with the complexity of the transportation system. However, the rise would be lower for the semi-empirical approach compared to the experiments and CFD.
Although the semi-empirical method is less accurate, it still offers the superior solution for the design of dynamic models of transportation systems, as its much lower computational cost and time-efficiency enable the treatment of multiple configurations with ease. Additionally, the performance of the semi-empirical method can be improved by applying corrections from the CFD solutions.

Future Work
The added mass terms used in the validation study were calculated by the potential flow theory based on the DNV standard. The potential flow theory is proved to be accurate for finding added mass terms for underwater vehicles in the translational directions. For instance, the added mass terms calculated using WAMIT potential flow software closely matched the terms found using free decay pendulum tests [68,69]. However, due to the spatial combination of multiple parts of various geometric shapes in the studied transportation system, the pressure gradient of added mass is likely to be affected by the arrangement of these parts in the combined system. For instance, accelerated flow over one part blocked by another would exhibit a pressure gradient different from the stand-alone rigid bodies of specific shapes. Therefore, a further experimental study is required for the accelerated motion of the system. This can be performed by attaching an Inertial Measurement Unit (IMU), which has accelerometers to measure linear acceleration and gyroscopes to measure angular velocities [70]. A commonly used IMU for underwater vehicles is MPU6050, which contains a triple-axis accelerometer and triple axis gyroscope [71,72]. After attaching MPU6050 to the transportation system, serial communication can be set to obtain the real-time acceleration reading in MATLAB via Arduino hardware. The linear acceleration profile can be integrated to obtain the velocity response and double integrated to obtain the displacement response of the system over time. The displacement response during the acceleration phase can be included in the experimental graph for the validation of the simulation response. This, in turn, validates the semi-empirical added mass terms which were used in the dynamic model for simulations.

Conclusions
The hydrodynamic parameters were estimated using a semi-empirical approach for a fabricated transportation system of two Seaperch UUVs. The hydrodynamic coefficients were worked out for each part of the system separately and they were summed together to obtain the net value in each direction. The semi-empirical drag results were then verified with a CFD analysis. It was found that the heave and sway drag forces are well estimated by the semi-empirical method, whereas the surge drag is overestimated. This is because the wake region generated by the front of the system, which affects the drag force on the aft bodies, was ignored in the semi-empirical method. Therefore, a corrective function needs to be derived and included in the semi-empirical equation to account for the interactions.
The semi-empirical approach can be improved by including a corrective function in the semi-empirical equation to account for the interactions. Once the semi-empirical equation is refined, it has the potential to facilitate the hydrodynamic assessment of multiple UUV systems. Although experiments and CFD are proved to be generally more accurate, they are relatively slow to complete and also costly. Thus, the rapid and less costly semi-empirical calculations are useful to develop a powerful tool in the design and controller optimisation of an underwater vehicle before it is built and tested. Added mass terms were calculated for the pipes and thruster motors of Seaperch UUV. 30.05 cm and 26.59 cm pipes: The flow passes over the circumference of the 30.05 cm and 26.59 cm pipes in surge and heave, i.e., perpendicular to the pipe. Therefore, the added mass terms in surge and heave were calculated using the data for the circular cylinder, as shown in Table A1. Table A1. Added mass data for the flow perpendicular to the circular cylinder (reproduced from [48] with permission from DNV, 2021).

Circular Cylinder
J. Mar. Sci. Eng. 2021, 9, 500 17 of 23 Table A1. Added mass data for the flow perpendicular to the circular cylinder (reproduced from [48] with permission from DNV, 2021). Conversely, in sway, the flow occurs horizontally over the 30.05 and 26.59 cm pipes, i.e., parallel to the pipe. However, the data are not available for the cylinder in [20]; therefore, the data for the parallel flow over a square prism were used, as shown in Table A2. Conversely, in sway, the flow occurs horizontally over the 30.05 and 26.59 cm pipes, i.e., parallel to the pipe. However, the data are not available for the cylinder in [20]; therefore, the data for the parallel flow over a square prism were used, as shown in Table A2. Table A2. Added mass coefficients for the flow parallel to the square prism (reproduced from [48], with permission from DNV, 2021).

Body Shape Direction of Motion b/a C A V R
Square Prism ∞ 1.00 Conversely, in sway, the flow occurs horizontally over the 30.05 and 26.59 cm pipes, i.e., parallel to the pipe. However, the data are not available for the cylinder in [20]; therefore, the data for the parallel flow over a square prism were used, as shown in Table A2. Table A2. Added mass coefficients for the flow parallel to the square prism (reproduced from [48], with permission from DNV, 2021). The flow parallel to a circular cylinder can be converted to the flow parallel to a square section using:

Body Shape Direction of Motion
where R is the radius of the circular cross-section of the cylinder, which is equal to 0.825 cm. Therefore, becomes 1.46 cm. 36.05 cm pipes: The flow over 36.05 cm pipes is perpendicular in sway and heave and parallel in surge. Therefore, the added mass coefficients for sway and heave were worked out using the data from Table A1 and for surge from Table A2. 13.55 cm pipes: Due to the position and orientation of 13.55 cm pipes, the flow passes over the circumference of the pipe in surge and sway, whereas it is parallel in heave. Therefore, the added mass coefficients for surge and sway were calculated using the data from Table A1 and for heave from Table A2.
Thruster motors: The added mass terms were worked out for the two axial and two vertical thruster motors. The flow passes over the circumference of the axial thrusters in sway and heave, whereas applying horizontally in surge. Therefore, data from Table A1 were used for calculating the added mass for sway and heave and from Table A2 for surge. For the vertical thruster motors, the flow passes over the circumference of the pipe in surge and sway, whereas it is horizontal in heave. Therefore, the added mass terms were worked out accordingly from Tables A1 and A2.

Appendix A.2. Added Mass Terms for Manipulator and Payload Plates
The manipulators and payload are the rectangular plates. The added mass terms for the manipulators and payload in heave were worked out using the added mass coefficients from Table A3. The flow parallel to a circular cylinder can be converted to the flow parallel to a square section using: a = √ πR, (A1) where R is the radius of the circular cross-section of the cylinder, which is equal to 0.825 cm. Therefore, a becomes 1.46 cm. 36.05 cm pipes: The flow over 36.05 cm pipes is perpendicular in sway and heave and parallel in surge. Therefore, the added mass coefficients for sway and heave were worked out using the data from Table A1 and for surge from Table A2. 13.55 cm pipes: Due to the position and orientation of 13.55 cm pipes, the flow passes over the circumference of the pipe in surge and sway, whereas it is parallel in heave. Therefore, the added mass coefficients for surge and sway were calculated using the data from Table A1 and for heave from Table A2.
Thruster motors: The added mass terms were worked out for the two axial and two vertical thruster motors. The flow passes over the circumference of the axial thrusters in sway and heave, whereas applying horizontally in surge. Therefore, data from Table A1 were used for calculating the added mass for sway and heave and from Table A2 for surge. For the vertical thruster motors, the flow passes over the circumference of the pipe in surge and sway, whereas it is horizontal in heave. Therefore, the added mass terms were worked out accordingly from Tables A1 and A2.

Appendix A.2. Added Mass Terms for Manipulator and Payload Plates
The manipulators and payload are the rectangular plates. The added mass terms for the manipulators and payload in heave were worked out using the added mass coefficients from Table A3.
The added mass terms in surge and sway were worked out by integrating the 2D added mass over the thickness. The 2D added mass was estimated using the coefficients from Table A4. The added mass terms in surge and sway were worked out by integrating the 2D added mass over the thickness. The 2D added mass was estimated using the coefficients from Table A4.
The drag term was calculated for each pipe on the Seaperch in surge, sway and heave. In Equation (A2), is the cross-sectional area of the pipe and is the drag coefficient of the pipe, which depends on the Reynolds number ( ).
To calculate the drag terms, was first worked out. This helps in deciding the data which are required to obtain the drag coefficient. The is calculated as: where is the density of fresh water, is the dynamic viscosity, is the diameter and is the length of the pipe. or depends on whether the flow is perpendicular or parallel to the pipe.
is the flow velocity. 30 The drag term B i was calculated for each pipe on the Seaperch in surge, sway and heave. In Equation (A2), A c is the cross-sectional area of the pipe and C D is the drag coefficient of the pipe, which depends on the Reynolds number (R e ).
To calculate the drag terms, R e was first worked out. This helps in deciding the data which are required to obtain the drag coefficient. The R e is calculated as: where ρ is the density of fresh water, µ is the dynamic viscosity, D is the diameter and L is the length of the pipe. D or L depends on whether the flow is perpendicular or parallel to the pipe. V is the flow velocity. 30.05 cm and 26.6 cm pipes: The R e for the flow in surge and heave was calculated from: where D is equal to 1.65 cm, V is equal to 0.167 m/sec in surge and 0.144 m/sec in heave.
The R e in both surge and heave came out to be below 10 5 . Therefore, the drag coefficients were calculated from the data in Table A5 for R e < 10 5 .
where is equal to 1.65 cm, is equal to 0.167 m/sec in surge and 0.144 m/sec in heave. The in both surge and heave came out to be below 10 . Therefore, the drag coefficients were calculated from the data in Table A5 for < 10 .
K is the reduction factor due to the finite length C ∞ D is the 2D steady drag coefficient C ∞ D for a circular cylinder of infinite length came out to be equal to 1 from Table A6. Therefore, K values can directly be taken as drag coefficients for the flow perpendicular to the pipe. 0.98 1.00 = is the reduction factor due to the finite length is the 2D steady drag coefficient for a circular cylinder of infinite length came out to be equal to 1 from Table A6. Therefore, values can directly be taken as drag coefficients for the flow perpendicular to the pipe. In sway, the flow is parallel to the axis of the pipe. Therefore, was calculated as = .
(A5) was calculated to be greater than 10 in sway for both 30.05 and 26.6 cm pipes. Therefore, the data from Table A7 were used to work out the drag coefficients. In sway, the flow is parallel to the axis of the pipe. Therefore, R e was calculated as R e was calculated to be greater than 10 3 in sway for both 30.05 and 26.6 cm pipes. Therefore, the data from Table A7 were used to work out the drag coefficients. Circular cylinder. Axis parallel to flow = is the reduction factor due to the finite length is the 2D steady drag coefficient for a circular cylinder of infinite length came out to be equal to 1 from Table A6. Therefore, values can directly be taken as drag coefficients for the flow perpendicular to the pipe. In sway, the flow is parallel to the axis of the pipe. Therefore, was calculated as = .
(A5) was calculated to be greater than 10 in sway for both 30.05 and 26.6 cm pipes. Therefore, the data from Table A7 were used to work out the drag coefficients.  Table A5 in sway and heave for R e < 10 5 due to the same outer diameter of the pipe. Moreover, R e was worked out to be greater than 10 3 in the surge; therefore, Table A7 was used. 13.55 cm pipes: Similarly, the drag coefficients were calculated from Table A5 in surge and sway for R e < 10 5 , whereas from Table A7 in heave for R e > 10 3 .
Thruster motors: The axial and vertical thruster motors on Seaperch UUV have a length and diameter of 2.10 cm. Therefore, R e remains the same in surge and sway for both the axial and vertical thrusters, and was worked out to be greater than 10 3 and less than 10 5 at the velocity of 0.167 m/sec. Therefore, Tables A5 and A7 were, respectively, used to work out the drag coefficients. On the other hand, R e was worked out at the velocity of 0.144 m/sec in heave and also came out to be between 10 3 and 10 5 . Therefore , Tables A5  and A7 were used, respectively, to work out the drag coefficients.

Appendix B.2. Drag Terms for Manipulator and Payload Plates
For the manipulator and payload plates, form drag is dominant in heave, whereas frictional drag is dominant in surge and sway. For both the manipulator and payload plates, the R e in heave was worked out to be greater than 10 3 . Therefore, the data from Table A8 were used.
In surge and sway, the friction drag terms were calculated from: where C f is the frictional coefficient and A s is the surface area of the plate. At the velocity of 0.167 m/sec in surge and sway, the R e was worked out to be lower than 10 5 over the flat plate of the manipulator and payload. Therefore, the following equation was used to work out the friction coefficient: Rectangular Plate Normal to Flow Direction locity of 0.144 m/sec in heave and also came out to be between 10 and 10 . Therefore ,  Tables A7 and A5 were used, respectively, to work out the drag coefficients.

Appendix B.2. Drag Terms for Manipulator and Payload Plates
For the manipulator and payload plates, form drag is dominant in heave, whereas frictional drag is dominant in surge and sway. For both the manipulator and payload plates, the in heave was worked out to be greater than 10 . Therefore, the data from Table A8 were used. In surge and sway, the friction drag terms were calculated from: where is the frictional coefficient and is the surface area of the plate. At the velocity of 0.167 m/sec in surge and sway, the was worked out to be lower than 10 over the flat plate of the manipulator and payload. Therefore, the following equation was used to work out the friction coefficient: