Method to Develop Legs for Underwater Robots: From Multibody Dynamics with Experimental Data to Mechatronic Implementation

Exploration of the seabed may be complex, and different parameters must be considered for a robotic system to achieve tasks in this environment, such as soil characteristics, seabed gait, and hydrodynamic force in this extreme environment. This paper presents a gait simulation of a quadrupedal robot used on a typical terrigenous sediment seabed, considering the mechanical properties of the type of soil, stiffness, and damping and friction coefficients, referenced with the specialized literature and applied in a computational multibody model with many experimental data in a specific underwater environment to avoi hydrodynamic effects. The requirements of the positions and torque in the robot’s active joints are presented in accordance with a 5R mechanism for the leg and the natural pattern shown in the gait of a dog on the ground. These simulation results are helpful for the design of a testbed, with a leg prototype and its respective hardware and software architecture and a subsequent comparison with the real results.


Introduction
Many devices have been modeled and manufactured for underwater exploration based on two main concepts: autonomy and remote operation. Autonomy is achieved with autonomous underwater vehicles (AUVs) [1] for survey tasks, seabed mapping, identification/inspection, rescue, etc., and remote operation with remotely operated vehicles (ROVs) [2] for observation, marine sampling, photography, etc. At present, other concepts have been deisnged for this task, such as legged devices, which offer other capabilities to achieve underwater exploration, including the two mentioned concepts and the idea of bio-inspired mechanisms, which are usually designed in accordance with underwater gait stability using hexapod anatomy.
In [3], a hexapod robot walking underwater considering hydrodynamics effects to improve the stability margin was considered. Picardi [4] reported a bioinspired hexapod robot with two locomotion modes: hopping and walking, where the leg motors reached the position of the next state, depending on which a momentum around a vertical axis could be generated, resulting in a rotation movement. A simulation of the mechanical modeling, seafloor environment setting, and gait planning for a hexapod robot was presented by Liu [5], who verified the correctness of movement using a virtual model. Wang Z. [6] reported a model and simulation of an underwater hexapod take in account the kinematic, dynamic, and hydrodynamics effects using an isolated body method for a virtual model. The principal parameters than influence the underwater gait were presented by Wang G. [7], with a method to compute join torque, motion parameters, hydrodynamic force, and modeling of foot-terrain force.
Many types of methods suitable for many tasks have been developed concerning quadrupedal robots. Katz [8] presented a system with back driveable modular actuators, which enabld high-bandwidth force control. Hutter [9] reported joint modules with integrated electronics that precisely control the torque. In both cases, the controllability, robustness, and importance of the motor-gearbox group during locomotion and the control of the angular position of the active joints were demonstrated, in addition to showing the kinematic, dynamic, and control strategies used. These terrestrial developments offers a different perspective for underwater exploration, with the respective technologies required, and performing (or improving) the same actions as AUVs and ROVs in inspection, sampling, and manipulation tasks in seafloor environments.
In this paper, a novel implementation based on quadruped anatomy is proposed. It has updates to its hardware, software, and movement control for terrestrial environment and offers an initial product for the underwater gait of a quadrupedal system. The quadruped robot can achieve different movement states (standby, walking, trotting, jumping, etc.) when performing its tasks, but in this study, only the normal gait state was analyzed.
Accordingly, this paper presents a simulation of the gait of a quadruped robotic system on seabed soil that must meet the required criteria regarding stiffness, contact damping, and static and dynamic friction in a specific environment with a terrigenous sediment seabed with its mechanical properties. We delimited the environment at a given depth of up to 10 m to avoid the effects of hydrostatic pressure [10], and incorporated the robotic mechanism in the soil model into the multibody simulation. We also included a previous analysis of the dynamic characteristics of the seabed with a simple mass-damping-spring-friction system, and then used these results to determine the torque required to overcome these seabed features and the quadrupedal gait state. The underwater forces' effects were not analyzed in this study, as the gait was achieved in a linear velocity lower than 0.2 m/s to mitigate the hydrodynamic damping effects (linear and quadratic), as mentioned by Fossen [11] for dynamic positioning speed regimes <2 m/s.
Often, a computational model depends on experimental data to complement and validate it. In the case in this study, through a method that shows all the phases for its implementation, from CAD design, kinematics, dynamic simulation of its motion and environment, and the implementation of control algorithms, we created a more realistic and easier-to-validate device that uses the appropriate hardware and software, with the model of its motor-gearbox group based on the simulation results for its prototyping and implementation with the Q 1 and Q 2 angular positions measured on the testbed.
The remainder of this paper is structured as follows: Section 2 presents the robot parameters, mechanical characteristics of the quadrupedal testbed robot, mechanism in each leg, geometry, and principal dimensions, as well as the inverse kinematics of the leg mechanism. Section 3 presents the seabed model, with the stiffness, damping, and friction coefficients calculated and represented as a block diagram in Simulink ® /SimscapeTM and another dynamic model created in Simulink ® to calculate the deep displacement in the same environment. Section 4 presents the prototype design with the gait simulation and results for the torque and angular position for each active joint in the simulation, gearbox design according to these results, and a leg prototype for a testbed, with the hardware and software architecture used for the low-level control of the motors. Section 5 presents the angular position and velocity (in PWM values) results for each active joint for the leg prototype in a testbed experimental gait, and the simulation results are compared with their equivalent mean squared error (MSE). Section 6 presents a discussion of the results. Section 7 presents the conclusions and recommendations for future studies.
All results will be used in future studies to develop a complete robotic system where other control algorithms can be tested and more easily validated in a real underwater test environment.

Mechanical Characteristics
The mechanical design was created in Autodesk Inventor ® and transferred to Simulink ® /Simscape in a multibody environment to apply kinematic control and achieve a simulated gait. The quadrupedal form was adopted to execute tasks and avoid underwater obstacles using the dynamic advantages of this form during gait [12] (high load-weight ratio, high precision, and high resistance).
The leg mechanism was 2 D.O.F.-5R planar with four links L i = 150 mm (i = 1, 2), one link L 3 = 300 mm, and another link L B = 380 mm in the base (see Figures 1a and 2). Each leg was connected to the body with a shoulder joint without movement; the main robot measurements are shown in Figure 1b. The proposed geometry offers the robot advantages for future developments for the lifting of loads, swimming, and navigation over the seabed.

Inverse Kinematics of 5R Leg Mechanism
The 5R mechanism had two active joints; the inverse kinematics indicated the articular positions of these from the Cartesian position P(x, y, z) [13,14] in the leg tip to the local frame origin O placed in its base (z = 0 for a planar system). In the corresponding literature, some kinematic models of planar 5R mechanisms have been reported [15][16][17]. Figure 2 shows the kinematic diagram for the computation; from the inventor geometry, the required angles and dimensions can be obtained. At point E, with trigonometric manipulations, we could obtain the following expressions: Moreover, for the α i and β i angles (i = 1, 2), we obtained: where E x is positive for α 1 and β 1 and negative for α 2 and β 2 . Then, from Equations (1) and (2), we obtained the Q 1 and Q 2 articular coordinates:

Seabed Model
The block diagram, for the seabed, was modeled with the Contact Force Library [18], with the sphere to plane force concept, including the penetration volume V P according to the sphere radius length (2R × 2R × R), knowing the value for the stiffness and damping, and for a continuous friction law, the kinetic friction and static friction coefficients. R was determined by the robot geometry in the leg tip (R = 0.031 m). In [19], the preconsolidation pressure for a terrigenous sediment seabed was P C = 34.47 kPa. We determined the largest overburden pressure that could be exerted on the soil without unrecoverable volume change. With these data, and the area of the robot leg tip, from its geometry A tip = 7.6 × 10 −4 m 2 , the limit applied force into the seabed soil could be found with Equation (7): then, this force is F limit = 26.2 N.

Dynamic Seabed Model: Mass-Spring-Damper-Friction
For the linear force law, the seabed can be modeled as a second-order system with an elastic spring and damper behavior, with friction in the leg tip traveling into the seabed, in accordance with the pseudo-rigid-body model approach [20]. Figure 3 shows the equivalent model for the mass-spring-damper-friction system (hereafter referred to as the M.S.D.F. system), and Equation (8) gives the mathematical model used in the subsequent analysis.  The force F(t) is applied with a sinusoidal input from the robot weight divided by four (from the Inventor model) with a value W Leg = 47.175 N (equivalent to the bias) and added force F e = 1 N for eventual external disturbance (equivalent to the wave amplitude). Comparing F(t) with F limit , we concluded that there will be deformation in the seabed soil, which can be obtained with the M.S.D.F. system modeling. The soil mass was found in [21], where the firm mud density is ρ mud = 1580 kg/m 3 . Then, the value for the penetration volume for R = 0.031 m is V P = 0.00012 m 3 , and the firm mud mass could be found: m mud = 0.1896 kg. For friction coefficient data, see Section 3.2. The static friction coefficient was not considered in the dynamical analysis.
Equation (8) was represented in a block diagram in Simulink to obtain its velocity response in time. Figure 4 shows the diagram model of the system with the respective variable and coefficient data. The velocity results had a sinusoidal form and were used with an integrator block to calculate the penetration depth in the seabed in the time domain, obtaining a maximum depth P D = 9.78 × 10 −8 m in firm mud. These results (Table 1) showed that the soil reference plane width in the simulation was >P D .

Contact Force Parameters
Contact stiffness coefficient (k): This is the reaction force that the seabed applies per unit penetration depth (displacement u) in the contact area. In [22], for a chain line body, the stiffness value per unit length of the line is (k/l) = 20 N/m 2 , and the corresponding coefficient for the contact depth of the sphere to plane concept for R is k = 0.62 N/m.
Contact damping coefficient (c): This is the relation between the normal force applied when the object travels into the seabed (not when it is coming out) and the velocity normal to the seabed. In [23], using the phase-shift method in a hardened soil model with small strain stiffness, the vertical damping coefficient, C (Ns/m), was calculated from 150 to 600 kN loads at 0.6 Hz (4 rad/s; see Section 4.1: Input signals). Then, with these data and a polynomial regression method, we obtained a value of C = 1.166 × 10 7 Ns/m for a leg weight W leg = 0.047175 N as a load. Kinetic and static friction coefficients µ k and µ s : As in the above paragraphs, the sliding friction coefficient was shown in [24]. Between a chain and seabed, the dynamic friction coefficient µ d = 0.62, and the static friction coefficient µ s = 1.01 for firm mud conditions. Table 2 summarizes all contact force parameters in firm mud conditions. The damping and contact stiffness coefficients were modeled in the respective references in a hardened soil model with slight deformation stiffness. These showed only one model for firm mud and was considered as such in this study. Table 2's data were input in the Simscape contact forces interface to define the seabed model with a block diagram. In Figure 5, the seabed model block shows four inputs, each corresponding to the sphere located at the leg tip to achieve sphere-plane contact. The velocity threshold in the friction tag was 0.001 m/s to ensure the lock or unlock friction contact [25], and the plane depth to the reference frame was R/2 = 0.0155 m > P D .

Quadrupedal Gait Simulation
The multibody model was simulated on a CPU with a sixth-generation Core i7 processor running Windows NT 10.0, and a 32 GB RAM memory board, with Simulink release 2020b and configuration parameters of a maximum step size of 0.008 and solver ode 15 s with a Simscape inventor integration. It organized the kinematic control elements throughout a block diagram (see Figure 6) with five principal parts, described below.
(1) Input signals: for x, y sinusoidal type with specific parameters. These signals were obtained after many experiments and complyiedwith two conditions: establishing an ellipsoidal trajectory in the leg tip [26], and avoiding a singular position [27] in the theoretical workspace [28]. Table 3 summarizes all signal parameters.  The experiments were performed with different values for amplitude, bias, and frequency, taking the local frame origin O placed in its base (Figure 2, z = 0 for a planar system) in each leg. The amplitude determined the stride length on the x-axis, the bias determined the leg tip height on the y-axis, and the frequency determined the linear velocity of the quadrupedal gait. These three parameters are mentioned by Catavitello et al. [29] for six dog types that were 57 cm at the withers in height and a gait speed Vc = 0.2 m/s. The experimental results reproduced these paramters. The phase parameter was out of phase in pi/2 to ensure that the gait alternation in each leg was in accordance with [29] regarding the kinematic pattern for a swim range of motion in dogs. Then, using the sensing tool in the Simscape 6-D.O.F. joint block, the linear translation velocity could be measured (V = 0.2 m/s) in the overall quadrupedal robot.
(2) MATLAB function: A code for the inverse kinematic computation that also showed their singular values in the Q 1 and Q 2 articular positions, shown as a function block.
(3) Initial positions: For previewing the simulation as a standby state, initial values in the active joints were necessary. The shoulder joints were always in zero value, and their analysis was not considered in this study.
(4) Multibody system: This was defined by exporting the inventor CAD to the Simscape interface in Simulink and verifying the corresponding active joints in each leg.
(5) Seabed model: This was used for defining the soil parameters, e.g., contact stiffness, contact damping, and static and kinetic friction coefficients (see Section 3).
With these five elements, the simulation showed the 3D gait ( Figure 7). The model was simplified from the original CAD to obtain a faster response in the simulation computing time (see the video called "Robot Gait.mp4" on the Supplementary Materials section link).

Gait Simulation Results
The gait pattern shown in the simulation had two principal situations: the visualization of the gait and the seabed model's influence on it. These results were conducive to the analysis of the angular position of the active joint and the required torque to achieve the gait and exceed the requirements of damping, stiffness, and friction for the seabed model. These requirements are described below.

Angular Position and Velocity
For the gait state, the robot has many requirements to achieve the angular position in active joints Q 1 and Q 2 . For applying the input signals (x, y) on the inverse kinematic block function, the angular position was computed and is presented in Figure 8a,b, respectively. For active joints' angular velocity, the results are shown in Figure 8c,d. Thus, from these results, we could obtain the equivalent data for the PWM signals in the motor.

Torque
From the multibody simulation, the torque required in the active joint Q 1 was as shown in Figure 9a, and for active joint Q 2 , as shown in Figure 9b. Torque values for the design were the maximum amplitude data in Figure 9b, with T Q2 = 16.02 Nm.

Gearbox
Torque (T Q2 = 16.02 Nm) and angular velocity (ω 2 = 4.32 rpm) requirements are described in Section 4.2, and the waterproofing could achieved with a marine electric motor, selected from BlueRobotics Inc., located in Los Angeles, U.S., of the brushless type, for operation at 20VCC, and using an ESC-R1 controller [30]. Table 4 shows the principal characteristics required.   [30]. rpm = 1120 is the angular velocity of the motor, but, according to Figure 8d, for a robot gait, it was necessary to set this at 4.32 rpm in a transmission relationship i = 1/6 with three similar stages. The angular output velocity from the gearbox was ω gb = 5.18 rpm ≈ ω 2 = 4.32 rpm from the simulation. For the torque, with a motor power P = 10 W and using the angular velocity selected ω gb , with dimensional analysis, the final calculated torque was T gb = 18.44 Nm > T Q2 .

PWM (µs) rpm Current (A) Voltage (V) Power (W) Force (kgf) Effic. (g/W)
For the design of a gear drive stage (with i = 1/6), a planetary gearbox concept was used for compact sizing. From Inventor, the Spur Gear toolbox was used to design the gearbox stage, with a central sun gear wheel, three planetary gears, and a crown (ring gear). The parameters for the design are shown in Table 5, and the dimensions are shown in Figure 10. Each gearbox stage was printed with ASA material, assembled with stainless steel bolts and nuts, and driven by the M200 brushless motor.

Leg Prototype
The leg prototype was created with ASA-printed parts for the links of the 5R mechanism and the base link with a high-density polycarbonate sheet with holes for the gearbox shafts.
The leg assembly was created on two linear carriages (lineal and vertical) and on an aluminium profile frame with their respective linear guides [31]; then, they were coupled on a steel frame. Figure 11a shows the mounted prototype for the testbed (see the video called "Leg system move.mp4" on the Supplementary Materials section link).

Hardware and Software Architecture
The angular position of motor θ was sensed by an encoder module consisting of a magnetic actuator and a separate sensor board, with a chip mounted on the motor shaft to achieve low-level control. The selected encoder was RMB20IC-Incremental, from RLS ® [32], with a square wave differential line driver to RS422; see Figure 11b. The motor had a basic speed controller ESC-R1 unit, connected to an Arduino Mega 2560 board output for processing the encoder signal and controlling its angular position. Figure 12 represents the hardware and software architecture for the leg position control. The application worked with IDE Arduino V1.8.57.0, and state estimation was managed from a laptop CPU with a sixth-generation Core i7 processor running Windows ® 10 Pro, which sent a signal from the serial COM port to the Arduino board at 115,200 Bps.

Low-Level Control of the Motor
From reference [30], the equivalent PWM data were used to calculate a ninth-order equivalent polynomial regression (with MSE = 4.91 error) for angular velocityθ (in revolutions per minute) converted to PWM data with the form: The generated curve is shown in Figure 13. Inverse kinematic block results sent the angular position signal multiplied by a constant 6 3 (by three gearbox stages). This signal was used to obtain the angular velocity (see Figure 14). With Equation (9), the PWM data were obtained for the setpoint. The final PWM curves obtained had an approximate sine waveform with the parameters in Table 6. Notably, the input PWM parameter to the controller ESC was a result of a conversion in the time domain. Using it as an input to the motor also resulted in a "PW M controller value" as a function of time.
The implemented algorithm for the angular position control used the PID concept. It sought to keep the input variable close to the desired setpoint signal by adjusting the output by tuning the PID parameters. Then, using the Simulink auto-tuning tool, we obtained the closer desired output and implemented it on the Arduino program for K P = 0.24, K I = 0.265 and K D = 0.0029 for Q 1 , and K P = 0.18, K I = 0.27 and K D = 0.0042 for Q 2 .
The motor-ESC group worked as a black box (plant model) and complied with the classical closed-loop controller. Thus, by introducing a PWM input signal to the motor, an approximately angular velocity response was obtained at the active joint that the encoder sensed to feedback the loop [33]. Path planning control at the leg tip was not considered in this study. Figure 15 shows the closed-loop control system.

Results
Angular velocity control followed by angular position control in active joints were tested in an underwater environment (see Figure 16 and the video called "Testbed.mp4" on the Supplementary Materials section link).     6. Discussion

Gait Simulation
The simulation model provided an adequate approximation for determining the influence of many variables, positions, and geometry in the quadrupedal gait and then for calculating torque. The multibody model needed an adequate computer system with many hardware resources to solve its dynamics. Referring to the specialized literature helped us to identify many characteristics of the real environment and place them in the virtual model. The multibody systems could be simulated with other input signals to obtain other gait and movement states. Different authors have reported different parameters for damping and contact stiffness coefficients that can determine a simulation with different position and torque results. An alternative would be to perform the experimental calculation of these coefficients using soil models with materials similar to those determined in the literature. Regarding the quadruped type of gait, other states of motion can be obtained and used to achieve other goals in path planning and high-level control considering the dynamic effects.

Angular Position Control of Motors
From the simulation results, the normal gait state presented different modes for movement from the eight active joints (the shoulder joint was not analyzed in this study), but it could be normalized by sending the correct input signals to the motors. The M200 motor was probed as a thruster, in use on thousands of marine robotic vehicles around the world, with a speed controller, but not for angular position control. The angular position control results from the velocity control showed a different application in these novel robot systems. Moreover, the testbed angular position results showed unusual behavior in the maximum wave peak for Q 1 and at the minimum peak for Q 2 . From the MSE information, higher error occurred in the Q 2 position, but distortion could be caused by the backlash in the angular position of the motor gearbox system at the time of landing the leg tip on the soil; this hypothesis should be verified by future experimentation.

Velocity Control of Motors
The designed PID controller for the motor achieved an adequate amplitude and frequency response and could be tuned with other methods. The backlash effect could be considered a mechanical impedance and could produce angular position and velocity errors, as can be seen in Figures 19 and 20, through the PWM velocity response. In addition, the disturbances caused in the angular position signals may have been due to harmonics produced by the variable-speed electronic drive. The possible implementation of filters to eliminate parasitic signals should be analyzed.

Conclusions and Future Works
The design of an alternative system for quadrupedal gait was developed and tested in a specific underwater environment in which the hydrodynamic damping effects did not need to be considered.A single method was constructed that was contrasted with a computational model.
Initially, we built a simulation model with the most realistic parameters for the soil and quadrupedal gait throughout a multibody model that provided initial information about the gait pattern produced to introduce an adequate input signal.
The sine wave signal required (Table 3) for reproducing the normal gait was found after many experiments to avoid singular positions of the 5R parallel mechanism. We observed the relationship between the ellipsoidal trajectory characteristics described by the leg tip and the gait type, as well as the linear velocity and torque required in each experiment to reproduce a pattern that approximated a dog's gait.
The seabed model was simulated by introducing it into a multibody system with the mechanical properties and characteristics of the soil in analysis via information obtained from the related literature. The deep penetration length was obtained by a secondorder mass-spring-damper-friction system (M.S.D.F.) and yielded a maximum depth P D = 9.78 × 10 −8 m in firm mud.
From the gait simulation results, the angular velocity required was ω 2 = 0.4524 rad/s (4.3201 rpm), as shown in Figure 8d, and the torque results in Figure 9 presented higher peaks in Q 2 displacement with T Q2 = 16.02 Nm.
A motor-gearbox group (see Section 4.3) was designed, made, and probed for the specific requirements of angular velocity and torque in the underwater environment and could be improved for decreasing the backslash between gears. Furthermore, it can be considered a device with patent potential for new developments in the underwater exploration given its characteristics of small sizing, high torque, and low angular velocity.
The hardware and software architecture (Section 4.5) were designed according to lowweight and high-performance requirements to implement low-level control in each motor. The angular position control was achieved through velocity control with the manufacturer's data for the equivalent PWM to velocity in revolutions per minute and introducing a PID controller in a closed-loop configuration, as well as the use of an encoder as a sensor to transform the velocity in the angular position, which we compared with a sinusoidal setpoint provided by the simulation.
The angular position testbed results, shown in Figures 17 and 18, indicated an unusual perturbation in the contact point between the leg tip and soil, in the maximum wave peak for Q 1 , with a mean squared error MSE Q 1 = 1.743 × 10 −7 , and a minimum peak for Q 2 and MSE Q 2 = 5.4184 × 10 −7 .
Velocity control resulted in an adequate amplitude and frequency response and could be tuned with other methods and probed with other control models. A complete underwater quadrupedal robot can be developed based on these initial results for many tasks in this environment.
The influence of the shoulder joint should also be analyzed in future studies to improve the gait pattern and extend the research toward robot trajectory control (high-level control).
The 5R geometry provides an alternative for mobility in underwater environments if a membrane is added between its links for swim tasks. Funding: This study was supported by the Spanish Government Projects under grants DPI2014-57220-C2-1-P and PGC2018-095939-B-I00; in part by the RoboCity2030 DIH-CM Madrid Robotics Digital Innovation Hub, S2018/NMT-4331, funded by the Programas Actividades I + D en la Comunidad de Madrid; and in part by Structural Funds of the EU.