On the theory and application of absolute coordinates-based multibody modelling of the rigid–flexible coupled dynamics of a deep-sea ROV-TMS (tether management system) integrated model

Depósito de Investigación de la Universidad de Sevilla https://idus.us.es/ This is an Accepted Manuscript of an article published by Elsevier in Ocean Engineering, Vol. 258, on august 2022, available at: https://doi.org/10.1016/j.oceaneng.2022.111748 Copyright 2022 Elsevier. En idUS Licencia Creative Commons CC BY-NC-ND On the theory and application of absolute coordinates-based multibody modelling of the rigid-flexible coupled dynamics of a deep-sea ROV-TMS (tether management system) integrated model Thant Zin Htuna,∗,1,Hiroyoshi Suzukiaand Daniel García-Vallejob aDept. of Naval Architecture and Ocean Engineering, Division of Global Architecture, Osaka University, Yamadaoka 2-1, Suita, 565-0871, Osaka Japan bDept. of Mechanical Engineering and Manufacturing, Universidad de Sevilla, Camino de los descubrimientos s/n , Seville, 41092, Sevelle Spain ARTICLE INFO Keywords: ROV Underwater cable dynamics Rigid-flexible multibody dynamics Absolute nodal coordinate formulation Natural coordinates ABSTRACT This paper presents a novel method to investigate the dynamics between the interconnected rigid and very flexible bodies of a deep-sea ROV-TMS with highly coupled dynamics. The aim is to develop an accurate yet effective system model, which is applicable to a wide range of deep-sea working-class ROV applications. The use of different sets of generalized coordinates in rigid-flexible multibody system usually leads to complex forms of constraints Jacobian and quadratic velocity vector when the interconnected rigid and flexible bodies share highly coupled dynamics. Thus, the compatibility of the formulations used in modelling is pivotal to a correct representation of the dynamics of an underwater multibody system with highly coupled dynamics. To this end, a mixed formulation, which is fully-based on absolute coordinates, is presented in modelling the dynamics of underwater tethered vehicle system. The development of system model of a basic ROV-TMS and theoretical implementations are presented in this paper. The applicability of the presented model is assessed through numerical simulations for the realistic scenarios during the deep-sea ROV operations. The numerical studies show that the proposed model is able to capture the nonlinear dynamics of cable and cable drum, and the nonlinear coupling dynamics between the tether cable and the ROV. 1. Introduction The deep regions of ocean space have always drawn substantial attention as the particular frontier explorations because of its abundant marine life and natural resources. Since the olden days, humans have been making continuous efforts to explore the deep oceans and to carry out subsea operations at deep sea either through manned submarines or through unmanned underwater vehicles (UUVs) such as autonomous underwater vehicles (AUVs) and remotely operated underwater vehicles (ROVs). ROVs are used for inspection, maintenance, and repair (IMR) operations while AUVs are mostly used for monitoring or search and rescue (SAR) applications due to autonomous capabilities with their reliance on sensor technology and their high maneuvering speeds. AUVs are normally operated with onboard power supplies whilst ROVs are operated with a physical-link providing the power to the vehicle and transferring data/signals between ROV and the pilot on surface vessel. Due to the distinguishing features of real-time data telemetry and providing the operator with telepresence, the tether-operated underwater systems have been widely applicable in deep sea exploration works. In deep water deployment of the working-class ROVs, they are usually deployed with substation, so-called ROV’s cage, inside of which a TMS is incorporated. The illustration of a working class ROV and its TMS during a deepwater deployment is shown in Fig. 1. During launching and retrieving the ROV, it is engaged to the TMS inside the cage to minimize the hazards and damage to the system through the splash zone. The TMS facilitates the ROV’s maneuverability and increases its operation range without moving the vessel or rig. In addition, it reduces the drag ⋆This document is the results of the research project funded by the JSPS. ∗Corresponding author [email protected] (T.Z. Htun); [email protected] (H. Suzuki); [email protected] (D. García-Vallejo) https://orcid.org/0000-0001-5516-4246 (T.Z. Htun); https://researchmap.jp/read0042800/ (H. Suzuki); https://orcid.org/0000-0002-2319-2688 (D. García-Vallejo) ORCID(s): 0000-0001-5516-4246 (T.Z. Htun) 1 T.Z. Htun, H. Suzuki and D. García-Vallejo: Page 1 of 39 Accepted Manuscript On the theory and application of absolute coordinates-based multibody modelling of a deep-sea ROV-TMS on the tether cable and alleviates the possible vibrations by reeving out the tether cable on demand to minimize the tension loads when operating with ocean currents. Though these systems are usually operated in calm conditions, the loads from marine-environments, and the motion of surface vessel, ROV’s cage and the ROV can still cause tension variations in tether cables. Also, they may tend to be easily excited due to their low weight in water and relatively high axial stiffness (Htun et al. (2020)). Therefore, the tether cable may increase the hydrodynamic drag, and exert inertia forces and nonlinear elastic forces on the ROV. These forces significantly affect the vehicle’s motion, complicating the job of a pilot. Normally, these forces can be counteracted by the human operator or PID-based control systems. For this purpose, the adequacy of proper training and evaluation assessments of requisite skills are mandatory for ROV’s operators. On the other hand, numerical models, which can correctly predict the system response over a wide range of applicability, are crucial and necessary means in advance towards the autonomous control systems. Therefore, accurate models of the ROV-TMS to capture the coupling kinematics and dynamics between the ROV, tether cable and TMS play a pivotal role in the advance towards the development of autonomous control systems and ROV training simulators. Some numerical ROV-tether systems have been reviewed in (Htun et al. (2019)). Vu et al. (2017) modelled an unmanned surface vessel (USV) connected to an ROV with a tether cable based on the multibody dynamics approach. Hong et al. (2020) studied the dynamic motion simulations of an unmanned ocean platform, in which an UUV is linked to an USV via a tether cable. The equations of motion are developed based on the lumped-mass method and Newton- Euler equations. Driscoll et al. (2000a,b) modelled a one-dimensional finite-element lumped-mass model of a deep-sea caged ROV system subject to surface excitations. In (Quan et al. (2015)), a deep-sea ROV system was modelled using a geometrically exact formulation to discretize the umbilical cable. However, their focus was on the dynamic response of the umbilical cable. Eidsvik and Schjølberg (2016) modelled ROV-umbilical system using Euler-Bernoulli beam equations for the umbilical cable. The umbilical cable was assumed to be neutrally buoyant (weightless in water) and the coupled motions between ROV and the cable were established by assuming that the cable is connected to the ROV at its center of gravity with moment-free condition. Deng et al. (2018) used a gradient-deficient ANCF beam to model a tether cable, and the tether was assumed to connect to the ROV with moment-free condition. A winch model was incorporated to the system as a rotating motion to reeve out the cable. Htun et al. (2019) modelled an ROV-tether system by utilizing fully parameterized ANCF beam elements to discretize the tether cable, and the coupling dynamics has been investigated with a specified boundary condition for the actual tether connection of the physical model. In the literature (Park and Kim (2015); Quan et al. (2015); Eidsvik and Schjølberg (2016); Suzuki et al. (2017); Nørve Eidsvik and Schjølberg (2018); Htun et al. (2019)), numerical studies on the dynamics of the underwater tethered vehicles are basically solved following the coupling algorithm shown in Fig. 2. The governing equations of each body are solved separately, and the velocity and position data are numerically integrated and updated iteratively for each body. The obtained data of the driving body are used as inputs to define boundary conditions of other interconnecting driven bodies, usually flexible bodies. This may come with some problems concerning the efficiency and accuracy of the numerical model. Firstly, despite of the strong theoretical foundation, direct numerical implementation of the governing equations without employing any constraint violation elimination technique generally leads to error accumulation with time due to constraint drift (Braun and Goldfarb (2009)). Secondly, the interpretation of the cable’s tension forces at the joint-defined node as reaction forces by the cable acting on the tethered body may lead to inaccurate representation of the coupling dynamics. This is because the said node is subjected to other external forces, inertia forces, internal viscous damping forces, and constrained forces if present. If the reaction forces are not accounted for in a correct manner, the model will poorly predict the interactions between the interconnecting bodies. Thirdly, the application of rotational parameters has numerical issues reported in the literature, relating to kinetic or kinematic singularities in the representation of rigid-body kinematics. Besides, the afore-mentioned works emphasized on the separate studies such as (i) dynamic responses of the underwater umbilical cables and TMS by hydrodynamic forces and motion of the surface vessel, and (ii) investigations on the coupling dynamics between the tether cable and the ROV. Therefore, limitations regarding the different operational scenarios of the ROV’s maneuvers, generality of the models, and the range of applications are therefore of concern. In the literature, the area that has not yet been explored is the coupling dynamics between the TMS, tether cable, and the ROV. This includes the investigation of dynamics of the system during disengaging the ROV from the cage, the operation of the ROV, and retrieving the ROV to the cage. To the best of the authors’ knowledge, this paper is the first work to investigate those problems and present new insights into modelling a deep-sea ROV-TMS in a systematic manner within a framework of rigid-flexible multibody system dynamics. In light of this extent, the T.Z. Htun, H. Suzuki and D. García-Vallejo: Page 2 of 39 Accepted Manuscript On the theory and application of absolute coordinates-based multibody modelling of a deep-sea ROV-TMS disturbances caused by the motion of the surface vessel and the waves-induced excitations to an umbilical cable and the TMS will not be covered by this study. One of the important challenges in modelling a multibody system, consisting of many different rigid and flexible bodies, is how to model the system accurately and efficiently. The choice of formulations to model rigid and flexible bodies in the system is a crucial first-step to accurate and efficient modelling. From the viewpoint of modelling and numerical implementation, the formulations adopted should be compatible to each other and easy to program. This may become a greater concern to many engineers and researchers when the system becomes complicated and contains very flexible bodies that normally undergo large deformation and large rotation. For very flexible bodies, a formulation known as ANCF has been widely used in the multibody dynamics for the last two decades. This method imposes no restrictions on the amount of deformation or rotation and, therefore, it is suited for large deformations and large rotation problems. More importantly, this formulation leads to a constant mass matrix, and it does not use any rotation parameter to describe the body’s orientation and deformation, and position gradient vectors are used instead. In addition, this method has a sturdy framework for the modelling of cables and ropes (Takehara et al. (2011); Htun et al. (2019,2020)), and the cable-pulley/winch applications (Lugrís et al. (2011); Takehara et al. (2016); Bulín et al. (2017); Westin (2018); Suzuki et al. (2020)). Therefore, the ANCF method is used for modelling flexible bodies in the system, which will be a tether cable in this study. Description of motion of the rigid bodies using Cartesian coordinates together with a set of Eulerian angles or rotational parameters leads to a coordinate-dependent mass matrix (García-Vallejo et al. (2003)). Consequently, centrifugal and Coriolis terms appear in the equations of motion. One of the most widely used methods in the dynamic analysis of rigid-flexible multibody system is the Floating Frame of Reference Formulation (FFRF), which employs two sets of generalized coordinates; absolute reference and local elastic coordinates. This formulation can be used to model rigid bodies and flexible bodies subjected to only small strains and small deformations. However, the use of a moving frame of reference, so-called floating frame, causes this method lead to geometric nonlinearities (Sugiyama et al. (2003)). As a result, when a flexible body modelled using ANCF is connected to another rigid body modelled using FFRF in a multibody system, the connectivity conditions are not simply linear, and the algebraic constraint equations are nonlinear even for the cases of simple three-dimensional joints such as revolute and rigid joints. When the types of connections between flexible and rigid bodies that employ two different sets of generalized coordinates, involve relative or constrained rotational motion as well as translational motion, the definitions of constraints and their time derivatives corresponding to velocity and acceleration constraints, which are crucial kinematic terms in describing multibody dynamics of a system, may become highly nonlinear and complex. When modelling a multibody system that consists of both rigid and flexible bodies, a formulation which utilizes fully-Cartesian coordinates, also known as 𝑛𝑎𝑡𝑢𝑟𝑎𝑙 𝑐𝑜𝑜𝑟𝑑𝑖𝑛𝑎𝑡𝑒𝑠, has a robust literature background and it has been used in combination with different formulations for modelling the rigid-flexible multibody systems (García de Jalón et al. (1986); García-Vallejo et al. (2003); Kraus (2006); García de Jalón et al. (2007); García de Jalón (2007); García- Vallejo et al. (2008); García de Jalón and Bayo (2012); Moseid (2017)). It is also easy to formulate different types of joint constraints. This is because the use of the natural coordinates in the formulation of joints strongly facilitates the joining of rigid and flexible bodies, and they also facilitate the interfaces of finite elements in flexible bodies (García de Jalón (2007)). It has been found that natural coordinate formulation (NCF) is very compatible and can be easily formulated in conjunction with the ANCF because both of these formulations use the globally-defined slope vectors to describe angular orientations in the case of rigid-body motions. In addition, only a minimal set of nonlinear constraint equations are required when this formulation and ANCF are used together in a mixed formulation to model, respectively, different rigid and flexible bodies if they are interconnected through kinematic joints in a system (García-Vallejo et al. (2008)). If the points and vectors are located at the joints, being shared by the interconnected bodies, some kinematic joints such as spherical, revolute and rigid joints can be readily defined using a set of linear constraint equations and one can even eliminate the need of evaluating the corresponding Lagrange multipliers (García-Vallejo et al. (2003,2008)). As a result, the generalized coordinates for those kinematic joints can be shared between interconnected bodies, reducing the total number of the system’s coordinates. This can significantly improve the numerical efficiency and reduce the system’s complexity. One of the outstanding features of using this mixed formulation is both of these formulations lead to global consistent mass matrices, resulting velocity-dependent inertia forces such as Coriolis and centrifugal forces are equal to zero. Since the velocity-dependent inertia forces do not appear in the equations of motion of a rigid body, for example, an underwater vehicle in this case, the Coriolis-centrifugal forces of the hydrodynamic added mass acting on the ROV can be distinctively expressed. Through the consistency of the kinetic energy expression, the equivalent T.Z. Htun, H. Suzuki and D. García-Vallejo: Page 3 of 39 Accepted Manuscript On the theory and application of absolute coordinates-based multibody modelling of a deep-sea ROV-TMS expression of the hydrodynamic added mass matrix in the NCF can be easily obtained from the conventional added mass matrix corresponding to the 6-DOFs motion in Newton-Euler equations. To capture the nonlinear behavior of the underwater vehicle, the vehicle’s kinematics as well as the hydrodynamic forces and moments are calculated in the vehicle’s body-fixed frame. Therefore, it does not apply any linearization on the hydrodynamic characteristics of the vehicle. Taking into account the problems mentioned above, this paper presents a novel approach to investigate the dynamics between the interconnected rigid and very flexible bodies of a deep-sea ROV-TMS with highly coupled dynamics. The proposed modelling method utilizes the ANCF to model very flexible bodies; tether cable, and the NCF to model rigid bodies such as the ROV and the cable drum of a basic ROV-TMS. By taking full advantages of the global definitions of position vector gradients, formulation of hydrodynamic forces acting on a flexible tether cable is presented. Based on the NCF framework, modelling rigid-body dynamics of an underwater vehicle, and detailed formulations of hydrodynamic added mass and forces acting on the vehicle are also presented. While the validation of the actual model is not directly addressed in this paper, the validation of the proposed method with the experimental results of a simple tethered body has been presented in (Htun et al. (2021)). In the present study, the applicability of the proposed model is confirmed through numerical studies on the realistic scenarios during the deep-sea ROV operations. By providing system model development and theoretical implementations within the framework of rigidflexible multibody dynamics, this study provides new insights into modelling the dynamics of the ROV-TMS using a mixed formulation fully-based on absolute coordinates. The remaining contents of the article are structured as follows: modelling of tether cable using ANCF beam elements is presented in section 2. Formulation of hydrodynamic forces acting on the tether cable is presented in section 3, followed by the modelling of equations of motion for the ROV in section 4. In section 5, formulation of contact forces between the tether cable and the cable drum is presented. Formulation of kinematic constraints for the tether connections are presented in section 6. Assembly of equations of motions and the integration of the system are presented in section 7. Lastly, numerical simulations and conclusion remarks on the proposed model are presented in section 8and 9, respectively. Figure 1: Schematic image of a working class ROV-TMS system during deep-water deployment. 2. Formulation of the nonlinear internal forces of the tether cable using ANCF The global position vector 𝒓of an arbitrary point on a three-dimensional ANCF beam element can be described using the global shape function 𝑺and a vector of absolute nodal coordinates 𝒆as 𝒓=𝑺(𝑥, 𝑦, 𝑧)𝒆(1) T.Z. Htun, H. Suzuki and D. García-Vallejo: Page 4 of 39 Accepted Manuscript On the theory and application of absolute coordinates-based multibody modelling of a deep-sea ROV-TMS Figure 2: Flowchart of general simulation procedure for ROV-tether system. where 𝑥, 𝑦, 𝑧 are beam’s local coordinates. The absolute velocity and acceleration vectors of the said point can also be written as  𝒓=𝑺 𝒆, 𝒓=𝑺 𝒆(2) where  𝒓and  𝒓are globally-defined velocity and acceleration vectors of a given point. The above kinematic expression allows one to derive the kinetic energy in a very straightforward manner. The isoparametric property of the ANCF elements allows one to integrate the kinetic energy of a given element in a regular domain as follows 𝑇=1 2∫Ω𝑖 𝜌 𝒓𝑇 𝒓dΩ𝑖=1 2∫𝑉𝑖 𝜌 𝒓𝑇 𝒓|𝑱𝑖 0|d𝑉𝑖(3) where 𝜌is the material density, Ω𝑖represents the volume of 𝑖𝑡ℎ element in the current configuration and 𝑉𝑖is the volume of the element in the reference configuration. |𝑱𝑖 0|is the Jacobian of the mapping between the arbitrary initial configuration and straight configurations of the 𝑖𝑡ℎ element at undeformed state. It is worth highlighting that the initial configurations of the element can be arbitrary and need not be straight when using the ANCF beam elements (García- Vallejo et al. (2004,2008); Htun et al. (2020)). If the element is initially straight and its local coordinate system is parallel to the global inertia frame, the matrix 𝑱0is equal to a 3 × 3 identity matrix (Sopanen and Mikkola (2003)). This mapping Jacobian allows one to evaluate the integrand at the arbitrary initial configuration over the reference integral volume at the undeformed straight configuration. The matrix 𝑱0can be defined using the position gradient vectors as 𝑱0=[𝒓,𝑥0𝒓,𝑦0𝒓,𝑧0]=[𝑺,𝑥𝒆0𝑺,𝑦𝒆0𝑺,𝑧𝒆0](4) where 𝑺,𝑥 =𝜕𝑺 𝜕𝑥 and 𝒆0represents a vector of nodal coordinates of an element at the initial reference configuration, which is arbitrary, that is to say, the initial reference configuration can be straight or curved. However, one must ensure that all the stresses in that configuration are zero when the initial stress-free curved elements, which are initially curved in the undeformed reference configuration, are used. From Eq. (3), the element mass matrix 𝑴𝑖 𝑒can be written as 𝑴𝑖 𝑒=∫𝑉𝑖 𝜌𝑺𝑖𝑇𝑺𝑖|𝑱𝑖 0|d𝑉𝑖(5) which is a constant mass matrix defined in the global inertia frame and as a result, Coriolis and centrifugal forces are zero in the equations of motion. However, little advantage can be gained from this feature in the dynamic analysis of T.Z. Htun, H. Suzuki and D. García-Vallejo: Page 5 of 39 Accepted Manuscript On the theory and application of absolute coordinates-based multibody modelling of a deep-sea ROV-TMS subsea cables in underwater applications. The total mass matrix is no longer a constant matrix in the global inertia frame because of the added mass matrix, which is dependent on the time and orientation of the element. Using a continuum mechanics approach, the evaluation of the nonlinear internal forces and their Jacobians can be implemented in a straightforward manner. The symmetric Green-Lagrange strain tensor 𝝐𝑚is calculated as 𝝐𝑚=1 2(𝑱−𝑇 0𝑱𝑇𝑱𝑱−1 0−𝑰3),𝑱=(𝑺,𝑥𝒆 𝑺,𝑦𝒆 𝑺,𝑧𝒆)(6) where 𝑱is the Jacobian of position gradient vectors, superscript (⋅)−𝑇means the transpose of the inverse of a matrix, and 𝑰3is a 3 × 3 identity matrix. Following the approach presented in (Htun et al. (2020)), a generalized force vector of the viscoelastic forces of an isotropic tether element is calculated as 𝑸𝑣𝑒 = − ∫𝑉𝐼 ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ (𝜆+ 2𝜇)𝜅𝑛 3 ∑ 𝑗=1 (𝜕𝜀𝑗𝑗 𝜕𝒆)𝑇 𝜀𝑗𝑗 +4𝜇𝜅𝑛 2 ∑ 𝑗=1 3 ∑ 𝑘>𝑗 (𝜕𝜀𝑗𝑘 𝜕𝒆)𝑇 𝜀𝑗𝑘 +𝜆𝜅2 3 ∑ 𝑗=1 3 ∑ 𝑘=1 𝑘≠𝑗(𝜕𝜀𝑘𝑘 𝜕𝒆)𝑇 𝜀𝑗𝑗 +(𝜆∗+ 2𝜇𝛾𝑑)3 ∑ 𝑗=1 (𝜕 𝜀𝑗𝑗 𝜕 𝒆)𝑇 𝜀𝑗𝑗 +4𝜇𝛾𝑑 2 ∑ 𝑗=1 3 ∑ 𝑘>𝑗 (𝜕 𝜀𝑗𝑘 𝜕 𝒆)𝑇 𝜀𝑗𝑘 +𝜆∗ 3 ∑ 𝑗=1 3 ∑ 𝑘=1 𝑘≠𝑗(𝜕 𝜀𝑘𝑘 𝜕 𝒆)𝑇 𝜀𝑗𝑗 ⎫ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎬ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎭ |𝑱0|d𝑉(7) where 𝜆∗=(𝐾𝛾𝑠)−2 3(𝜇𝛾𝑑), 𝜅𝑛=𝜅1, 𝜅2(8) where 𝜆and 𝜇are Lame’s constants, 𝐾is the bulk modulus, and 𝜺and  𝜺are the Green-Lagrange strain vector and its time-derivative, respectively. A detailed derivation of 𝜸𝑑,𝜸𝑠and the definition of 𝜅𝑛can be found in (Htun et al. (2020)). In this study, for the sake of computational load, a fully-parameterized beam element with 24-DOFs (Yakoub and Shabana (2001)) is used with a null Poisson’s ratio to avoid locking problems, and the internal damping model is incorporated in the same manner as in (Htun et al. (2020)). The averaged effective elastic properties of the cable are obtained from the cantilever bending test (Suzuki et al. (2018)). 3. Formulation of fluid forces acting on the cable A generalized force vector for the weight of the submerged tether cable element is calculated as 𝑸𝑖 ℎ𝑠 =∫𝑺∗𝑇𝒇𝑖 𝑔d𝜉, 𝒇𝑖 𝑔=[0 0 (𝜌𝑐−𝜌𝑤)𝑔𝐴𝑐𝑙]𝑇 (9) where 𝐴𝑐is the cross-sectional area of the tether cable, 𝑙is the underformed length of the tether element and 𝑔, the acceleration due to gravity. 𝜌𝑐and 𝜌𝑤denote the density of the tether cable and density of water, respectively. The shape function: 𝑺∗=𝑺(𝑥, 0,0) is evaluated along the centerline of the ANCF tether element. Since the above integral expression is constant for the submerged tether cable, it can be evaluated at the preprocessing stage prior to the simulation. T.Z. Htun, H. Suzuki and D. García-Vallejo: Page 6 of 39 Accepted Manuscript On the theory and application of absolute coordinates-based multibody modelling of a deep-sea ROV-TMS 3.1. Nonlinear drag force In marine applications such as underwater cables and mooring line applications, it is customary to calculate the drag force components in beam’s local coordinate system by mapping the inflow velocity vector into the local frame and then these forces are transformed to the global inertia frame for the assembly of the forces. In the ANCF framework, some formulations of the hydrodynamic forces of underwater cables for two-dimensional applications can be found in the literature (Takehara et al. (2011); Sheng et al. (2020); Westin (2018)). In this study, the hydrodynamic forces in Morison’s equation are directly calculated in the global inertia frame when one uses a nonlinear finite element method such as ANCF to discretize the cable or mooring lines. The intrinsic definitions of ANCF’s gradient slope vectors are fully utilized instead of using a local transformation matrix, which usually uses Euler’s rotation angles or quaternions, and beam’s local frame. Using such kind of transformation matrices unnecessarily increases the computational cost, and therefore it could be computationally ineffecient. Another advantage of using ANCF beam elements, which employ higher order polynomials in the longitudinal direction, is that obtaining the information of the beam’s orientation or beam’s lay angle in relation to the surroundings are not limited at the nodal points only and can be obtained at an arbitrary point belonging to the element. Therefore, nonlinear hydrodynamic forces can be more accurately predicted when the curvature of the cable segment becomes small and the buckling deformations of the cable becomes large; which are the actual phenomena in the case of ROV’s tether cables. Consequently, a considerably small number of elements is required to accurately represent the hydrodynamic forces along the cable. The unit vector tangent to the cable’s centerline can be written as 𝒆𝑥=𝒓,𝑥 ‖𝒓,𝑥‖= 𝑺∗ ,𝑥𝒆 √𝒆𝑇𝑺∗𝑇 ,𝑥 𝑺∗ ,𝑥𝒆 (10) where 𝒆𝑥is a column unit vector, the subscript (, 𝑥)means 𝜕(⋅) 𝜕𝒙, and ‖⋅‖represents a 𝐿2norm. Hence, the tangential and normal components of the relative velocity vector 𝑽𝑅can be calculated as follows; ⎧ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎩ 𝑽𝑅= 𝒓−𝑼=𝑺∗ 𝒆−𝑼 𝑽𝑅,𝑡 =𝒆𝑥𝒆𝑇 𝑥𝑽𝑅= 𝒓,𝑥𝒓𝑇 ,𝑥 ‖𝒓,𝑥‖2𝑽𝑅=𝑫𝑡𝑽𝑅 𝑽𝑅,𝑛 =(𝑰3−𝒆𝑥𝒆𝑇 𝑥)𝑽𝑅=1 ‖𝒓,𝑥‖2(‖𝒓,𝑥‖2𝑰3−𝒓,𝑥𝒓𝑇 ,𝑥)𝑽𝑅=𝑫𝑛𝑽𝑅 (11) where 𝑼is the inflow velocity vector of the surrounding fluid, and the projection matrices 𝑫𝑡and 𝑫𝑛can be written explicitly in terms of the position gradient vectors as follows 𝑫𝑡(𝑥, 𝑡) = 1 ‖𝒓,𝑥‖2⎡⎢⎢⎢⎢⎣ 𝑟2 1,𝑥 𝑟1,𝑥𝑟2,𝑥 𝑟1,𝑥𝑟3,𝑥 𝑟1,𝑥𝑟2,𝑥 𝑟2 2,𝑥 𝑟2,𝑥𝑟3,𝑥 𝑟1,𝑥𝑟3,𝑥 𝑟2,𝑥𝑟3,𝑥 𝑟2 3,𝑥 ⎤⎥⎥⎥⎥⎦ (12) 𝑫𝑛(𝑥, 𝑡) = 1 ‖𝒓,𝑥‖2⎡⎢⎢⎢⎣ 𝑟2 2,𝑥 +𝑟2 3,𝑥 −𝑟1,𝑥𝑟2,𝑥 −𝑟1,𝑥𝑟3,𝑥 −𝑟1,𝑥𝑟2,𝑥 𝑟2 1,𝑥 +𝑟2 3,𝑥 −𝑟2,𝑥𝑟3,𝑥 −𝑟1,𝑥𝑟3,𝑥 −𝑟2,𝑥𝑟3,𝑥 𝑟2 1,𝑥 +𝑟2 2,𝑥 ⎤⎥⎥⎥⎦ (13) It can be seen that the projection matrices in Eqs. (12) and (13) are defined without imposing any assumption on cable’s deformations. Since the transverse gradient vectors associated with the cross-sectional deformation modes are not included, this can be conveniently used with gradient deficient beam elements for two-dimensional or threedimensional applications as well. Therefore, the definitions are uniquely given and applicable even in the presence of shear and torsional deformations in the case of a fully-parameterized beam element. T.Z. Htun, H. Suzuki and D. García-Vallejo: Page 7 of 39 Accepted Manuscript On the theory and application of absolute coordinates-based multibody modelling of a deep-sea ROV-TMS Norton et al. (1981) reported that the flow-independence principle (Hoerner (1965)) is valid only up to incline angles of about 50◦(the angle between the cable’s centerline and the incident flow). Therefore, the empirical loading functions proposed in (Driscoll and Nahon (1996)) are applied to account for the nonlinear breakup of the tangential and normal drag components as follows; 𝒇𝐷=𝒇𝐷,𝑡 +𝒇𝐷,𝑛 =(𝐷𝑓𝑓𝑡) 𝒕+(𝐷𝑓𝑓𝑛) 𝒏(14) where 𝑓𝑡and 𝑓𝑛are loading functions and their expressions are given in the Appendix. The definitions of drag 𝐷𝑓, 𝒕 and  𝒏can be given by 𝐷𝑓= −1 2𝜌𝑤𝐶𝑑(𝑅𝑒𝑛)𝑑𝑐‖𝑽𝑅‖2, 𝒕=𝑽𝑅,𝑡 ‖𝑽𝑅,𝑡‖, 𝒏=𝑽𝑅,𝑛 ‖𝑽𝑅,𝑛‖(15) where 𝑑𝑐is the diameter of the cable, and semiempirical relations of the drag coefficient 𝐶𝑑and the Reynolds number 𝑅𝑒𝑛, which corresponds to the relative velocity normal to the cable’s centerline, can be given by (Choc and Casarell (1971)); 𝐶𝑑=⎧ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎩ 8𝜋 𝑅𝑒𝑛𝑆(1 − 0.87 𝑆2),(0< 𝑅𝑒𝑛≤1) 1.45 + 8.55 𝑅𝑒0.90 𝑛 ,(1< 𝑅𝑒𝑛≤30) 1.1 + 4 √𝑅𝑒𝑛 ,(30 < 𝑅𝑒𝑛≤105) (16) where 𝑆= −0.077215665 + ln (8 𝑅𝑒𝑛), 𝑅𝑒𝑛=𝑑𝑐‖𝑽𝑅,𝑛‖ 𝜈(17) where 𝜈is the kinematic viscosity of the fluid. The definitions of velocity vectors in relation to the lay of the cable are shown in Fig. 3. It is noteworthy that since the global definitions of velocity components are used, 𝑽𝑅,𝑛 automatically accounts for both the normal and binormal components of the relative velocity vector, which normally need to be calculated by defining and using the local reference frame in conventional approach for three dimensional applications, for example, (Buckham et al. (2003); Nørve Eidsvik and Schjølberg (2018); Quan et al. (2015); Hong et al. (2020); Wu et al. (2022)). Finally, the generalized drag force vector 𝑸𝑖 ℎ𝑑 of an 𝑖𝑡ℎ element can be calculated as 𝑸𝑖 ℎ𝑑 =∫𝑙 𝑺∗𝑇𝒇𝑖 𝐷d𝑥=𝑙∫1 0 𝑺∗𝑇𝒇𝑖 𝐷d𝜉(18) The drag forces are evaluated at each Gaussian quadrature point and then mapped onto the nodal points using the shape function matrix, thereby accounting for the local curvature and local lay angle at each point with respect to the incident flow along the length of the tether element. The definitions of tangential and normal components of the relative inflow velocity vector in relation to the cable’s centerline are illustrated in Fig. 4. 3.2. Hydrodynamic inertia force and added mass For a moving body in an oscillatory flow, the non-drag terms such as hydrodynamic mass force and Froude-Krylov force in Morison’s equation can be given by 𝒇𝐼=𝜌𝑤𝐴𝑐𝐶𝑎( 𝑼𝑛− 𝒓𝑛) ⏟⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏟⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏟ hydrodynamic mass force +𝜌𝑤𝐴𝑐 𝑼𝑛 ⏟⏞⏟⏞⏟ Froude-Krylov force (19) where 𝐶𝑎is the added mass coefficient, 𝐴𝑐is the cable’s cross-sectional area, and  𝑼represents the acceleration vector of the incident flow. Since the amount of fluid displaced along the tangential direction is negligibly small, only the normal component is considered. The expressions in Eq. (19) can be rewritten as follows for further developments: 𝒇𝐼= −𝜌𝑤𝐴𝑐𝐶𝑎 𝒓𝑛+𝜌𝑤𝐴𝑐𝐶𝑚 𝑼𝑛(20) T.Z. Htun, H. Suzuki and D. García-Vallejo: Page 8 of 39 Accepted Manuscript On the theory and application of absolute coordinates-based multibody modelling of a deep-sea ROV-TMS where the first term on the right side in Eq. (45) is an acceleration-dependent inertial force vector while the second and third terms are Coriolis and centrifugal forces, which are quadratic in velocities. In marine hydrodynamics, these Coriolis and centrifugal forces can be obtained in a six-dimensional vector space by using the same hydrodynamic derivatives in Eq. (36) and are presented in the following section. At this stage, one can augment the added mass matrix to the rigid-body mass matrix of the ROV given in Eq. (35) as 𝑴(𝑅𝑂𝑉 )=𝑴𝑅+𝑴𝑁𝑟 𝐴∈ℝ12×12 (46) where 𝑴(𝑅𝑂𝑉 )is the total inertia matrix of the ROV defined in the global inertia frame. 4.4. Formulation of external forces acting on the ROV The external forces acting on the ROV, in general, include propulsive forces, hydrostatic (or) restoring forces by dry weight and buoyancy of the ROV, and hydrodynamic forces. These hydrodynamic forces usually include three types of forces; (i) added mass forces caused by the attached mass of the fluid to an accelerating body, (ii) the nonlinear Coriolis and centrifugal forces associated to the added mass and moment of inertia due to the rotation of ROV’s body-fixed frame relative to the inertial reference frame, and (iii) the hydrodynamic damping forces. The added mass forces have been defined previously in terms of the added mass matrix. Hence, the total external force vector of the tether-free ROV, denoted by  𝑭𝑒𝑥𝑡 ∈ℝ6, is expressed in the body-fixed coordinate system as  𝑭𝑒𝑥𝑡 =[ 𝒇𝑒𝑥𝑡  𝑴𝑒𝑥𝑡]= 𝑛𝑝 ∑ 𝑖[ 𝑻𝑖  𝒓𝑖× 𝑻𝑖]+[ 𝒇𝑔 𝟎]+[ 𝒇𝑏  𝒓𝑏× 𝒇𝑏]−𝑪𝐴( 𝒗) 𝒗−𝑫( 𝒗) 𝒗∈ℝ6(47) where 𝑛𝑝denotes the total number of propellers mounted on the vehicle,  𝒇𝑒𝑥𝑡 and  𝑴𝑒𝑥𝑡 represent an external force vector and moment vector defined in the body-fixed coordinate system, respectively. In Eq. (47), 𝒓𝑏represents the local positioning vector of the center of buoyancy while  𝒓𝑖, the local positioning vector of a point located on the line of action of 𝑖𝑡ℎ propeller with its thrust  𝑻𝑖.𝑪𝐴( 𝒗)and 𝑫( 𝒗)represent the Coriolis-centrifugal matrix, and the damping matrix,  𝒇𝑔and  𝒇𝑏represent the gravity and buoyant forces, respectively, defined in the body-fixed coordinate system as follows  𝒇𝑔=𝑨𝑇 𝑅[00−𝑚𝑔]𝑇 , 𝒇𝑏=𝑨𝑇 𝑅[0 0 𝜌𝑤𝑔▽]𝑇∈ℝ3(48) 𝑪𝐴( 𝒗) = ⎡⎢⎢⎢⎢⎢⎢⎢⎢⎣ 0 0 0 0 −𝑐𝑓3𝑐𝑓2 000𝑐𝑓30 −𝑐𝑓1 0 0 0 −𝑐𝑓2𝑐𝑓10 0 −𝑐𝑓3𝑐𝑓20 −𝑐𝑚3𝑐𝑚2 𝑐𝑓30 −𝑐𝑓1𝑐𝑚30 −𝑐𝑚1 −𝑐𝑓2𝑐𝑓10 −𝑐𝑚2𝑐𝑚10 ⎤⎥⎥⎥⎥⎥⎥⎥⎥⎦ ,𝑫( 𝒗) = 𝑫𝑙+𝑫𝑛( 𝒗) ∈ ℝ6×6 (49) where {𝑐𝑚1=𝐾𝑣𝑣+𝐾𝑝𝑝+𝐾𝑟𝑟, 𝑐𝑚2=𝑀𝑢𝑢+𝑀 𝑤𝑤+𝑀𝑞𝑞, 𝑐𝑚3=𝑁𝑣𝑣+𝑁𝑝𝑝+𝑁𝑟𝑟 𝑐𝑓1=𝑋𝑢𝑢+𝑋 𝑤𝑤+𝑋𝑞𝑞, 𝑐𝑓2=𝑌𝑣𝑣+𝑌𝑝𝑝+𝑌𝑟𝑟, 𝑐𝑓3=𝑍𝑢𝑢+𝑍 𝑤𝑤+𝑍𝑞𝑞(50) where ▽is the submerged volume of the ROV, 𝑫𝑙and 𝑫𝑛( 𝒗)represent the linear and quadratic damping coefficient matrices, respectively. The values of the hydrodynamic derivatives for the added mass and damping matrices of CUROV are obtained from the previous studies (Suzuki et al.,2013;Ohbiki,2014). T.Z. Htun, H. Suzuki and D. García-Vallejo: Page 15 of 39 Accepted Manuscript On the theory and application of absolute coordinates-based multibody modelling of a deep-sea ROV-TMS 4.4.1. Propulsive forces The dynamics of the thrusters mounted on the underwater vehicles are highly nonlinear and complex (Yuh (2000)). Yet, an accurate thruster model is required to adequately represent the thrust force to the underwater vehicles for a given propeller’s revolution speed (Blanke et al. (2000)). Nonetheless, the induced thrust force of a single-screwed propeller can be estimated as 𝑇=𝜌𝑤𝐷4 𝑝𝐾𝑇(𝐽)|𝑛|𝑛(51) where 𝑛is the no. of revolution per second of the propeller, 𝐾𝑇is a thrust coefficient, which is expressed as a function of the advance ratio 𝐽. The advance ratio 𝐽is related to the advanced speed at the propeller 𝑉𝑎, propeller rotation 𝑛 and diameter of the propeller 𝐷𝑝as follows 𝐽=𝑉𝑎 𝑛𝐷𝑝 (52) where 𝑛is defined positive for clockwise rotation and negative for counterclockwise rotation. The relationship between 𝐽and 𝐾𝑇of a RHD-3-bladed propeller (shown in Fig. 8) is given in Fig. 9. The characteristic of the propeller is found Figure 8: RHD-3-bladed propeller. Figure 9: Relation between advance ratio and thruster constant obtained from open water test (Takatsu (2015)). to be approximately linear within the range of positive advance ratio. It is important to highlight that when a propeller is working in the range of negative advance ratio, thrust and torque fluctuate irregularly; which suggests that the propeller is working in a non-steady condition such as vortex ring state. Consequently, the stability and controllability of the ROV can be greatly affected. To avoid this problem, the propeller’s thrust forces are calculated in the simulation based on the linear characteristic of the propeller within the range of positive advance ratio. When the propeller has a positive T.Z. Htun, H. Suzuki and D. García-Vallejo: Page 16 of 39 Accepted Manuscript On the theory and application of absolute coordinates-based multibody modelling of a deep-sea ROV-TMS rotation, it will generate a forward propulsive force, i.e., ROV is driven fowards. Conversely, when the propeller has a negative rotation, a backward propulsive force will exert on the ROV to drive it backwards. It is important to point out that the propeller’s forward and backward thrusts are not symmetrical. Thus far, all the required terms in Eq. (47) have been obtained. In the hydrodynamics of underwater vehicles such as ROVs and AUVs, it is always customary to calculate the propulsive forces, hydrodynamic forces and moments caused by the surrounding fluid in the body-fixed coordinate system. In the classical NCF (García de Jalón et al. (1986); García de Jalón and Bayo (2012); García de Jalón et al. (2007); García de Jalón (2007)), each concentrated force applied at a certain point, which needs to be referred to the natural coordinate system, can be represented by means of an equivalent force vector congruent with the natural coordinates. This equivalent force vector can be obtained using a constant transformation matrix by making use of the principle of virtual work. For an arbitrarily applied concentrated force at point 𝑃, denoted by 𝒇𝑃, the virtual work done can be expressed as 𝛿𝑊𝑃=𝛿𝒓𝑇 𝑃𝒇𝑃=𝛿𝒒𝑇 𝑅𝑪𝑇 𝑃𝒇𝑃(53) where 𝛿𝒓𝑃=𝑪𝑃𝛿𝒒𝑅is the virtual displacement defined in the global inertia frame of an arbitrary point 𝑃, and 𝑪𝑃is the constant transformation matrix which refers the said point 𝑃to the natural coordinate system. From Eq. (53), one can obtain the expression of an equivalent generalized external force vector as follows 𝑸𝑃=𝑪𝑇 𝑃𝒇𝑃(54) However, using this expression for each applied force increases the arithmetic operations and matrix multiplications in the computation when the ROV is actuated by many external forces. In this study, all the external forces excluding the tether coupling forces are therefore first calculated in a six-dimensional vector space using the expression in Eq. (47). An equivalent generalized force vector is then calculated for the resultant force and moment obtained with respect to the center of mass of the ROV. In the formulation of external forces, by so doing, any extra computational load that comes from the increase in the vehicle’s DOFs expressed in the natural coordinate system can be avoided. It is also worth highlighting that since the vehicle’s kinematics and Newton-Euler equations are expressed in the ROV’s local frame, this method does not apply any linearization on the hydrodynamic characteristics of the ROV and therefore it can capture the nonlinear dynamics of the ROV in water. 4.5. Obtaining an equivalent generalized force vector of external moments and torques acting on the ROV in natural coordinate system One should note that in the NCF, imposing concentrated moments and torques, which are applied externally to a rigid body, is not as straightforward as in a formulation that uses rotational parameters, yet not complicated. In natural coordinates-based formulation, any applied torque is replaced by an equivalent pair of two equal and opposite forces, which are separated by a displacement vector, and acting on a plane perpendicular to the direction of the applied torque (García de Jalón and Bayo,2012). Different methods to formulate the moments in natural coordinate system can be found in the literature (García de Jalón and Bayo (2012); Liu et al. (2012)). However, the original method presented in (García de Jalón and Bayo (2012)) can lead to a singular representation of the torque/moment about one local axis in the case of the NCF using two basic points and two unit vectors as pointed out in (Liu et al. (2012)). It was also found that the singularity behavior is present in the case of one basic point and three basic vectors when the resultant moment is acting along the projected axis only. This can be demonstrated by using Fig. 10 (left) and the following equations: ⎧ ⎪ ⎨ ⎪ ⎩ 𝝉1=(𝝉⋅𝒆𝑅1)𝒆𝑅1,𝝉⟂1=𝝉−𝝉1 𝒇1=𝝉⟂1× 𝒖⟂1,𝒇⟂1=𝝉1×𝒆𝑅3, 𝒖⟂1=𝒆𝑅1×𝝉⟂1 ‖𝒆𝑅1×𝝉⟂1‖ (55) where 𝝉is a vector representing a concentrated torque/moment defined in the global inertia frame, 𝝉⟂1and  𝒖⟂1both lie on the plane formed by two unit vectors 𝒆𝑅2and 𝒆𝑅3while 𝒇1is acting in a direction parallel to 𝒆𝑅1axis. As can be seen in Eq. (55), when the value of 𝝉⟂1approaches to zero; this can happen at any time during the normal ROV operation, this leads to the singular representation of the resultant torque/moment. This can be circumvented by using the local definitions of the applied torque/moment instead in finding the equivalent pair of forces straightforwardly. T.Z. Htun, H. Suzuki and D. García-Vallejo: Page 17 of 39 Accepted Manuscript On the theory and application of absolute coordinates-based multibody modelling of a deep-sea ROV-TMS Figure 10: Illustration of an arbitrarily applied concentrated torque on a spatial rigid-body described in the natural coordinate system (García de Jalón and Bayo (2012)) (left) and proposed method (right). Once the forces and moments are obtained from Eq. (47), one can calculate the equivalent generalized force vector 𝑸𝜏 as 𝑸𝜏=[𝒇𝑇 𝑒𝑥𝑡 𝒇𝑇 ⟂1𝒇𝑇 ⟂2𝒇𝑇 ⟂3]𝑇∈ℝ12 (56) where {𝒇𝑒𝑥𝑡 =𝑨𝑅 𝒇𝑒𝑥𝑡 𝒇⟂1=(𝑀𝑒𝑅3𝒆𝑅3)×𝒆𝑅1,𝒇⟂2=(𝑀𝑒𝑅1𝒆𝑅1)×𝒆𝑅2,𝒇⟂3=(𝑀𝑒𝑅2𝒆𝑅2)×𝒆𝑅3 (57) where  𝑴𝑒𝑥𝑡 =(𝑀𝑒𝑅1𝑀𝑒𝑅2𝑀𝑒𝑅3)𝑇 can be obtained from Eq. (47). The definitions of pairs of forces are illustrated in Fig. 10 (right). When the information of the applied moment is only available in the global coordinate system, the local definitions can be simply obtained as  𝑴𝑒𝑥𝑡 =𝑨𝑇 𝑅𝝉(58) Note that the procedure mentioned above can be conveniently used for other NCF-based formulations with different numbers of basic points and basic vectors whenever the body orientation matrix is obtainable. 5. Modelling three-dimensional flexible beam to rigid surface contacts (contacts between the cable and the cable drum) In this study, the Hunt-Crossley contact model is used to describe the contact behavior of the tether cable and cable drum. This contact model has been widely used for the cable and pulley contact models in the literature (Lugrís et al. (2011); Bulín et al. (2017); Westin (2018)). Hunt-Crossley contact model represents the contact surface as a nonlinear spring-damper pair, by which the force per unit length acting on the cable is given by 𝑭𝑛=𝐹𝑛 𝒏𝑖, 𝐹𝑛={𝐾𝑛𝛿𝑛(1 + 𝐷 𝛿)for 𝛿≥0(in contact) 0for 𝛿 < 0(out of contact) (59) where  𝒏𝑖is the normal vector to the drum at any given point along the beam length in contact with the drum (see Fig. 11), 𝛿is the relative radial penetration between the contact surfaces, and  𝛿is the relative velocity of radial penetration, respectively, defined as 𝛿=(𝑅𝑤+𝑑𝑐 2)−√𝒓𝑇 𝑂𝑖 (𝑰3−𝒆𝑤2𝒆𝑇 𝑤2)𝒓𝑂𝑖, 𝛿= − ( 𝒓𝑂𝑖)𝑇 𝒏𝑖(60) T.Z. Htun, H. Suzuki and D. García-Vallejo: Page 18 of 39 Accepted Manuscript On the theory and application of absolute coordinates-based multibody modelling of a deep-sea ROV-TMS Figure 11: Illustration of contacts between the tether and the cable drum. where 𝑑𝑐is the diameter of the cable, 𝑅𝑤is the radius of the cable drum, and 𝒓𝑂𝑖 and  𝒏𝑖are defined below; 𝒓𝑂𝑖 =𝒓𝑖−𝒓𝑂, 𝒏𝑖=(𝑰3−𝒆𝑤2𝒆𝑇 𝑤2)𝒓𝑂𝑖 ‖(𝑰3−𝒆𝑤2𝒆𝑇 𝑤2)𝒓𝑂𝑖‖=(𝑰3−𝒆𝑤2𝒆𝑇 𝑤2)𝒓𝑂𝑖 √𝒓𝑇 𝑂𝑖 (𝑰3−𝒆𝑤2𝒆𝑇 𝑤2)𝒓𝑂𝑖 (61) since ‖(𝑰3−𝒆𝑤2𝒆𝑇 𝑤2)𝒓𝑂𝑖‖=√𝒓𝑇 𝑂𝑖 (𝑰3− 2𝒆𝑤2𝒆𝑇 𝑤2+𝒆𝑤2𝒆𝑇 𝑤2𝒆𝑤2𝒆𝑇 𝑤2)𝒓𝑂𝑖 =√𝒓𝑇 𝑂𝑖 (𝑰3−𝒆𝑤2𝒆𝑇 𝑤2)𝒓𝑂𝑖 (62) where the property of the basic unit vector such that 𝒆𝑇 𝑤2𝒆𝑤2= 1 has been used. The positive exponent 𝑛depends on the geometry of the contact bodies, and the normal contact stiffness per unit length 𝐾𝑛depends on the geometry of bodies, shape of contact surfaces, and material properties. In the case of external contacts, Pereira et al. (2015) suggested a value of 𝑛= 1.094 for most values of the geometric and material properties of the cylinders. In addition, according to the literature (Hunt and Crossley (1975); Flores et al. (2008)), a value of 𝑛 in the range of 1≤𝑛≤1.5is adequate to describe the cylinder to cylinder and cylinder to plane contacts. Another challenge is to determine the values of 𝐾𝑛and damping factor 𝐷. As suggested in (Pereira et al.,2015), the normal contact stiffness per unit length can be estimated as 𝐾𝑛=(0.39 + 0.85 △𝑅)𝐸∗,△𝑅=𝑑𝑐+𝑑𝑤 2(63) where 𝑑𝑐and 𝑑𝑤are the diameters of the cable and the cable drum, respectively. As the contact bodies have different material properties, 𝐸∗is taken as a composite modulus given by 𝐸∗=𝐸𝑐𝐸𝑤 𝐸𝑐(1 − 𝜈2 𝑐)+𝐸𝑤(1 − 𝜈2 𝑤)(64) where 𝐸and 𝜈represent the elastic modulus and Poisson ratio, respectively. As can be seen in Eq. (60), 𝛿can be negative during the contact (𝛿≥0)at each contact point. Therefore, the damping factor 𝐷must be chosen in the sense that 𝒇𝑁is always acting in the direction of  𝒏, for all possible values of  𝛿. Otherwise, the normal contact force will pull the cable towards the center of the cable drum at that contact point, which is impractical in the contact phenomena. In contrast to the cable-pulley contacts, where the travel of the cable is limited only along the sheave groove, the surface of the cable drum can be idealized, for simplicity, to be as the circumference of an infinite cylinder. Finally, the generalized contact forces acting on an 𝑖𝑡ℎ element which comes into contact with the cable drum can be obtained as follows 𝑸𝑖 𝑐𝑜𝑛𝑡 =∫𝑙 𝑺∗𝑇𝑭𝑛d𝑥=𝑙∫1 0 𝑺∗𝑇𝑭𝑛d𝜉(65) T.Z. Htun, H. Suzuki and D. García-Vallejo: Page 19 of 39 Accepted Manuscript On the theory and application of absolute coordinates-based multibody modelling of a deep-sea ROV-TMS It should be mentioned that the cable drum is assumed, for this underwater application, to have a frictionless surface and no friction at the bearing. 6. Formulation of kinematic constraints for tether connections The coupled dynamics between the tether and the ROV depends on the type of tether connection to the ROV, which, in turn, can vary based on the installation mechanism. Htun et al. (2021) presented the formulation of different kinematic joints regarding the common types of tether connections to the ROV. It has also been reported that the conventional modelling method used for underwater tetherd systems with highly coupled dynamics leads to the complex forms of constraints Jacobian and quadratic velocity vectors when the ANCF gradient vectors are used together with rotational parameters. However, as presented in this section, modelling this kind of connection using a mixed formulation fully-based on absolute coordinates is straightforward and easy for computer implementation. Figure 12: Tether connection to ROV and ROV’s cage. In the numerical ROV-TMS model, the tether connections to the cable drum and the ROV can be defined by developing different kinematic joints, respectively. The connection of the tether to the CU-ROV model and the winch system inside the ROV’s cage is shown in Fig. 12. In deep-sea working class ROVs, two types of tether connections are common; (i) a rigid connection: the tether is connected via a water-tight tether socket, and (ii) a spherical joint: the tether is connected via a rotating mechanism to mitigate some coupled forces. For the interested readers, different types of tether connections can be found in (Christ and Wernli Sr (2013)). In the latter case, the coupled kinematic terms between the tether cable and the ROV are linear in most commonly used formulations. In other words, they do not need to capture the coupled dynamics due to the change of orientation of the ROV. An interesting case, which is worth investigating and in-depth study should be given, is therefore a fully-coupled connection between the tether and the ROV. 6.1. Tether connection to the cable drum The type of tether connection to the cable drum depends on the installation, and it, in most cases, can be connected via a spherical joint or a rigid-connection. It should be stressed that, in the case of a rigid-connection to the drum, the element’s length should be small enough to obtain a smooth tether profile around the drum, in particular, if the diameter of the drum is small. Due to the fixed end constraints, the cable will act as a cantilever beam, which can raise stiff response of the system in numerical simulations for a shorter cable length. Therefore, the tether connection to the drum is described as a spherical joint for the convenience. The description of tether connection to the cable drum is shown in Fig. 13. T.Z. Htun, H. Suzuki and D. García-Vallejo: Page 20 of 39 Accepted Manuscript On the theory and application of absolute coordinates-based multibody modelling of a deep-sea ROV-TMS Provided no residual stresses and no cross-sectional deformations are initially present at the tether’s top node, the orientation of the said node can be described by using the ANCF gradient vectors as follows 𝑨(0) 𝑐=[𝒓(0) ,𝑥 𝒓(0) ,𝑦 𝒓(0) ,𝑧 ](66) where superscript (0) denotes the 0𝑡ℎ node (top node) of the tether. Note that the use of the transformation matrix defined in Eq. (66) would require the inclusion of quadratic constraints to ensure that vectors 𝒓(0) ,𝑥 ,𝒓(0) ,𝑦 and 𝒓(0) ,𝑧 remain orthogonal and with unit norms at all times. Hence, the constant transformation matrix that relates the orientation of the tether’s top node, described by frame 𝑨(0) 𝑐and cable drum/winch’s orientation matrix 𝑨𝑤can be written as 𝑨𝑐𝑤 =[𝑨(0) 𝑐]𝑇𝑨𝑤,𝑨𝑤=[𝒆𝑤1𝒆𝑤2𝒆𝑤3](67) If the tether is connected to the cable drum via a spherical joint, the kinematic constraint equation for this connection can be given as 𝚽(0) =𝒓(0) −𝒓𝑃=𝒓(0) −𝑪( 𝒘)𝒒𝑤=𝟎{𝚽(0) ∈ℝ3}(68) where 𝑪( 𝒘)is a constant matrix which references the cable connection point 𝑃to the natural coordinate system in terms of local positioning vector  𝒘, and 𝒒𝑤is a vector of natural coordinates of the cable drum. Figure 13: Tether connection to the cable drum. 6.2. Tether connection to the ROV (socket-type connection) Without any loss of generality, different types of ROV-tether connections can be formulated using the absolute nodal coordinates of the ANCF tether elements, and basic points and basic vectors of the NCF belonging to the ROV. In this type of connection, the tether cable is assumed to be rigidly connected at point 𝐶, which, in general, does not coincide with the center of mass of the ROV, and the ROV is coupling all DOFs with the tether cable at the connection point. Then, the algebraic constraint equations can be written as 𝚽(𝑛)=𝚽(𝑅)=𝑹𝑩𝑅𝒆(𝑛)−𝒒𝑅=𝟎{𝚽(𝑛)∈ℝ12}(69) where 𝑩𝑅is a Boolean matrix that extracts the 12 coordinates of the bottom node where the ROV is attached, and the constant transformation matrix 𝑹is defined as 𝑹=⎡⎢⎢⎣ 𝑰3 𝒙𝑇 𝑐⊗𝑰3 𝟎9×3 𝑨𝑐𝑅 ⊗𝑰3⎤⎥⎥⎦ ∈ℝ12×12 (70) where  𝒙𝑐is a local positioning vector which defines the tether connection point in reference to the vehicle’s coordinate system. Since Eq. (69) is linear, it can be used to express 𝒒𝑅as a function of the nodal coordinate vector, 𝒆. By so doing, T.Z. Htun, H. Suzuki and D. García-Vallejo: Page 21 of 39 Accepted Manuscript On the theory and application of absolute coordinates-based multibody modelling of a deep-sea ROV-TMS the vector 𝒒𝑅will not be part of the generalized coordinate vector of the system, 𝒒. Provided no residual stresses and no cross-sectional deformations are initially present at the tether’s bottom node, a constant rotation matrix that relates the orientation of the ROV and the orthogonal basis at the joint-defined point of the tether can be written as 𝑨𝑐𝑅 =[𝑨(𝑛) 𝑐]𝑇𝑨𝑅(71) where the ROV’s orientation matrix 𝑨𝑅has been defined after Eq. (41) and 𝑨(𝑛) 𝑐=[𝒆1𝒆2𝒆3]=[𝒓(𝑛) ,𝑥 𝒓(𝑛) ,𝑦 𝒓(𝑛) ,𝑧 ](72) where 𝒓(𝑛) ,𝑥 ,𝒓(𝑛) ,𝑦 and 𝒓(𝑛) ,𝑧 are restrained by constraint equations to remain orthogonal to each other and to have unit norms. On the other hand, one can also take the advantage of the similar attributes between the ANCF gradient vectors and the basic vectors of the NCF at the joint-defined point for this tether connection. Provided the orthonormality is satisfied between the ANCF gradient vectors, the kinematics and dynamics of the ROV can be seamlessly augmented to the tether’s bottom node. This can be done by describing the kinetic energy of the ROV in terms of the absolute nodal coordinates of the tether cable at the joint-defined point. To do so, the cross-sectional deformation modes of the tether associated to the transverse shear and torsional stresses should be completely eliminated so that the tether connection satisfies the rigid-connection. Then, an orthogonal basis (𝒆1,𝒆2,𝒆3)is defined at the joint-defined point by using three gradient slope vectors of the bottom node of the tether mesh, provided the algebraic relations, which maintains the unit module of the position gradient vectors and the orthogonality between them, are satisfied. Hence, a corresponding vector of the holonomic constraints 𝚽(𝑂)(𝒒, 𝑡)can be defined as follows: 𝚽(𝑂)(𝒒, 𝑡) = ⎡⎢⎢⎢⎢⎢⎢⎢⎢⎣ ‖𝒓,𝑥‖− 1 ‖𝒓,𝑦‖− 1 ‖𝒓,𝑧‖− 1 𝒓𝑇 ,𝑥𝒓,𝑦 𝒓𝑇 ,𝑥𝒓,𝑧 𝒓𝑇 ,𝑦𝒓,𝑧 ⎤⎥⎥⎥⎥⎥⎥⎥⎥⎦ =𝟎,{𝚽(𝑂)∈ℝ6},𝒒=[𝒒𝑇 𝑤𝒆𝑇]𝑇 (73) where 𝒒is a vector of system’s generalized coordinates. It can be seen that the constraint equations are quadratic in the system coordinates and therefore, it is not cumbersome to obtain the first and second time derivatives of the holonomic constraints. Thus, it is computationally cheap to get the corresponding Jacobian matrix and quadratic velocity vector. Hence, the corresponding Jacobian matrix can be written as 𝜕𝚽(𝑂)(𝒒, 𝑡) 𝜕𝒒= ⎡⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣ ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... 𝒓𝑇 𝑥 ‖𝒓,𝑥‖𝒓𝑇 𝑦 ‖𝒓,𝑦‖𝒓𝑇 ,𝑧 ‖𝒓,𝑧‖ 𝒓𝑇 ,𝑦 𝒓𝑇 ,𝑥 𝒓𝑇 ,𝑧 𝒓𝑇 ,𝑥 𝒓𝑇 ,𝑧 𝒓𝑇 ,𝑦 ⎤⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦ (74) where all the other entries are zeros. It is noteworthy that the six-dimensional subspace of the algebraic relations for rigid-body constraints defined in Eq. (25) are no longer required to set for the ROV separately as one can take the advantage of the above constraint equations to satisfy the rigid-body condition of the ROV. Thus, all the kinematic terms, inertia properties and external forces of the ROV will be conveniently transferred to the cable connection point. The acceleration-independent parts of the acceleration constraints, also known as a quadratic velocity vector, 𝜸can be T.Z. Htun, H. Suzuki and D. García-Vallejo: Page 22 of 39 Accepted Manuscript On the theory and application of absolute coordinates-based multibody modelling of a deep-sea ROV-TMS given as 𝜸(𝑂)= − ⎡⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣  𝒓𝑇 ,𝑥  𝒓,𝑥 ‖𝒓,𝑥‖−( 𝒓𝑇 ,𝑥𝒓,𝑥)2 ‖𝒓,𝑥‖3  𝒓𝑇 ,𝑦  𝒓,𝑦 ‖𝒓,𝑦‖−( 𝒓𝑇 ,𝑦𝒓,𝑦)2 ‖𝒓,𝑦‖3  𝒓𝑇 ,𝑧  𝒓,𝑧 ‖𝒓,𝑧‖−( 𝒓𝑇 ,𝑧𝒓,𝑧)2 ‖𝒓,𝑧‖3 2 𝒓𝑇 ,𝑥  𝒓,𝑦 2 𝒓𝑇 ,𝑥  𝒓,𝑧 2 𝒓𝑇 ,𝑦  𝒓,𝑧 ⎤⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦ ∈ℝ6(75) Provided that the above constraints are imposed on the position gradient vectors, the position gradient vectors: 𝒓,𝑥, 𝒓,𝑦 and 𝒓,𝑧 form the orthogonal basis at all times as follows 𝑨𝐶=[𝒓(𝑛) ,𝑥 𝒓(𝑛) ,𝑦 𝒓(𝑛) ,𝑧 ]=[𝒆1𝒆2𝒆3](76) The frame formed by (𝒆𝟏,𝒆𝟐,𝒆𝟑)will, in general, have a different orientation to that of the ROV’s local frame 𝑨𝑅(see Fig. 14). Since the ROV and the tether share a common point of interest at the tether connection point (point Figure 14: Illustration of tether connection to the ROV using tether’s cross-section frame at the connection point. C, Fig. 14), the kinetic energy of the ROV can be expressed in terms of absolute nodal velocities. By so doing, the number of the system’s total coordinates are reduced by 12. It also eliminates the need of setting joint constraints and of evaluating the corresponding Lagrange multipliers as well. Using a linear relation, the origin of the ROV’s local frame, as represented in Fig. 14, can be expressed in the global inertia frame as follows 𝒓𝐺=[𝑰3 𝜉𝐶1𝑰3 𝜉𝐶2𝑰3 𝜉𝐶3𝑰3]⎡⎢⎢⎢⎢⎣ 𝒓𝐶 𝒆1 𝒆2 𝒆3 ⎤⎥⎥⎥⎥⎦ =𝑪( 𝝃)𝒑𝑅(77) where 𝑪( 𝝃)is a constant matrix containing the components of  𝝃𝑪, which is defined in the frame at point C (Fig. 14). 𝒑𝑅is a vector of natural coordinates, which are solely the absolute nodal coordinates of the bottom node of the tether T.Z. Htun, H. Suzuki and D. García-Vallejo: Page 23 of 39 Accepted Manuscript On the theory and application of absolute coordinates-based multibody modelling of a deep-sea ROV-TMS element. 𝒑𝑅=[𝒓𝐶𝒆1𝒆2𝒆3]𝑇∈ℝ12 (78) Using the relation defined in Eq. (71), the following 12 × 12 constant transformation matrix can be obtained as 𝑹𝐶=⎡⎢⎢⎣ 𝑰3 𝝃𝑇 𝑐⊗𝑰3 𝟎9×3 𝑨𝑐𝑅 ⊗𝑰3⎤⎥⎥⎦ ∈ℝ12×12 (79) One can express the following relations by using the above constant transformation matrix 𝑹𝐶as follows 𝒒𝑅=𝑹𝐶𝒑𝑅, 𝒒𝑅=𝑹𝐶 𝒑𝑅(80) where  𝒑𝑅is a vector of natural velocities, containing the absolute nodal velocities of the bottom node of the tether cable. Therefore, the kinetic energy of the ROV can then be expressed in terms of absolute nodal velocities as 𝑇=1 2 𝒑𝑇 𝑅𝑹𝑇 𝐶𝑴(𝑅𝑂𝑉 )𝑹𝐶 𝒑𝑅(81) This leads to 𝑇=1 2 𝒆𝑩𝑇 𝑅𝑹𝑇 𝐶𝑴(𝑅𝑂𝑉 )𝑹𝐶𝑩𝑅 𝒆(82) where  𝒑𝑅=𝑩𝑅 𝒆,𝑩𝑅is a Boolean connectivity matrix, and  𝒆is a vector of absolute nodal velocities of the bottom node of the tether. To this end, the mass matrix and the external force vector 𝑸𝜏(defined in Eq. (56)) of the ROV can be seamlessly augmented to the tether cable through the Boolean matrix and the constant transformation matrix 𝑹𝐶. 7. Assembly of equations of motion and system integration of an ROV-TMS model Assembling the equations of motion for each respective body belonging to the ROV-TMS model yields the complete system of equations as a monolithic system, which fully describes the ROV-TMS model as an integrated system of flexible and rigid bodies. The equations of motion can be written as an index-1 DAE (differential-algebraic equations) system as ⎡⎢⎢⎣ 𝑴𝚽𝑇 𝑞 𝚽𝑞𝟎⎤⎥⎥⎦⎡⎢⎢⎣  𝒒 𝝀⎤⎥⎥⎦ =⎡⎢⎢⎣ 𝑸 𝜸− 2𝛼 𝚽−𝛽2𝚽⎤⎥⎥⎦ (83) where 𝝀is a vector of Lagrange multipliers and 𝜸represents a vector containing the velocity-dependent terms of the kinematic acceleration constraints. 𝛼and 𝛽are the 𝐵𝑎𝑢𝑚𝑔𝑎𝑟𝑡𝑒’s stabilization parameters. 𝑴is the global assembly of the total mass matrix of the system and 𝒒is a vector of generalized coordinates containing different sets of the absolute coordinates. 𝑸is a generalized force vector containing all external and internal forces of the system. The definitions of kinematic constraints of the system are illustrated in Fig. 15. After employing the constraint-elimination technique presented in section 6.2 for the socket-type connection, the expression of each term in Eq. (83) becomes 𝑴=⎡⎢⎢⎣ 𝑴(𝑤)𝟎 𝟎𝑴(𝑡𝑒𝑡ℎ𝑒𝑟)+𝑩𝑇 𝑅𝑹𝑇 𝐶𝑴(𝑅𝑂𝑉 )𝑹𝐶𝑩𝑅⎤⎥⎥⎦ ,𝒒=[𝒒𝑤 𝒆],𝑸=[𝑸(𝑤) 𝑸(𝑡𝑒𝑡ℎ𝑒𝑟)+𝑩𝑅𝑹𝑇 𝐶𝑸(𝑅𝑂𝑉 )](84) where 𝑴(𝑤)and 𝑸(𝑤)represent the global mass matrix and generalized external force vector of the cable drum, respectively, and 𝑸(𝑡𝑒𝑡ℎ𝑒𝑟)=𝑸ℎ𝑠 +𝑸ℎ𝑑 +𝑸ℎ𝑖 +𝑸𝑐𝑜𝑛𝑡𝑎𝑐𝑡 −𝑸𝑣𝑒 (85) T.Z. Htun, H. Suzuki and D. García-Vallejo: Page 24 of 39 Accepted Manuscript On the theory and application of absolute coordinates-based multibody modelling of a deep-sea ROV-TMS Figure 24: ROV-tether coupled motions at certain selected time steps (spiral turning maneuver, 𝑥𝑦-plane). Figure 25: ROV-tether coupled motions at certain selected time values (spiral turning maneuver, 3D view). where the ramping factor, 𝑅𝑓given in (Orcina Ltd. (2000)), is calculated as 𝑅𝑓=𝑓3(6𝑓2− 15𝑓+ 10)(90) T.Z. Htun, H. Suzuki and D. García-Vallejo: Page 31 of 39 Accepted Manuscript On the theory and application of absolute coordinates-based multibody modelling of a deep-sea ROV-TMS Figure 26: Trajectory of ROV (spiral turning maneuver). where 0≤𝑓≤1is the proportional factor, which is an explicit function of time, given by 𝑓(𝑡)=⎧ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎩ 𝑡 𝜏1 1 𝑡−𝜏2 𝜏1 1 if 0< 𝑡 ≤𝜏1 if 𝜏1< 𝑡 ≤𝜏2 if 𝜏2< 𝑡 ≤𝜏3 if 𝑡>𝜏3 (91) The angular displacement and the angular acceleration can be obtained as follows; 𝜃(𝑡)=𝜔𝑐 ⎧ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎩ ∫𝑡 𝑅𝑓d𝑡 (𝑡−𝜏1 2) [(𝜏2−𝜏1 2)+(𝑓−∫𝑡 𝑅𝑓d𝑡)] 𝜏2 for 0< 𝑡 ≤𝜏1 for 𝜏1< 𝑡 ≤𝜏2 for 𝜏2< 𝑡 ≤𝜏3 for 𝑡>𝜏3 (92) 𝛼(𝑡)=𝜔𝑐 𝜏1 ⎧ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎩ 𝑓2(30𝑓2− 60𝑓+ 30) 0 −𝑓2(30𝑓2− 60𝑓+ 30) 0 for 0< 𝑡 ≤𝜏1 for 𝜏1< 𝑡 ≤𝜏2 for 𝜏2< 𝑡 ≤𝜏3 for 𝑡>𝜏3 (93) Once the solution reaches the static equilibrium, the converged results of the absolute nodal coordinates through Newton-Raphson iterations are stored for the later use. Since the tether is initially winding around the cable drum, it is expected to have the initial stresses. The effect of the winding-induced initial stresses is introduced to the tether cable by calculating the strain and stress tensors using Eq. (6). 8.2.2. Launching the ROV Motion simulation of lunching the ROV is performed as the first example in this section to assess the model’s applicability in simulations that closely match the real operating scenario of the ROV. The cable is reeled out by the T.Z. Htun, H. Suzuki and D. García-Vallejo: Page 32 of 39 Accepted Manuscript On the theory and application of absolute coordinates-based multibody modelling of a deep-sea ROV-TMS Figure 27: Launching simulation of the ROV. Table 5 Simulation conditions for launching the ROV Description Symbols Values Simulation time 𝑇𝑠(s) 60 Overall tether length 𝑙(m) 7 Diameter of cable drum 𝑑𝑤(m) 10 Mass of cable drum 𝑚(kg) 10 Propellers’ rotation 𝑛1, 𝑛2, 𝑛3, 𝑛4, 𝑛5(rps) 0 winch rotational velocity: 0.25rps steadily by using Eqs. (89∼ 93). Simulation parameters are given in Table 5. The density of water is set to 999.7kg/m3so that the ROV itself has the slight gravity to give rise the axial tension of the tether cable. The simulation results after every 5s are given in Fig. 27, and the effect of the winding-induced initial stresses can be observed in 𝑦𝑧-plane motion. However, it should be pointed out that since the present study does not consider the cable-to-cable contacts and frictions between them, the tether is winding around the cable drum as if each winding was overlapping on another around the drum when obtaining its initial configuration. 8.2.3. Recovering the ROV In the deep-sea explorations and inspection works, when the designated mission is done, the ROV is recovered to the TMS or ROV cage and then the TMS is retrieved back to the surface ship. To mimic this scenario, motion simulation for the recovering the ROV has been studied as second example. The same cable length and number of elements are used as in the launching and the simulation time is 70s. In the recovering simulation, the cable is assumed to hang vertically along the 𝑧-axis from the drum at the beginning of simulation. By using the same ramping functions mentioned before, the reeving velocity is steadily increased from zero to the saturated value of 0.3rps, and then decreased to zero velocity. The simulation results at every 5s apart are given in Fig. 28. As expected, small movement of the ROV in the 𝑥𝑦-plane can be observed. This can be explained by the externally-excited pendulum motion of the hanging object and the effect of the fluid forces acting on the open-frame ROV. From the numerical experiments, it has been found that the excited pendulum motion phenomenon strongly depends on the span length of the cable and reeving velocity. Therefore, when the rotation of the cable drum is actuated by applying an external torque, the values of the applied torque should be carefully chosen to obtain the optimum reeving velocity. T.Z. Htun, H. Suzuki and D. García-Vallejo: Page 33 of 39 Accepted Manuscript On the theory and application of absolute coordinates-based multibody modelling of a deep-sea ROV-TMS Figure 28: Recovery simulation of the ROV. Table 6 Thrusters’ revolutions (miscellaneous maneuvers with overall tether length L = 15m) Time (s) motion thruster’s revolution (rps) 𝑛1𝑛2𝑛3𝑛4𝑛5 T=0 ∼22 Reeving in (0.3rps) 0 0 0 0 0 T=22 ∼30 backing+ascending -5 -3 -3 8 8 T=30 ∼40 advan.+ascen.+turning 10 -3 8 10 10 T=40 ∼60 advancing+ascending 5 3 3 2 2 8.2.4. Combined motion simulation of ROV In the preceeding sections, the coupled motion simulations of the tether cable and the ROV, and launching and recovering the ROV have been presented. However, in those simulations, the ROV is driven forward and to ascend only, and the cable’s travel on the drum is small. Therefore, a study case is selected to investigate the model’s applicability in describing the miscellaneous motions of the ROV-tether system. The operating conditions of the ROV are given in Table 6. The simulation results for every 2s are given in Fig. 29. The cable’s travel on the drum can be observed in Fig. 29, which is different from the above simulations where the cable’s travel on the drum occurs only in two dimensions. It can also be confirmed that the proposed model can correctly describe the normal contacts between the tether cable and the drum when the cable’s travel on the drum occurs in three dimensions. As in the case of recovery motion, small movements in the 𝑥𝑦-plane are present while reeving the cable in. For 22 ∼ 30s, the ROV is moving backward while ascending, and only small pitch and yaw motions, as expected, can be observed. During 30 ∼ 40s, the ROV is driven with three-combined motions in such a way that the ROV is ascending and turning while moving forward. It can be observed that the model can execute the prescribed motions. During 40 ∼ 60s, the ROV’s back thrusters and vertical thrusters are, respectively, given for the advancing and ascending motion, but the inputs for the propellers’ revolutions are reduced. Therefore, it can be observed that while the ROV is moving forward, the ROV is descending with pitching angles instead of ascending. This is because of the submerged cable weight is acting on the ROV. This effect is strongly influenced by the deployed length of the tether cable. Through performing the numerical studies for different ROV maneuvers, it has been confirmed that the proposed model can be applicable for the deep-sea underwater ROV-tethered systems. However, it must be pointed out that the T.Z. Htun, H. Suzuki and D. García-Vallejo: Page 34 of 39 Accepted Manuscript On the theory and application of absolute coordinates-based multibody modelling of a deep-sea ROV-TMS Figure 29: Combined motion simulation of the ROV (in 𝑥𝑧-plane and 𝑦𝑧 plane). proposed model has limitations regarding the applications in some operation scenarios in which the tether and ROV may collide each other, and in the worst scenario, the tether is wrapping around the ROV. 9. Conclusions In this study, a redundant set of absolute coordinates-based generalized coordinates, which do not have the kinetic and kinematic singularity problems, are employed to describe the attitude of the underwater vehicle with respect to the global inertia frame. A general form of the equations of motion of a spatial underwater vehicle is presented by taking account of the intrinsic rigid-body constraints for the non-minimal sets of generalized coordinates. The definitions of the inertia forces, propulsive forces and other external forces actuating on the underwater vehicle are given in the NCF framework. Definition of the kinematic joint associated with the socket-type tether connection is presented by taking full advantage of the absolute coordinates. The constraint Jacobian and quadratic velocity vectors are explicitly given for the fully-coupled tether connection. Some distinguishing features of using the NCF in the modeling of underwater vehicle dynamics are summarized as follows; 1. A single set of natural coordinates define the position and orientation of the vehicle directly in the global inertia frame. T.Z. Htun, H. Suzuki and D. García-Vallejo: Page 35 of 39 Accepted Manuscript On the theory and application of absolute coordinates-based multibody modelling of a deep-sea ROV-TMS 2. Rigid-body inertia matrix is directly defined in the global inertia frame and it is a constant matrix, and no velocity-dependent inertial forces therefore appear. Since the governing equation is directly solved in the global inertia frame, it is very convenient to use jointly with other global-coordinates based formulations such as ANCF for flexible bodies. 3. One can also observe that the orthogonal transformation matrix can be easily obtained and is easier to implement when using the NCF. In addition, the velocity transformation matrix is simply composed by the unit basic vectors and therefore it is not quadratic in coordinates like in the case of rotational parameters. 4. When two adjoining bodies share a common point of interest, some constraint equations are not required and the associated Lagrange multipliers can be completely eliminated. As a result, the total number of system variables can be reduced. 5. An equivalent generalized force vector corresponding to hydrostatic restoring forces/moments, hydrodynamic drag and inertia forces/moments of the vehicle can be systematically expressed in natural coordinate system. 6. The need to evaluate the constraint Jacobian and time-derivatives of the constraint equations can be completely avoided when the constraint equations are linear. Therefore, this advantage can be consequential in the model implementation and for computational efficiency when the rigid-flexible multibody system contains many interconnected bodies via such kinematic joints. On the other hand, the number of coupled DOFs is proportional to the increase in the number of system coordinates and nonlinearity in the formulation of joint constraints when the rotational parameters are used. To sum up, the coupling dynamics of the ROV-tether system has been investigated. The effectiveness of the proposed model in launching and recovering has also been studied by performing coupled-motion simulations. Through various motion simulations for different operating scenarios of the ROV, the applicability of the proposed ROVTMS model has been confirmed. By providing system model development and theoretical implementaions within the framework of rigid-flexible multibody dynamics, this study provides new insights into modelling the dynamics of deep-sea ROV systems with higher pair DOFs by using the alternative approach, which is fully-based on the absolute Cartesian coordinates and singularity-free. Declaration of Competing Interests The authors declare that there is no conflict of interest regarding the research materials or methods presented in the manuscript. Acknowledgments This research was supported by the JSPS KAKENHI Grant Number JP16K06903 (Grant-in-Aid for Scientific Research (C): (FY2016-FY2018) and JP19H02359 (Grant-in-Aid for Scientific Research (B): (FY2019-FY2021). We appreciate everyone who may concern. References Blanke, M., Lindegaard, K.P., Fossen, T.I., 2000. Dynamic model for thrust generation of marine propellers. IFAC Proceedings Volumes 33, 353–358. Braun, D.J., Goldfarb, M., 2009. Eliminating constraint drift in the numerical simulation of constrained dynamical systems. Computer Methods in Applied Mechanics and Engineering 198, 3151–3160. doi:10.1016/j.cma.2009.05.013. Buckham, B., Nahon, M., Seto, M., Zhao, X., Lambert, C., 2003. Dynamics and control of a towed underwater vehicle system, part i: model development. Ocean Engineering 30, 453 – 470. doi:10.1016/S0029-8018(02)00029-X. Bulín, R., Hajžman, M., Polach, P., 2017. Nonlinear dynamics of a cable–pulley system using the absolute nodal coordinate formulation. Mechanics Research Communications 82, 21–28. doi:10.1016/j.mechrescom.2017.01.001. Choc, Y.I., Casarell, M.J., 1971. Hydrodynamic resistance of towed cables. Journal of Hydronautics 5, 126–131. doi:10.2514/3.62882. Christ, R.D., Wernli Sr, R.L., 2013. The ROV manual: a user guide for remotely operated vehicles. Butterworth-Heinemann. Deng, F., Su, X., Dong, H., Zheng, D., Zhou, J., Xiong, W., 2018. The motion modeling of rov system using absolute nodal coordinate formulation, in: 2018 3rd International Conference on Advanced Robotics and Mechatronics (ICARM), IEEE. pp. 262–267. Driscoll, F., Lueck, R., Nahon, M., 2000a. Development and validation of a lumped-mass dynamics model of a deep-sea rov system. Applied Ocean Research 22, 169–182. doi:10.1016/S0141-1187(00)00002-X. Driscoll, F., Lueck, R., Nahon, M., 2000b. The motion of a deep-sea remotely operated vehicle system: Part 2: Analytical model. Ocean Engineering 27, 57–76. doi:10.1016/S0029-8018(98)00077-8. Driscoll, R., Nahon, M., 1996. Mathematical modeling and simulation of a moored buoy system, in: OCEANS 96 MTS/IEEE Conference Proceedings. The Coastal Ocean - Prospects for the 21st Century, pp. 517–523 vol.1. T.Z. Htun, H. Suzuki and D. García-Vallejo: Page 36 of 39 Accepted Manuscript On the theory and application of absolute coordinates-based multibody modelling of a deep-sea ROV-TMS Eidsvik, O.A., Schjølberg, I., 2016. Time domain modeling of rov umbilical using beam equations. IFAC-PapersOnLine 49, 452–457. doi:10.1016/j.ifacol.2016.10.447. Flores, P., Ambrósio, J., Claro, J.P., Lankarani, H.M., 2008. Kinematics and dynamics of multibody systems with imperfect joints: models and case studies. volume 34. Springer Science & Business Media. doi:10.1007/978-3-540-74361-3. Fossen, T.I., 2011. Handbook of marine craft hydrodynamics and motion control. John Wiley & Sons. doi:10.1002/9781119994138. García-Vallejo, D., Escalona, J., Mayo, J., Dominguez, J., 2003. Describing rigid-flexible multibody systems using absolute coordinates. Nonlinear Dynamics 34, 75–94. doi:10.1023/B:NODY.0000014553.98731.8d. García-Vallejo, D., Mayo, J., Escalona, J., Dominguez, J., 2004. Efficient evaluation of the elastic forces and the jacobian in the absolute nodal coordinate formulation. Nonlinear Dynamics 35, 313–329. doi:10.1023/B:NODY.0000027747.41604.20. García-Vallejo, D., Mayo, J., Escalona, J., Dominguez, J., 2008. Three-dimensional formulation of rigid-flexible multibody systems with flexible beam elements. Multibody System Dynamics 20, 1–28. doi:10.1007/s11044-008-9103-9. Hoerner, S., 1965. Fluid-Dynamic Drag: Theoretical, Experimental and Statistical Information. Hoerner Fluid Dynamics. URL: https: //books.google.co.jp/books?id=cJs0mwEACAAJ. Hong, S.M., Ha, K.N., Kim, J.Y., 2020. Dynamics modeling and motion simulation of usv/uuv with linked underwater cable. Journal of Marine Science and Engineering 8, 318. doi:10.3390/jmse8050318. Htun, T.Z., Suzuki, H., García-Vallejo, D., 2020. Dynamic modeling of a radially multilayered tether cable for a remotely-operated underwater vehicle (rov) based on the absolute nodal coordinate formulation (ancf). Mechanism and Machine Theory 153, 103961. doi:10.1016/j. mechmachtheory.2020.103961. Htun, T.Z., Suzuki, H., García-Vallejo, D., et al., 2021. On singularity-free multibody modelings of rigid-flexible coupled dynamics of an underwater tethered system, in: The 31st International Ocean and Polar Engineering Conference, International Society of Offshore and Polar Engineers. Htun, T.Z., Suzuki, H., Kuwano, A., et al., 2019. Numerical motion analysis of rov coupled with tether applying 24-dofs absolute nodal coordinate formulation, in: The 29th International Ocean and Polar Engineering Conference, International Society of Offshore and Polar Engineers. Hunt, K.H., Crossley, F.R.E., 1975. Coefficient of restitution interpreted as damping in vibroimpact doi:10.1115/1.3423596. García de Jalón, J., 2007. Twenty-five years of natural coordinates. Multibody System Dynamics 18, 15–33. doi:10.1007/s11044-007-9068-0. García de Jalón, J., Bayo, E., 2012. Kinematic and dynamic simulation of multibody systems: the real-time challenge. Springer Science & Business Media. doi:10.1007/978-1-4612-2600-0. García de Jalón, J., Shimizu, N., Gómez, D., 2007. Natural coordinates for teaching multibody systems with matlab, in: International Design Engineering Technical Conferences and Computers and Information in Engineering Conference, pp. 1539–1548. doi:10.1115/DETC2007-35358. García de Jalón, J., Unda, J., Avello, A., 1986. Natural coordinates for the computer analysis of multibody systems. Computer Methods in Applied Mechanics and Engineering 56, 309–327. doi:10.1016/0045-7825(86)90044-7. Kraus, C., 2006. Efficient Object-Oriented Modelling, Simulation and Parameter Estimation for Biomechanical Problems. Ph.D. thesis. Liu, C., Tian, Q., Hu, H., García-Vallejo, D., 2012. Simple formulations of imposing moments and evaluating joint reaction forces for rigid-flexible multibody systems. Nonlinear Dynamics 69, 127–147. doi:10.1007/s11071-011-0251-8. Lugrís, U., Escalona, J.L., Dopico, D., Cuadrado, J., 2011. Efficient and accurate simulation of the rope-sheave interaction in weight-lifting machine, SAGE Publications Sage UK: London, England. pp. 331–343. doi:10.1177/1464419311403224. Moseid, J.A., 2017. Mathematical Modelling of Cable and Pulley Systems. Master’s thesis. NTNU. Norton, D.J., Heideman, J.C., Mallard, W.W., et al., 1981. Wind tunnel tests of inclined circular cylinders, in: Offshore Technology Conference, Offshore Technology Conference. Nørve Eidsvik, O.A., Schjølberg, I., 2018. Finite element cable-model for remotely operated vehicles (rovs) by application of beam theory. Ocean Engineering 163, 322 – 336. doi:10.1016/j.oceaneng.2018.06.012. Ohbiki, M., 2014. numerical study on motion analysis of underwater vehicle with tether cable. Master’s thesis. Osaka University. Orcina Ltd., 2000. Visual orcaflex user manual, version 9.1a, daltongate, ulverston, cumbria, la12 7aj uk. URL: https://www.orcina.com/ webhelp/OrcaFlex/Content/html/Dynamicanalysis,Ramping.htm. Park, J., Kim, N., 2015. Dynamics modeling of a semi-submersible autonomous underwater vehicle with a towfish towed by a cable. International Journal of Naval Architecture and Ocean Engineering 7, 409–425. doi:10.1515/ijnaoe-2015-0029. Pereira, C., Ramalho, A., Ambrosio, J., 2015. An enhanced cylindrical contact force model. Multibody System Dynamics 35, 277–298. doi:10.1007/s11044-015-9463-x. Quan, W.c., Zhang, Z.y., Zhang, A.q., Zhang, Q.f., Tian, Y., 2015. A geometrically exact formulation for three-dimensional numerical simulation of the umbilical cable in a deep-sea rov system. China Ocean Engineering 29, 223–240. doi:10.1007/s13344-015-0016-0. Sheng, F., Zhong, Z., Wang, K.H., 2020. Theory and model implementation for analyzing line structures subject to dynamic motions of large deformation and elongation using the absolute nodal coordinate formulation (ancf) approach. Nonlinear Dynamics 101, 333–359. doi:10.1007/s11071-020-05783-4. Sopanen, J.T., Mikkola, A.M., 2003. Description of elastic forces in absolute nodal coordinate formulation. Nonlinear Dynamics 34, 53–74. doi:10.1023/B:NODY.0000014552.68786.bc. Sugiyama, H., Escalona, J.L., Shabana, A.A., 2003. Formulation of three-dimensional joint constraints using the absolute nodal coordinates. Nonlinear Dynamics 31, 167–195. doi:10.1023/A:1022082826627. Suzuki, H., Sakaguchi, J., Inoue, T., Watanabe, Y., Yoshida, H., 2013. Evaluation of methods to estimate hydrodynamic force coefficients of underwater vehicle based on cfd. IFAC Proceedings Volumes 46, 197–202. doi:10.3182/20130918-4-JP-3022.00026. Suzuki, H., Tomobe, H., Kuwano, A., Htun, T.Z., et al., 2018. Numerical motion analysis of rov applying ancf to tether cable considering its mechanical property, in: The 28th International Ocean and Polar Engineering Conference, International Society of Offshore and Polar Engineers. Suzuki, H., Tomobe, H., Kuwano, A., Takasu, K., Inoue, T., Ishiwata, J., et al., 2017. Numerical analysis of the motion of rov applying anc method to the motion of tether cable, in: The 27th International Ocean and Polar Engineering Conference, International Society of Offshore and Polar Engineers. T.Z. Htun, H. Suzuki and D. García-Vallejo: Page 37 of 39 Accepted Manuscript On the theory and application of absolute coordinates-based multibody modelling of a deep-sea ROV-TMS Suzuki, H., Yamazoe, A., Htun, T.Z., 2020. Numerical modeling of cable-winch system for rov launching and recovering processes based on the finite element approach, in: The 30th International Ocean and Polar Engineering Conference, International Society of Offshore and Polar Engineers. Takatsu, K., 2015. numerical analysis of the motion of ROV applying ANCF method to the tether cable. Master’s thesis. Osaka University. Takehara, S., Kawarada, M., Hase, K., 2016. Dynamic contact between a wire rope and a pulley using absolute nodal coordinate formulation. Machines 4, 4. doi:10.3390/machines4010004. Takehara, S., Terumichi, Y., Sogabe, K., 2011. Motion of a submerged tether subject to large deformations and displacements. Journal of System Design and Dynamics 5, 296–305. doi:10.1299/jsdd.5.296. Tang, S.C., 1999. Modeling and simulation of the autonomous underwater vehicle, Autolycus. Ph.D. thesis. Massachusetts Institute of Technology. Vu, M.T., Choi, H.S., Nhat, T.Q.M., Ji, D.H., Son, H.J., 2017. Study on the dynamic behaviors of an usv with a rov, in: OCEANS 2017-Anchorage, IEEE. pp. 1–7. Westin, C., 2018. Modelling and simulation of marine cables with dynamic winch and sheave contact. Ph.D. thesis. Carleton University. doi:10.22215/etd/2018-13246. Wu, J., Yang, X., Xu, S., Han, X., 2022. Numerical investigation on underwater towed system dynamics using a novel hydrodynamic model. Ocean Engineering , 110632doi:10.1016/j.oceaneng.2022.110632. Yakoub, R.Y., Shabana, A.A., 2001. Three dimensional absolute nodal coordinate formulation for beam elements: implementation and applications. J. Mech. Des. 123, 614–621. doi:10.1115/1.1410099. Yuh, J., 2000. Design and control of autonomous underwater robots: A survey. Autonomous Robots 8, 7–24. doi:10.1023/A:1008984701078. Appendix Loading functions According to (Driscoll and Nahon (1996)), the loading functions are given by ⎧ ⎪ ⎨ ⎪ ⎩ 𝑓𝑡= 0.01(2.008 − 0.3858𝜂+ 1.9159𝜂2− 4.1615𝜂3+ 3.5064𝜂4− 1.1873𝜂5) 𝑓𝑛= 0.5+0.1(sin 𝜂− cos 𝜂)−0.4 cos 2𝜂− 0.11 sin 2𝜂 (94) where the relative angle 𝜂is calculated as 𝜂= cos−1 (𝑽𝑅,𝑡 ⋅𝑽𝑅 ‖𝑽𝑅,𝑡‖‖𝑽𝑅‖),0≤𝜂≤𝜋 2(95) Jacobians of hydrodynamic forces The use of implicit time-stepping algorithms requires the evaluation of Jacobian of the forces. The superscripts are removed here for simplicity. The Jacobains of the drag forces with respect to the nodal coordinate vector 𝒆can be calculated as follows; 𝜕𝑸ℎ𝑑 𝜕𝒆=𝑙∫1 0 𝐷𝑓𝑺∗𝑇[𝑓𝑡⋅𝜕 𝜕𝒆(𝑽𝑅,𝑡 ‖𝑽𝑅,𝑡‖)+𝑓𝑛⋅𝜕 𝜕𝒆(𝑽𝑅,𝑛 ‖𝑽𝑅,𝑛‖)]d𝜉 =𝑙∫1 0[𝐷𝑓𝑺∗𝑇𝑓𝑡[𝑽𝑅,𝑡 𝜕 𝜕𝒆(1 ‖𝑽𝑅,𝑡‖)+1 ‖𝑽𝑅,𝑡‖ 𝜕𝑽𝑅,𝑡 𝜕𝒆] +𝐷𝑓𝑺∗𝑇𝑓𝑛[𝑽𝑅,𝑛 𝜕 𝜕𝒆(1 ‖𝑽𝑅,𝑛‖)+1 ‖𝑽𝑅,𝑛‖ 𝜕𝑽𝑅,𝑛 𝜕𝒆]]d𝜉 (96) where the gradient of the L2 norms and partial derivatives of 𝑽𝑅,𝑡 with respect to the nodal coordinate vector 𝒆can be given respectively as follows; ⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ 𝜕‖𝑽𝑅,𝑡‖−1 𝜕𝒆= − 𝑽𝑇 𝑅,𝑡 ‖𝑽𝑅,𝑡‖3(𝜕𝑽𝑅,𝑡 𝜕𝒆) 𝜕‖𝑽𝑅,𝑛‖−1 𝜕𝒆= − 𝑽𝑇 𝑅,𝑛 ‖𝑽𝑅,𝑛‖3(𝜕𝑽𝑅,𝑛 𝜕𝒆)(97) T.Z. Htun, H. Suzuki and D. García-Vallejo: Page 38 of 39 Accepted Manuscript On the theory and application of absolute coordinates-based multibody modelling of a deep-sea ROV-TMS 𝜕𝑽𝑅,𝑡 𝜕𝒆=𝜕(𝑫𝑡𝑽𝑅) 𝜕𝒆=𝜕 𝜕𝒆(𝑺𝑥𝒆𝒆𝑇𝑺𝑇 𝑥𝑽𝑅 𝒆𝑇𝑺𝑇 𝑥𝑺𝑥𝒆)=𝒆𝑇𝑺𝑇 𝑥𝑽𝑅𝑺𝑥+𝑺𝑥𝒆𝑽 𝑇 𝑅𝑺𝑥 𝒆𝑇𝑺𝑇 𝑥𝑺𝑥𝒆≡𝒓𝑇 𝑥𝑽𝑅𝑺𝑥+𝒓𝑥𝑽𝑇 𝑅𝑺𝑥 ‖𝒓𝑥‖2(98) Recalling the relation in Eq. (11), partial derivatives of 𝑽𝑅,𝑛 with respect to the nodal coordinate vector 𝒆can be given as 𝜕𝑽𝑅,𝑛 𝜕𝒆= −𝜕𝑽𝑅,𝑡 𝜕𝒆(99) The partial derivatives of the drag forces with respect to the nodal velocity vector  𝒆can be calculated as follows; 𝜕𝑸ℎ𝑑 𝜕 𝒆=𝑙Λ∫1 0 𝑺∗𝑇[𝑓𝑡 𝜕 𝜕 𝒆(‖𝑽𝑅‖2𝑽𝑅,𝑡 ‖𝑽𝑅,𝑡‖)+𝑓𝑛 𝜕 𝜕 𝒆(‖𝑽𝑅‖2𝑽𝑅,𝑛 ‖𝑽𝑅,𝑛‖)]d𝜉 =𝑙Λ𝑓𝑡∫1 0[2𝑓𝑡𝑺∗𝑇𝑽𝑅 𝑽𝑇 𝑅,𝑡 ‖𝑽𝑅,𝑡‖𝑺∗+𝑓𝑡𝑺∗𝑇‖𝑽𝑅‖2(𝑫𝑡𝑺∗ ‖𝑽𝑅,𝑡‖−𝑽𝑅,𝑡 𝑽𝑇 𝑅,𝑡 ‖𝑽𝑅,𝑡‖3𝑫𝑡𝑺∗) + 2𝑓𝑛𝑺∗𝑇𝑽𝑅 𝑽𝑇 𝑅,𝑛 ‖𝑽𝑅,𝑛‖𝑺∗+𝑓𝑛𝑺∗𝑇‖𝑽𝑅‖2(𝑫𝑛𝑺∗ ‖𝑽𝑅,𝑛‖−𝑽𝑅,𝑛 𝑽𝑇 𝑅,𝑛 ‖𝑽𝑅,𝑛‖3𝑫𝑛𝑺∗)]d𝜉 (100) where Λ = −1 2𝜌𝑤𝐶𝑑(𝑅𝑒𝑛)𝑑𝑐, and the partial derivatives of the L2 norms with respect to the nodal velocity vector  𝒆 are obtained as follows; ⎧ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎩ 𝜕‖𝑽𝑅‖2 𝜕 𝒆 𝒕=𝜕𝑽𝑇 𝑅𝑽𝑅 𝜕 𝒆 𝒕= 2𝑽𝑅 𝜕𝑽𝑇 𝑅 𝜕 𝒆 𝒕= 2𝑽𝑅 𝒕𝑇𝑺∗ 𝜕‖𝑽𝑅,𝑡‖−1 𝜕 𝒆= − 𝑽𝑇 𝑅,𝑡 ‖𝑽𝑅,𝑡‖3(𝜕𝑽𝑅,𝑡 𝜕 𝒆)= − 𝑽𝑇 𝑅,𝑡 ‖𝑽𝑅,𝑡‖3𝑫𝑡𝑺∗ 𝜕‖𝑽𝑅,𝑛‖−1 𝜕 𝒆= − 𝑽𝑇 𝑅,𝑛 ‖𝑽𝑅,𝑛‖3(𝜕𝑽𝑅,𝑛 𝜕 𝒆)= − 𝑽𝑇 𝑅,𝑛 ‖𝑽𝑅,𝑛‖3𝑫𝑛𝑺∗ (101) It must be pointed out here that according to Eqs. (95) and (94),𝑓𝑡and 𝑓𝑛are functions of the nodal velocity vector  𝒆. However, considering this dependency in the evaluation of Jacobian of the hydrodynamic drag forces leads to a more complicated expression. Therefore, it is assumed that the Jacobian expresion in Eq. (100) is accurate enough. In a similar manner, 𝜕𝑸ℎ𝑖 𝜕𝒆can be calculated in the case of an unsteady flow. Jacobians of normal contact forces The Jacobian of the contact forces with respect to the nodal coordinates and nodal velocities of the element is given by ⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ 𝜕𝑸𝑐𝑜𝑛𝑡 𝜕𝒆=𝑙∫1 0 𝑺∗𝑇𝜕𝑭𝑛 𝜕𝒆d𝜉 𝜕𝑸𝑐𝑜𝑛𝑡 𝜕 𝒆=𝑙∫1 0 𝑺∗𝑇𝜕𝑭𝑛 𝜕 𝒆d𝜉 (102) where ⎧ ⎪ ⎨ ⎪ ⎩ 𝜕𝑭𝑛 𝜕𝒆=𝐾𝑛 𝜕 𝜕𝒆(𝛿𝑛 𝒏+𝛿𝑛𝐷 𝛿 𝒏) 𝜕𝑭𝑛 𝜕 𝒆= − 𝐾𝑛𝛿𝑛𝐷 𝒏𝑖 𝒏𝑇 𝑖𝑺∗ (103) T.Z. Htun, H. Suzuki and D. García-Vallejo: Page 39 of 39 Accepted Manuscript View publication stats