Decoupled Hydrodynamic Models and Their Outdoor Identiﬁcation for an Unmanned Inland Cargo Vessel with Embedded Fully Rotatable Thrusters

: Expanding the automation level of the freshly introduced ﬂeet of self-propelled Watertruck + barges, which house fully-rotatable embedded thrusters, might increase their ability to compete with their less sustainable but dominating road-based alternatives. Hydrodynamic motion models, which reveal the manoeuvring capabilities of these barges, can serve as inputs for many pieces of this automation puzzle. No identiﬁed motion models or hydrodynamic data seem to be publicly available for the hull design and the novel actuation system conﬁguration of these barges. Therefore, this study offers: (i) decoupled motion model structures for these barges for surge, sway, and yaw, with a focus on the thruster and damping models; (ii) two identiﬁcation procedures to determine these motion models; (iii) all the experimental data, generated outdoors with a scale model barge to identify (i) based on (ii). In addition, the identiﬁed surge models were compared with both computational and empirical data. These comparisons offer more physical insights into the identiﬁed model structures and can aid in the model selection for which the desired complexity and accuracy evidently depend on their envisaged application. Finally, this methodology need not be limited to the vessel and actuation types utilised by us.


Introduction
The European Watertruck + project [1] aims to revivify the inland waterway transport sector by introducing a novel modular fleet of push vessels and active or passive barges of CEMT hull types I and II [2], which enables the decoupling between sailing and transshipment time. With the purpose of studying the automation potential of this new fleet, Peeters et al. [3] constructed a scale model of a self-propelled Watertruck + barge in conjunction with a complementary inland shore control centre [4] to monitor or control this barge. Several links in this envisaged automation chain can benefit from a hydrodynamic model which describes the motion of this vessel. For example, these models can serve as plant models to verify or tune control algorithms, or as control models to construct model predictive controllers [5], which can also be used for collision avoidance algorithms [6,7]. Furthermore, Kotzé et al. [8] and Peeters et al. [4] suggested that these motion models could be applied to construct dynamic proximity zones surrounding the vessel hull, which could be shared with other ambient vessels or could trigger internal control or monitoring events.
The literature handling these hydrodynamic models appears to hold four main approaches: (i) Transfer functions can produce response models, e.g., relating yaw-rates with rudder changes [9,10]; (ii) Abkowitz [11] introduced a third-order Taylor expansion model to mathematically describe the forces acting on a vessel around its equilibrium state [12]; (iii) artificial neural networks can model the ship manoeuvring motion [13]; and (iv) physically oriented modular models-such as the modular mathematical model [14] and the robot-like vectorial model [15,16]-can offer more physical insights in the vessel dynamics. This study uses the robot-like vectorial model given that it leverages system properties such as symmetry, skew-symmetry, and positiveness of matrices [17], which align with the envisaged future utilisation of the models in this work.
No public hydrodynamic data appear to be available for the CEMT type I-II hulls, and little for inland vessels in general [18]. In addition, the novel self-propelled barges have a non-conventional propulsion system consisting of two fully nested and 360-degree-steerable actuators: a steering grid thruster in the bow and a four-channel thruster in the stern, for which little or no models or data seem to be available too. This propulsion configuration enables a high level of manoeuvrability, crucial for inland vessels which need to navigate independently [19]. Therefore, this study investigates decoupled hydrodynamic motion models for surge, sway, and yaw. This decoupling offers insights in the vessel manoeuvrability on the one hand, and enables a focused investigation of the thruster models on the other hand. More precisely, this study aimed to use the scale model vessel to: (i) Compare different model structures for the damping forces and non-conventional thrusters, with an emphasis on the potential speed-dependent propulsive forces of these thrusters. (ii) Provide two identification methods to determine the coefficients of these or other manoeuvring models: one method uses the instantaneous force balance, and the other one uses the integrated system dynamics by solving the differential equations of motion. (iii) Benchmark the results of the surge models with data from both a computational fluid dynamics (CFD) study [20] and an empirical approach [21]. (iv) Provide all the experimental data of this study in S1 Experimental Data (see Supplementary Materials), together with the details of the utilised data logging and processing approaches.
The achievement of these aims provides more understanding of the manoeuvrability of these vessels and the potential hydrodynamic models to predict such movements in surge, sway, and yaw. In addition, these decoupled models present a stepping stone for subsequent coupled models in the future, and the supplemental experimental data unlock the possibility for other researchers to investigate this novel Watertruck + fleet. This paper continues as follows: Section 2 details the design of the scale model vessel, named the Cogge. Thereupon, Section 3 lists the methods used to achieve the results of Section 4. Finally, Section 5 provides an overarching discussion of these results, and Section 6 concludes this work.

Materials
Section 2.1 handles the geometrical similarities between the Cogge and the real-size Watertruck + barges, and Section 2.2 details their unconventional actuation configurations. Afterwards, Section 2.3 lists the most relevant onboard components used during the outdoor experiments with the Cogge. Figure 1a depicts a push boat pushing four coupled Watertruck + barges and Figure 1b shows the scale model of a Self-Propelled CEMT type-I (SP-CEMT-I) Watertruck + barge, for which the full design details can be found in [3]. Table 1 presents the geometrical similarity-by using a scale factor of λ = 8-between the Cogge and the SP-CEMT-I vessels. The resultant size of the Cogge makes it physically impossible to board this vessel, turning it into an unmanned surface vessel.

360-Degree-Rotatable Thrusters
Figure 2a outlines a bottom view of the hull geometry of the Cogge, demonstrating the water inlets for both thruster systems: on the left hand side for the four-channel stern thruster, and on the right hand side for the steering-grid bow thruster. Both thrusters draw in water via these inlets and exhaust a water stream according to the orientation of their steering mechanism. To illustrate the working principles of these steering mechanisms, Figure 2b,c shows longitudinal cross-sections through the centres of the four-channel and steering-grid thrusters, respectively. Figure 2d,e similarly presents the internal angle-conventions-α i c and α i g -of these thrusters and the working principle of their steering mechanisms, by means of an abstract functional drawing. The steering mechanism of the four channel thruster consists of half a sphere with an opening of approximately 85°which can be rotated 360°. However, the water flow can only exit the vessel through one or two of its four channels, introducing a non-linear mapping between α i c and the resultant orientation of its thrust force. Note that the two longitudinal channels, see Figure 2b, have a small downwards angle, whereas the two-not shown-transversal channels have no such bending. The steering grid thruster can also exhaust its water flow in a full 360 • range. Figure 2c indicates that this systems bends the flow by approximately 180°in the zx-plane before it exits the system. The grid itself has a small angle β to orient the flow underneath the hull.
The bollard pull forces of both thrusters can be found in Figure A1, and more design details together with their main discrepancies with conventional actuation systems in [22]. Since the bow of the Cogge houses the steering-grid and the stern the four-channel thruster, α i b and α i s will denote their internal angles in this study, and n b and n s their propeller speeds. Furthermore, in [3] a flow straightener was placed inside the stern-oriented exhaust channel of the stern thruster to increase the straight line stability of the Cogge. This self made straightener was not installed during the bollard pull tests from Figure A1.
were adapted or reproduced from [22] with permission from MDPI, 2020. Figure 3 summarises the main onboard components and their communication links, for which more details can be found in [3,4]. In summary, a programmable logic controller (PLC) controls the internal thruster angles and rotational speeds of the actuators based on their desired values, which it receives from an industrial computer (I-PC)-or remote controller or Web interface. During the experiments, this I-PC ran the Mission Oriented Operating Suite (MOOS) middleware [23] and received data from the GNSS and IMU sensors. An AsteRx-U Marine Global Navigation Satellite System (GNSS) receiver with two PolaNt* MC antennas [24] was pulse per second (PPS) synchronised with an Ekinox-2E Inertial Measurement Unit (IMU) [25]. The GNSS sensor received real-time kinematic (RTK) corrections via the Flemish Positioning Service, which currently has 45 RTK base stations installed [26]. With these RTK corrections, this GNSS sensor claims to have a horizontal accuracy-depending on the distance to the nearest base station-of 0.6 cm + 0.5 ppm, expressed in twice the distance root mean square. The main GNSS antenna, placed at the stern, expressed the position data and was separated by approximately 4.44 m from the auxiliary antenna, at the bow, which enabled accurate pitch and yaw data [3].

Main Onboard Components
The IMU sensor ran a real-time extended Kalman filter at 200 Hz, where it received corrective steps from the GNSS unit at 5 Hz. A calibration procedure-to find the position and orientation differences between the IMU and GNSS sensors-was conducted using the IMU software; more details and a benchmark of this procedure can be found in [27]. In this study, both the GNSS sensor and the IMU internally logged-PPS synchronised-data at 50 Hz. The GNSS sensor also sent the Coordinated Universal Time (UTC) to the I-PC, where the MOOS middleware logged these data together with the desired and measured actuation system states which were requested at 20 Hz. Accordingly, the UTC stamps provided the connection for the post-processing synchronisation of all the relevant data; see Section 3.2.2.

Methods
Section 3.1 starts with deriving the decoupled equations of motion for surge, sway, and yaw, which will be used throughout this study. Section 3.2 discusses the design of the experiments to generate data in order to identify these models, for which Section 3.3 offers two identification methods. Finally, Section 3.4 provides an overview of the used approaches to compare the results of these identification procedures.

Modelling Assumptions
Equation (1) describes the equations of motion for a vessel in calm water, and thus with frequency-independent coefficients, according to the robot-like vectorial model of [15,17]: (1) Here ν = [u, v, w, p, q, r] denotes the generalised velocity vector expressed in a body-fixed reference frame with origin, O b , and η = [x, y, z, φ, θ, ψ] the generalised position and orientation vector expressed in an inertial reference frame with origin, O i ; Figure 4 illustrates their convention throughout this study. The matrix M RB represents the rigid-body mass matrix, C RB the rigid-body Coriolis and centripetal matrix, M A the added mass matrix, C A the Coriolis and centripetal matrix due to M A , D(ν) the damping matrix, g(η) the gravitational/buoyancy vector, g o the ballast control vector, and finally, τ external accumulates all the external forces acting on the vessel.  The following assumptions were made to model the decoupled hydrodynamic equations of motion for the barge in this study: i.
A local north-east-down reference frame approximates an inertial reference frame to describe the position and orientation of the vessel navigating over short distances; see Figure 4. ii. O b is positioned at the centre of gravity and its axes align with the principal directions of the barge. Furthermore, given the rather prism-shaped geometry of the vessel and its block coefficient of 0.95, the horizontal positions of the centre of gravity, buoyancy, and floatation are assumed to coincide in the middle of the vessel in both length and width. iii. The experimental design of this study aimed at solely studying, and thus exciting, the horizontal motions of the vessel, i.e., neglecting pitch, roll, and heave motions, meaning: η = [x, y, ψ] and ν = [u, v, r] . Furthermore, the remaining motions in the water plane were investigated separately for each degree of freedom. Subsequently, no Coriolis, centripetal, buoyancy, or ballast forces were modelled. iv. Only thruster related external forces were explicitly modelled, i.e., τ external = τ thrust , and their dynamic behaviour was neglected. Do note that for the surge and sway models a bias term was added to the external forces, which was intended to implicitly capture any constant wind forces; see Section 3.3.
The accumulation of these four main assumptions results in the following set of equations which represent the decoupled hydrodynamic motions for the barge in surge, sway, and yaw motion: Accordingly, three hydrodynamic phenomena need to be modelled and identified: the added mass, damping, and actuation forces, which are respectively detailed in Sections 3.1.2-3.1.4

Added Mass Models
In the manoeuvring theory framework [17], the added masses [28] can be modelled by a constant term. Several approaches exist to determine these terms in the absence of well-conditioned towing tank data: (i) numerical approaches (strip theory [29] or panel methods [30]) which are often based on (ii) theoretical values [29,31] for common shapes, (iii) empirical methods [32,33], and (iv) less-conditioned outdoor experiments. Liu et al. [18] successfully implemented the empirical methods of Zhou et al. [32] and Clarke et al. [33] on two scale model reference inland vessels for the Yangtze River. Therefore, this study applied the method of Zhou et al. [32] for the transversal and rotational added masses, and the method of Clarke et al. [33] for the longitudinal added mass, Xu, which equates to: where M denotes the mass of the vessel. Given that our vessel has a higher C B than the reference vessels used to derive (3) [33], the upper limit of (3) was be selected for Xu. Additionally, Zhou et al. [32] suggested that the transversal, Yv, and rotational, Nṙ added masses can be approximated based on the main ship dimensions, listed in Table 1, by the following two relations:

Damping Models
The arising damping forces have many known contributors, e.g., potential damping (due to wave forming), skin friction, vortex shedding, and lifting forces (both due to linear circulation of the water around the hull and due to cross-flow drag). Different modelling assumptions exist to model this complex array of contributors [17,[34][35][36]. Given that: (i.) it is not convenient to independently identify all these contributors, and (ii.) this study investigates the decoupled planar motions of the barge, the complexity of the one degree of freedom damping models was chosen to be cubic, quadratic, or linear, respectively represented by the following three equations:

Propulsion System Models
In the spectrum of the more conventional propulsion systems, tunnel thrusters [37][38][39] and azimuth thrusters [40][41][42] might appear to hold the largest similarities with the non-conventional thrusters of this study [22]. The modelling literature for these tunnel and azimuth thrusters does remain scarce and often uses the open-water propeller characteristics model of [43] as the starting point. Furthermore, the present actuation system might look similar to a water jet system, but it was not designed to perform at the conventional higher speed ranges for these systems [44,45], nor does the overall inlet and outlet design seem to be optimal for such purposes [46]. Do note that some water jet models [47] use the findings of [43] as a modelling starting point, for consistency with the existing hydrodynamic literature. Subsequently, the present work will also start with the first order approximation of the lift forces induced by the propeller blades [43]. Accordingly, the theoretical thrust force, T th , depends quadratically on the shaft speed, n, and bi-linearly on n and the axial inflow speed V A : Introducing the non-dimensional thrust coefficient K T and advance ratio J: where ρ represents the water density and D p the propeller diameter, in conjunction with a linear approximation of K T : gives rise to their connection with (9): where Appending a propeller behind a ship hull changes (9) due to propeller-hull and propeller-wake interactions which can be modelled by the thrust deduction coefficient, t, and the wake fraction, w, respectively: with T the actually available thrust force. Combining (9) and (14) results in: for a conventionally-positioned propeller on a vessel with forward speed. Given (i) the present 360-degrees-rotatable actuators (see Section 2.2), (ii) their angle-dependent bollard thrust forces (see Section 2.2 and Figure A1), and (iii) the manoeuvrability of the vessel, t could be made angle-dependent, and w velocity-dependent, as follows: In [22], the authors demonstrated that T(ν = 0) = [1 − t(α i )]T nn n 2 can serve as a robust theoretical thrust force model for the speed-independent part of (16), by studying the data from Figure A1. Note that t and w denote thrust differences between fully open-water tests (only a propeller) and tests with a propeller attached to a hull. Therefore, the static parts of t(α i ) and w(ν) cannot be explicitly identified in this study, but they will be implicitly incorporated in the thrust coefficients. Nevertheless, their respective angle-and speed-dependency can be investigated.
For the experiments with decoupled vessel motions, i.e., ν = [u, 0, 0] , [0, v, 0] , or [0, 0, r] , three thruster angles suffice, α i ∈ [−90, 0, 90]°(see Section 3.2.1), and assuming vessel and actuator symmetry around the zx-plane narrows these angles down to α i ∈ [0, 90 = −90]°. Furthermore, the water speed at the inlets of the propellers will be assumed to be purely longitudinal or transversal, where for the rotational case the inlet transversal speed at the bow or stern thruster depends on their longitudinal distance L x b/s from the centre of gravity (v inlet = L x b/s r). In addition, given the axisymmetric water inlet at the bottom of the ship hull and the position of the propeller shaft which stands orthogonal to the calm water plane, it will be further assumed that the differences in wake effects for a longitudinal or transversal inlet flow will be negligible. This assumption means, w(ν) can be modelled by a single term, w(ν) = w ν ν, where ν will be assumed to be either pure u or v. Combining these assumptions allows for a further simplification of the thruster models: The thruster model of (17) and (18) will be coined the conventional thrust model in this study. Two other models will be explored too: (19) the bollard thrust model and (20) the direct speed loss thrust model, in order to offer more insights into the potential speed dependency of the thrust forces:

Decoupled Surge, Sway, and Yaw Motion Experiments
In the comparative review of identification papers on ship manoeuvring models in [48], it can be noted that most studies use standard manoeuvres as inputs. In this study, the decoupled motions of surge, sway, and yaw need to be excited independently over a variety of speeds, in order to generate data for the models of Section 3.1.5. Table 3 shows that the actuation system of the Cogge enables three submodes to generate surge motion data: (i) bow thruster only, named "Bow", (ii) stern thruster only, named "Stern", and (iii) both thrusters simultaneously, named "Dual". In addition, the Cogge can sway to the l"Left" or to the "Right", not unlike the crabbing test manoeuvres [49]. Finally, a counter thrust test as suggested by [50] can be performed clockwise or counter-clockwise, named "CW" or "CCW", to induce the yaw motion data.
For each of the aforementioned submodes, three desired time-dependent propeller speed missions, n d b/s (t), were requested from the thrusters named "Ramp," "Stairs," and " Step" in Table 3, and each mission was conducted twice. The first batch of missions, denoted by "m 1 ", were conducted in calm weather conditions, whereas a slight 1-2 Beaufort wind picked up during the second batch, "m 2 ", which can also be seen in Section 4.3.1. For the surge missions, both thrusters had an equal internal angle, α i b/s = 0°, for the sway missions both angles were oriented sideways, α i b/s = 90°or = −90°, and for the yaw missions both thrusters had opposite sideway angles, α i b = −α i s = 90°or = −90°. Figure 5 illustrates the three different desired and measured, n m b/s (t), time-dependent propeller speed profiles. The requested and measured rotational speeds for all missions can be found in the Supplementary Materials S1 Experimental Data. For the sway and yaw motions, the n d b/s were manually determined during the experiments and set to n d b = 2000 rpm and n d s = 500 rpm for all their submodes. This speed ratio gave good visual results-by watching the vessel motion from inside a supply vessel-for both sway and yaw, but it introduces a small discrepancy: since a no-yaw requirement needs equal but counteracting moments of the thrusters, and a no-sway requirement needs equal but counteracting forces of these thrusters. Given the different lever arms of both thrusters, see Figure 2a, it might be assumed that one cannot achieve both situations, i.e., no yaw during sway and vice versa, with the same propeller speed ratio, although this depends on the potential speed-dependency of the thrust forces which was unknown during the experiments. Table 3. Overview of all experiments.

Surge
Sway Yaw

Data Post Processing
Due to the different data logging locations and frequencies (see Section 2.3), the relative time vector of the MOOS database needed to be converted to UTC stamps, and the MOOS data needed to be merged with the other sensor data. Consequently, the following manipulations were done to the MOOS data using the Pandas package [51] in Python version 3 [52]: (i) the time vector was converted to a UTC time vector, (ii) the data were clipped to their relevant mission time range, (iii) the offsets in relative time and position were removed, (iv) the variables were linearly interpolated, given that MOOS allows for asynchronous data logging, and (v) the resultant data were upsampled to 50 Hz. Afterwards, these data were merged with the IMU and GNSS data based on the UTC time vector. The potential time delays between a measurement and its log were assumed negligible and irrelevant compared to the studied vessel dynamics.
The Supplementary Materials S1 Experimental Data holds the results of these manipulations for the variables used in this study in a file per mission, structured according to the taxonomy of Table 3. As a result, each mission file contains: (i) a relative time vector, (ii) the position data of the main GNSS mushroom, which were calculated in the MOOS in real-time using the driver from [53], (iii) the desired and measured α i b/s and n b/s from the MOOS, (iv) the roll, pitch, and yaw angles from the IMU, (v) the yaw-rate from the IMU, (vi) the northern and eastern speeds from the IMU, and (vii) the surge and sway accelerations from the IMU. Note that the IMU publishes these accelerations in a horizontal ship reference frame where the gravity vector has been taken into account. This horizontal frame aligns with the assumptions made in Section 3.1.1. In combination with the decoupled design of the experiments, and thus by neglecting other motions, one can assume these accelerations to be equal to those at the centre of gravity for the surge and sway experiments.
Based on the logged northern and eastern speeds, the course over ground was calculated which allowed the determination of the drift vector. Thereupon, the surge and sway speeds were calculated based on the total speed and drift vectors. In addition, the yaw-accelerations were calculated using the gradient function from Numpy [54]. Finally, all data-to have similar distortions on each variablewere filtered forward and backward-to avoid phase distortions -using a low pass Butterworth filter of order 2 with desired cutoff frequency of 0.3 Hz and Nyquist frequency of 25 Hz, by using the Scipy package [55]. Note that these calculations and filters were not added or performed on the Supplementary Material S1 Experimental Data, in order to provide the original as-measured data sets.

Identification Methods
The offline identification of hydrodynamic vessel models can be roughly divided into frequency and time domain methods [48,56,57]. This study used two different time domain approaches of which the usage need not be restricted to the models suggested in this study. The first method uses the instantaneous force balance of the dynamic equations, given that all inputs and outputs were measured during the experiments; see Section 3.3.1. The second approach integrates the system dynamics by solving its differential equations; see Section 3.3.2. For each model structure x ∈ [a, · · · , i] and motion d ∈ [u, v, r], both methods called the bounded nonlinear least squares function of Virtanen et al. [55], which uses a thrust-region reflective method based on Branch et al. [58] to minimise a cost function, , for k sets of data by altering the values of their parameter vector θ x d : To each translational training mission, i.e., d ∈ [u, v], a bias term β d i was added as an external force, which was intended to capture the static parts of the wind forces. Accordingly, the vector . . , β d k ] was added to the θ x d vectors of the surge and sway cases. In other words, for these two cases, the k training missions share the unknown parameter vector θ x d but each mission has an individual bias term β d i . To facilitate the notation of the cost functions in Sections 3.3.1 and 3.3.2, let the matrices ] capture all the measured thruster states. For x u i and x v i , the outputs u i and v i should be added too, given that some thruster models depend on these speeds. Similarly, x r i needs to be appended with the calculated transversal speeds at the thrusters: v b,i and v s,i . Consequently, each parameter vector θ x d can be similarly divided into two parts which accumulate the parameters relevant for y d i and x d i respectively: θ x y d and θ x x d . Finally, all "left" and "CCW" missions were transformed to "right" and "CW" missions by changing the signs of the measurements.

Force Balance Method (FBM)
Using the knowledge from (2), i.e., the inertial and damping forces should equal the thrust forces, and by adding the individual potential wind bias term for each translation mission, the following three equations show the cost functions for k training sets for the surge, sway, and yaw methods respectively using the FBM approach:

Differential Equation Method (DEM)
In order to integrate the differential equations of motion, only the measured α i b/s and n b/s are needed in combination with the initial vessel speed, which provides the start condition. These integrations were done using the "odeint" function from the SciPy python package [55] which calls the "lsode" solver from the FORTRAN library "odepack" [59]. To further clarify this approach, Equation (28) expresses such a differential equation for M(θ e v ) of (22): Integrating (28) for each time period of 20 ms, where the associated propeller speed and angle measurements are updated accordingly, generates a vector,v, with the predicted sway speeds for that mission. The vectorsû andr can be calculated analogously. The DEM cost functions aim to minimise the difference between these predicted speeds and the measured speeds by changing the values of θ x d and β d , as respectively shown for surge, sway, and yaw speeds: Thus, this DEM method uses an implementation of the prediction error method philosophy [60]. In comparison to methods using the known predictor equation, e.g., [61], which only predicts one step ahead based on the previous measurement, the DEM method can predict the speed profile for the whole time range based on θ x d and β d . Furthermore, the equations of motion need no linearisation, which often is the case (e.g., [62,63]), and no direct accelerations need to be used. The extended Kalman filter of the IMU does use these accelerations to calculate its speeds, so they are implicitly nested in the cost functions.

Model Structure Comparison and Selection
Section 3.4.1 offers the comparison potential for the surge identification results with external data sets from a CFD study and an empirical approach. Afterwards, Section 3.4.2 briefly lists the comparison possibilities solely using the data of this study.

Comparisons with External Data
Peeters et al. [20] applied CFD to investigate the surge resistance of the Cogge. Accordingly, they simulated its hydrodynamic behaviour at five surge speeds, u ∈ [0.2, 0.4, 0.6, 0.8, 1.0] m s , with the OpenFOAM [64] software, using the volume of fluid approach [65] and a k − ω turbulence model [66]. In order to verify their modelling approach, they first conducted the same methodology on a scale model KVLCC2 hull for which towing tank data were present [67]. For this vessel, six speed simulations in the range of u = 0.944 to 1.283 m s resulted in a maximal error of 1.55% compared with the towing tank data. To examine their spatial grid convergence, they used only two grids which resulted in a grid convergence index (GCI) estimation [68] for their coarser grid of around 7%, which depends on the applied safety factor and assumed grid convergence. This GCI suggests an error band around the asymptotic value of the simulations, in which the current solution lies. The calculations for the Cogge were done on a finer grid refinement compared to this coarser grid for the KVLCC2 hull. Hence, perhaps it might be assumed that its GCI would lie in the same order of magnitude. Figure 6a shows the relative pressure field of the water surface around the Cogge, which was longitudinally cut in the middle, had a draft of T CFD = 0.23 m, and had a simulated speed of u = 1.0 m s . For the same configurations, Figure 6b depicts the water speed and streamlines for a longitudinal section of water at the middle of the vessel, i.e., at the top of Figure 6a. Figure 6b demonstrates that, during the initial design phase, the steering-grid thruster at the bow protruded the vessel hull. The physical Cogge has no such protrusion, since its thrusters are fully embedded. In addition, these simulations were done with a total vessel height of 0.35 m, whereas the final hull had a height of approximately 0.43 m, if sensor extensions were neglected. Bearing these assumptions and discrepancies in mind, the CFD calculations resulted in the following second order approximation of the surge damping forces: (a) (b) Figure 6. CFD results with u = 1.0 m s , reproduced from [20] with permission from IEEE, 2020: (a) relative pressure field on the water's surface, and (b) speed and streamlines for a longitudinal section of the water at the bow.
Kristensen et al. [21] extensively reported on an empirical approach based on the ITTC-1957 method [69] to estimate the damping forces of a vessel. Accordingly, the total resistance coefficient C T can be approximated by (33), with R T being the total resistance, S the wetted surface, and the following coefficients: C F the frictional resistance, C A the incremental resistance, C AA air resistance, and C R the residual resistance. Applying the method from Kristensen et al. [21] the following estimations and calculations were made: (i) S according to (34), where ∇ denotes the displacement, (ii) C F by (35), where Re represents the Reynolds number and ν the kinematic viscosity, (iii) C A = 0.0004, given that the hull was not produced according to towing tank standards, and its surface roughness is more in line with larger vessels, (iv) C AA = 0.00007 given that the vessel has plenty of extensions on its deck, and finally (v) for the estimation of C R , [70] provided diagrams based on the length-displacement ratio M d , prismatic coefficient C P where C M denotes the midship coefficient, and the Froude number F r -see (36) for their definitions. This diagram estimation, C R d , should be corrected for several factors, ∑ C R corr : the position of the centre of buoyancy, hull shape and form, bulbous bow shapes, and the B T ratio-see (37). Based on Appendix I from [21], C R d = 0.001, and only a correction for B T was added by the following equation: C R corr, B T = 0.16( B T − 2.5) × 10 −3 . Table 4 lists the additional measured, calculated, or estimated parameters of the geometry of the Cogge. The measured draft of 0.22 m would equal a vessel weight of 623 kg, and vice versa, the measured weight of 590 kg would generate a draft of 0.21 m. The latter weight and draft have been used for the calculations in this study. The geometry of the Cogge, with a C B = 0.95, lies outside the boundaries of some parts of the empirical method of Kristensen et al. [21]. Consequently, it is not straightforward to read, or extrapolate, the right C R d from the diagrams, or to apply all these suggested C R corrections; hence, only the B T correction was added. Considering these limitations, made assumptions, and the neglected additional corrections, the predicted damping forces by this empirical method can thus be reasonably assumed to serve as an underestimation of the actual damping forces. Calculating R T via (33) with u = 1.0 m s , ρ = 997 kg m 3 , ν = 1.15 × 10 −6 m 2 s was used to define (38) which provides an under limit for the surge damping forces of the Cogge. D Emp (u) = 11.82u 2 (38)

Comparisons Based on Experimental Data of this Study
In addition to the two external data comparisons of Section 3.4.1, this section briefly handles three approaches to compare the model structures in this data, based solely on the present experimental data. First, one can use both identification methods of Section 3.3 for the same model structures and see whether both methods independently achieve similar results, which might hint at a good underlying model structure. Second, one can compare the final cost of (24) between model structures for an identification method. Note that one should not compare these costs between identification methods, given their different cost function declarations. Third, for both identification methods, all the data have been split into a verification and a validation set. The former groups all the training missions, whereas the latter gathers the missions which were not used during the identification procedure. Table 5 lists these manually selected validation missions. These missions were chosen to be the missions with the highest initial conditions, i.e., vessel speeds. Accordingly, we aimed to offer a robust test set of data to investigate the model structures and identification methods of this study. This selection should also increase the performance of the FBM, which is more sensitive to these start conditions.

Results
Sections 4.1-4.3 respectively discuss the results for the surge, sway, and yaw motion models. Each section starts with plotting experimental data, continues with the identification and comparison of their model structures, and ends with a suggested model selection. Figure 7 summarises the relevant measured data, which were processed according to Section 3.2.2, from the first-step mission with both thrusters on. Figure 7a reveals that the measured n b and n s overshot before they reached their requested values. This behaviour offers a richer data set, and can also be noted in the measuredu. Figure 7b,c confirms that the sway speeds remained small compared to the surge speeds, and that the heading did not change much during the mission.

Surge Model Structures' Identification and Comparison
The surge motion submode experiments enabled data to be recorded where the thrusters did not cooperate. These independent operation modes, in combination with the dual modes, might offer the most versatile data sets in this study. For this reason, all the surge models, M(θ x u ), have been identified by both the FMB and DEM once using physical boundaries on the thruster parameters, and once with unbounded thruster coefficients, denoted by the superscripts B and U respectively. For both the bounded and unbounded training runs, the initial guesses T 0,b nν 0 = T 0,b nnν 0 = T 0,s nν 0 = T 0,s nnν 0 = 0. The initial guesses for T 0,b nn 0 = 13.7 . These boundaries aim to: (i) take into account the data uncertainty of both thrusters [22], (ii) allow some slack, given that higher rotational speeds were used compared to the towing tank data and the steering grid was now inside a hull, see Figure A1, and (iii) investigate the installed flow straightener, see [3], in the stern thruster which was not installed during the towing tank experiments. Finally, the T 0,b nν , T 0,b nnu , T 0,s nν , T 0,s nnu had an upper limit of zero, and an lower limit which would induce a thrust deduction the size of the available thrust which was taken to be 85% of the initial guess. More precisely, these lower limits were calculated using the following data of two independent step missions: n b = 2700 rpm, u = 1.2 m s for the bow, and n s = 1400 rpm, u = 0.8 m s for the stern. The initial guesses and the boundaries for the added mass and damping models stayed the same throughout all model structures for both the bounded and unbounded methods. The initial guess of Xu = 35.4 Ns 2 m equalled the upper limit of (3). To potentially compensate for the weight measurement accuracy and the high C B of the Cogge, this added mass received the following training boundaries Xu ∈ [25,75]. The present damping parameters received an initial value of seven and the training limits X uuu , X uu , X u ∈ [0, 100]. Evidently, if a damping or thruster parameter was not in θ x u it was not used during its model fit. Table 6 summarises the residuals of the bounded and unbounded FBM and DEM cost functions for all the surge model structures, i.e., (25) and (29) for x ∈ [a, · · · , i]. Table 7 zooms in on the DEM B by listing its identified coefficients. The coefficients for DEM U , FBM B , and FBM U , can be respectively found in the Appendices Tables A1-A3. Based on these identified coefficients, Figure 8 plots the predicted speed profiles, simulated by the differential equations, for all the model structures with a quadratic damping model, i.e., M(θ b u ), M(θ e u ), and M(θ h u ), for three validation and three verification cases. The bias term for each mission was determined by retraining (25) or (29) for only one mission using the identified coefficients and only allowing β u I to vary. For both the validation and the verification missions, the FBM B seems to better capture the data compared to its unbounded alternative. For the bounded and unbounded DEM, there seems to be no significant difference, which can also be noted by the final costs of Table 6. In general, if we judge the model structures on their ability to capture the shape of the data sets, the speed-independent thruster model, M(θ e u ), looks to perform the worst, which Table 6 might be confirming, given that this model structure consistently has a higher final cost compared to M(θ b u ) and M(θ h u ) for all the listed identification approaches.

Surge Model Selection
The identification results of Section 4.1.2 can also be compared with the external data from Section 3.4.1. Accordingly, Figure 9a-d compares the identified cubic damping models for the three different thruster models, i.e., M(θ a u ), M(θ d u ), and M(θ g u ), with the curves defined by (32) and (38). The damping models for the speed-independent thruster model seem to be consistently too large for both the bounded and unbounded FBMs and DEMs. Figure 9e-h similarly contrasts the identified bollard thrust forces of the bow thruster together with its experimental data. Note that this comparison was not made for the stern thruster due to its installed flow straightener which would make it an unfair benchmark. Evidently, the bollard thrust forces for the bounded FBM and DEM stay close to the experimental data, however, both unbounded methods do stay within a close range too, which might hint at a good model structure and/or identification methodology. Finally, given that the conventional thruster model seems to perform well according to Table 6 and Figures 8 and 9, in both DEM B and DEM U , Figure 10 illustrates the three different damping models for this thruster model. In particular, Figure 10a-c shows the simulated the model responses for the three validation cases, Figure 10d compares the damping forces with the curves from (32) (38); (e) predicted bollard thrust bow compared with towing tank data from Peeters et al. [22], and (f) predicted bollard stern thrust with installed flow straightener; hence the towing tank data was not plotted. Figure 11 outlines the relevant captured data, which were processed according to Section 3.2.2, from the first ramp mission where the vessel swayed to the right. This figure displays that both propellers have a lower limit on their rotational speed, hence the straight parts in the ramp profiles of the propeller speeds. Figure 11c reveals that although the heading oscillated a bit during the mission the vessel stayed rather straight. Finally, Figure 11a shows a rather small surge speed compared to the sway speed of Figure 11b.

Sway Model Structures' Identification and Comparison
Considering the good DEM B performance in Section 4.1, these results for the sway model identifications will be exclusively discussed here. The determination of the initial parameter estimations and their training boundaries happened similarly to the surge model approach. Thus, according to (4), the added mass Yv = 383 Ns 2 m , with the limits Yv ∈ [373, 583]. The damping terms started at 500 and had the boundaries of Y vvv , Y vv , Y v ∈ [0, 1.0 × 10 4 ]. Taking into consideration that (i) the surge experiments had the most versatile data sets, (ii) Section 3.1.4 assumed the surge or sway wake effect differences to be negligible, and (iii) both the bounded and unbounded FBMs and DEMs identified speed-dependent thruster coefficients in the same order of magnitude, the start estimations for the coefficients, T 90,b nν 0 , T 90,b nnν 0 , andT 90,s nν 0 , T 90,s nnν 0 , were selected to equal the values of their DEM b identified T(0 • ) equivalents for surge; see Table 7. In order to consider the angle dependency of the thrust deduction (see Figure A1), t, lower and upper limits of 90% and 110% of their initial value were given to these parameters. The initial guesses for T 90,b nn 0 = 14.52 N rpm 2 were calculated based on the towing tank data by assuming a purely quadratic fit to the highest data points closest to the highest propeller speed of the experiments; see Figure A1. Here too, upper and lower limits of 90% and 110% were implemented for both thrusters. These boundaries aimed to: (i) take into account the data uncertainty of both thrusters [22], and (ii) allow some slack, given that higher rotational speeds were used compared to the towing tank data and the steering grid was now inside a hull, see Figure A1. Table 8 states the DEM B identified coefficients for M(θ x v ). Thereupon, Figure 12 predicts the sway speed profiles for three validation and three verification missions for the quadratic damping models, i.e., M(θ b v ), M(θ e v ), and M(θ h v ). The bias terms for each missions were determined similar to the surge speed predictions of Section 4.1.2, and the negative initial sway speeds of two of the validation missions were taken into account in their differential equations. Both Table 8 and Figure 12 indicate minor differences between the different model structures.

Sway Model Selection
With the currently available training data, no model structure seems to outperform its competitors. Figure 13a-c depicts the sway speed predictions for three verification missions for the different damping models using the conventional thruster model. These figures indicate little differences between the linear, quadratic, or cubic damping models. Figure 13d predicts these damping forces, and Figure 13e,f evaluates the bollard pull forces for the bow and stern thruster in comparison with their towing tank data. (f) Thrust force stern.

Yaw Model Structures' Identification and Comparison
Since the yaw experiments used the same thruster angles and propeller speed configurations as the sway experiments, the exact same initial conditions and parameter boundaries were used for these thruster coefficients. The yaw moment of inertia estimation, I z = m(0.2536L) 2 = 878 kgm 2 , followed the approach documented in [18,71]. The initial guess Nṙ = 845 Nms 2 rad , based on (5), received the boundaries Nṙ ∈ [845, 1245]. In addition, all the damping parameters started at 500 and were bounded to N rrr , N rr , N r ∈ [0, 1.0 × 10 04 ]. Table 9 lists the coefficients for M(θ x r ), identified by the DEM B . Subsequently, Figure 15 shows simulations of the yaw rates for three validation and three verification missions for the quadratic damping models, i.e., M(θ b r ), M(θ e r ), and M(θ h r ). Note that the yaw missions have no bias terms. In the same vein as the sway results of Section 4.2.2, both Table 9 and Figure 15 suggest minor differences between the different yaw model structures.

Yaw Model Selection
Comparable with Section 4.2.3, no M(θ x r ) seems to capture the yaw-rate data better than its alternatives. Figure 16a-c compares the conventional thruster model for the three damping models on the three validation missions. These comparisons indicate a better curve fitting for the non-linear damping models. Finally, Figure 16d illustrates the different predicted damping moments, and Figure 16e,f plots the computed bollard thrust forces for the bow and stern thruster.

Discussion
When investigating the results of Section 4, the modelling assumptions and limitations of Section 3.1.1 need to be kept in mind. The provided exemplary mission data sets-see Sections 4.1.1, 4.2.1, and 4.3.1-already aided to judge the adequacy of these assumptions. In addition, Figure 17 sheds more light on the neglected degrees of freedom, i.e., roll, pitch, and heave, for three-step missions of the studied motions of surge, sway, and yaw. This figure reveals that the pitch angle did not vary much for these missions. The roll angle varied more compared to this pitch angle, and even changed almost 3°during the rotational manoeuvre. The plotted heave represents the motions of the IMU which was positioned approximately 1.3 m behind the middle of the vessel [27]. Evidently, these neglected degrees of freedom do impact the hydrodynamics of the vessel. Heave (c) ccw-stairs-m 1 Figure 17. Roll, pitch, and heave for a three-step mission in the three different motion modes: (a) dual-step-m 1 , (b) right-step-m 1 , and (c) ccw-stairs-m 1 .
Bearing the just-mentioned assumptions and limitations in mind, the conventional thruster models, i.e., M(θ a u ), M(θ b u ), and M(θ c u ), seem to outperform their alternative thruster models for the surge motion according to: (i) their cost function results for both the FBM and the DEM-see Table 6; (ii) their predictions of both validation and verification missions-see Figure 8; and (iii) their comparison with the external CFD and empirical methods-see Figures 9 and 10. This combination of evidence might hint that these model structures best capture the occurring hydrodynamics. However, for the yaw and sway motion, there seems to be no real significant performance difference between the different thruster models. Accordingly, Figure 18 attempts to provide more insights into why this might be the case, by plotting the longitudinal or transversal pseudo-advanced speed ratios for surge, sway, and yaw step missions. The vessel speed was taken as advance speed to calculate these ratios; however, the propeller itself does not move in this speed direction, hence the name pseudo-advanced speed ratio. It can be seen that both thrusters achieved higher pseudo-advanced speed ratios in surge than in sway or yaw motion. bow stern (c) ccw-stairs-m 1 Figure 18. Pseudo longitudinal, J u , or transversal, J v , advance speed ratios, at both thrusters for a three-step mission in the three different motion modes: (a) dual-step-m 1 , (b) right-step-m 1 , and (c) ccw-stairs-m 1 .
Whether or not the conventional thruster models of this study capture the conventional advance speed effect of a normal propeller configuration might not be determinable with the current data. It might also be the case that this conventional model structure succeeds in approximating another physical phenomena. Nevertheless, it can be concluded that both speed-dependent thruster models of this study, i.e., (17) and (20), perform better than the speed-independent thruster model, (19), for the surge missions, which did induce higher vessel speeds compared to the sway and yaw cases. For the latter two motion cases, the simpler bollard thrust model (19) seems to suffice to capture the vessel behaviour. Be noted that, running the sway and yaw models without limitations on their speed-dependent thruster parameters, i.e., T 90,b/s nν , T 90,b/s nnν ∈ [0, −∞], resulted in similar identification results where the speed-dependent thrust losses became assigned to only one thruster. This assignment might have been caused due to the fact that both propellers received the same shapes of input signals. Future work could add the no-yaw or no-sway constraints to the cost functions of these sway or yaw modes, which might reduce or avoid this identification effect. Finally, for all motions, the damping forces or moments seem to be capturable with only two parameters: a combination of a linear and a quadratic or cubic term. Consequently, it was demonstrated that a variety of model structures can capture the decoupled hydrodynamic behaviour of the Cogge. Evidently, the selection of a model structure depends on the envisaged requirements of its end user.
Future work could also focus on investigating the coupled equations of motion, where the addition of the roll motion might be beneficial, considering the data of Figure 17. Likewise, the propeller dynamics [39] or the potential thruster-thruster interactions [72] could be investigated. These interactions might be present when the outflow of the steering-grid points towards the stern. Adding a wind sensor and explicitly modelling this disturbance should also offer a better performance than the currently implemented bias terms which were included to capture these effects. It should also be noted that, based on Figure 10, the installed flow straightener seems to have significantly decreased the thrust force of the four-channel thruster. Hence, it might be better to replace this configuration with a rudder in the future. Finally, it should be pointed out that the identification methods used need not be restricted to the suggested model structures or the utilised vessel and its actuators. For example, a similar approach could be applied to the real-size Watertruck + vessels, which would eliminate the scale effects [73] which would occur when one extrapolates the presently identified models. It is noteworthy that the small scale factor of eight should help in decreasing these scale effects [74].

Conclusions
This study detailed the geometry of a scale model of a self-propelled CEMT-I barge which houses a novel actuation system configuration consisting of two fully-rotatable and embedded thrusters. Thereupon, this work derived several model structures for the decoupled equations of motion in surge, sway, and yaw, with a focus on the thruster and damping models. The supplementary materials hold all the outdoor experimental data captured to identify and compare these model structures. Accordingly, this work discussed the onboard sensors which fetched these data sets, and detailed the implemented data post-processing. Furthermore, this work offered two identification methods-one based on the instantaneous force balance, and the other one integrated the differential equations of motion. Both methods were used to identify the derived motion models but they could also be used to determine other manoeuvring models. The identified model structures for the surge motion were compared with an external CFD study and with an empirical approach. Subsequently, the conventional thruster model seems to be a good model structure for the surge motion. Speed-independent thrust models seem to suffice for the yaw and sway motions. Finally, cubic or quadratic models seemed to suit to predict the damping forces or moments for surge, sway, and yaw. . Figures reproduced from [22], with permission from MPDI, 2020. Note that the stern thruster data were fetched with this thruster installed in half a ship hull, whereas the steering-grid was tested as a stand-alone device. Consequently, it may be assumed that the full-internal and external-angle-dependent thrust deductions were captured for the stern thruster, whereas only the, probably dominant, internal angle-dependent deductions were measured for the bow thruster; see [22].