Numerical Analysis of a Drop-Shaped Aquatic Robot

: Finite-dimensional equations constructed earlier to describe the motion of an aquatic drop-shaped robot due to given rotor oscillations are studied. To study the equations of motion, we use the Poincaré map method, estimates of the Lyapunov exponents, and the parameter continuation method to explore the evolution of asymptotically stable solutions. It is shown that, in addition to the so-called main periodic solution of the equations of motion for which the robot moves in a circle in a natural way, an additional asymptotically stable periodic solution can arise under the influence of highly asymmetric impulsive control. This solution corresponds to the robot’s sideways motion near the circle. It is shown that this additional periodic solution can lose stability according to the Neimark–Sacker scenario, and an attracting torus appears in its vicinity. Thus, a quasiperiodic mode of motion can exist in the phase space of the system. It is shown that quasiperiodic solutions of the equations of motion also correspond to the quasiperiodic motion of the robot in a bounded region along a trajectory of a rather complex shape. Also, strange attractors were found that correspond to the drifting motion of the robot. These modes of motion were found for the first time in the dynamics of the drop-shaped robot.


Introduction
To describe the free and controlled motion of rigid bodies in a fluid, finite-dimensional mathematical models based on ordinary differential equations are often used.It should be noted that a full description of motion requires simultaneously solving the Navier-Stokes equations and the equations of motion of the body [1,2].However, finite-dimensional models have a number of advantages.To study them, the well-developed analytical and computational methods of the theory of dynamical systems can be used, which makes it possible to study the dynamics of the system entirely.In addition, finite-dimensional models allow a real-time prediction of system dynamics, compared to infinite-dimensional models based on partial differential equations.
Almost all finite-dimensional dynamical systems describing the motion of rigid bodies in a fluid have forms similar to those of the Kirchhoff equations [3] or the Chaplygin equations [4].The above-mentioned equations of motion are supplemented with terms describing forces and torque due to the viscous resistance of the medium, and the presence of any singularities of the flow such as sinks/sources or vortices.In the latter case, the equations of motion of a rigid body are usually supplemented with equations for the transfer of point vortices.We recall that a point vortex is a mathematical abstraction that assumes that the vorticity is concentrated at one point in the plane [5], in contrast to a physical vortex that occupies a certain region of space, for example, a tornado [6].This makes it possible to obtain, instead of a partial differential equation governing the transport of vorticity, ordinary differential equations for the coordinates of point vortices.Note that the equations of motion of point vortices are first-order equations [5], in contrast to the equations of motion of material particles [7].
In the case of a smooth body, it is assumed a priori that the vortices are always present in a fluid and have constant intensities; see, for example [8][9][10][11][12].For bodies with a sharp edge, it is possible to describe the generation of vorticity based on the point vortex model and the Kutta-Chaplygin condition [13,14].Such mathematical models were considered, for example, in [15][16][17][18][19][20].
Note that the model with the shedding of point vortices is easier to compute than the joint numerical solution of the Navier-Stokes equations and the equations of body motion.However, adding a new vortex increases the dimension of the phase space of the system by 2, thus making it impossible to study the dynamics of the system qualitatively.Therefore, models that have a small fixed dimension of phase space and describe the influence of vorticity indirectly in the coefficients of the equations are more convenient for qualitative analysis of motion.See, for example [21,22].
In addition to point vortices, one considers the influence of other types of point singularities on the motion of a rigid body.In [23], the motion of a balanced and an unbalanced circular foil in the field of a fixed-position point source was studied.
The study of finite-dimensional models is also of interest for the theory of dynamical systems.Despite their simplicity, such models can exhibite various dynamical effects: asymptotic stability with respect to part of variables [24,25], Fermi acceleration under the periodic controls [26][27][28][29], stabilization or loss of stability due to periodic disturbances [30], chaotic scattering of trajectories [23,31], the transition to chaos according to Feigenbaum's scenario [32,33] and through a series of doublings of the invariant tori [34][35][36], etc.
In this paper, we continue the study of the mathematical model developed in [22], which describes the motion of an aquatic robot under the periodic rotation of a rotor.In [37], it was shown that for this model the fastest locomotion is achieved when control becomes impulsive.Note that periodic controls close to impulsive ones also occur in practical implementations.For example, the work [38] describes a fish-like robot with an elastic tail that has two stable positions.The transition between these positions occurs almost impulsively and generates propulsion (https://www.youtube.com/watch?v=DZ8SU_00 BLU&ab_channel=ZechenXiong, accessed on 9 December 2023).
Section 2 presents the equations of motion and the control law.Section 3 shows that, in the case of a highly asymmetric impulsive control law, additional stable and unstable periodic and quasiperiodic modes of motion and strange attractors can occur in the system.These modes of motion do not occur in the case of admissible controls specified by the technical data of the engine.

Equations of Motion
Let us consider the equations governing the plane-parallel motion of an aquatic robot [22] shown in Figure 1a.
In this paper, we keep the notation used earlier in [22,37]: v 1 , v 2 are the projections of the linear velocity of point C on the axis of the moving coordinate system attached to the robot (see Figure 1b), ω is the angular velocity of the robot, k = I r Ω(t) is the angular momentum of the rotor (the control), Ω is the angular velocity of the rotor, m = 0.905 kg is the mass of the robot, I = 0.00844 kg • m 2 is the central moment of inertia of the robot, and I r = 0.00058 kg • m 2 is the central moment of inertia of the rotor.The values of the moments of inertia I and I r were calculated based on a 3D models built in a CAD system.The coefficients λ k ij and c k were calculated in [22] by processing the data of experimentally determined trajectories of the robot and have the following values: The technique for processing data and calculating coefficients (3) is described in Appendix A.

The Form of Control
The equations of motion (1) contain the derivative of the angular momentum of the rotor k or its angular acceleration Ω.In the previous work [22], it was noted that rotation of the rotor with constant velocity is equivalent to the motion of a system with a fixed (frozen) rotor.In addition, a shift in the angular velocity of the rotor by a constant Ω → Ω + const does not change the form of the equations of motion (1).
As before, we will control the motion of the robot by alternating between uniform and uniformly accelerated rotations of the rotor.In particular, we will use the following dependence of the angular velocity of the rotor: where Ω max is the constant that specifies the maximum angular velocity of the rotor, t 1 and t 3 are the lengths of time intervals of uniform rotation of the rotor, t 2 and t 4 are the lengths of time intervals of uniformly accelerated rotation of the rotor, and A graph of the dependence ( 4) is shown in Figure 2a. Figure 2b shows a graph of the derivative of the angular momentum of the rotor k(t).Remark 1.As a rule, piecewise constant controls turn out to be convenient for practical implementation and research (ignoring the transient processes due to switches between pieces of the control).
In particular, such an approach was used in [39] to study parametric resonance.
In [37], in studying Equations ( 1) and ( 2) we restricted ourselves to the values Ω max ⩽ 50 s −1 and the time of uniformly accelerated rotation of the rotor equal to no less than 0.3 s.These restrictions are due to the technical data of the engine.Since this work is theoretical, we ignore restriction on the time of accelerated rotation of the rotor, tending it to zero, and we limit ourselves to studying the motion at We consider a question that is important from a theoretical point of view: how does impulsive control influence the dynamics?In this case, t 2 = t 4 → +0 and a δ-singularity appears in the third equation of motion (1): Thus, the values v 2 and ω have a jump when the rotor changes the direction of rotation impulsively.In this case it is necessary to integrate the second Equation ( 1) and the Equation ( 7) in the neighborhood of the δ-singularity over the interval tending to zero.For times t = τ n 1 = τ n 2 , the changes of v 2 and ω are given by the relations , and for times t , Since the system (1) includes dissipative terms, the following statement can be proven similar to one in [40]: For large values of v 1 , v 2 , and ω, the kinetic energy of the system decreases, that is, the vector field specified by Equation ( 1) is directed inside some compact region.Thus, there is always an attractor in the system.

Limit Cycles
In the previous work [37], it was shown that under technically feasible control actions the system (1) has a single asymptotically stable limit cycle.This limit cycle corresponds to directed propulsion in the case of a symmetric control law (t 1 = t 3 , t 2 = t 4 ) and to the motion near a circle in the case of an asymmetric control law (t 1 ̸ = t 3 , t 2 = t 4 ).
The question arises: do limit cycles different from those mentioned above and more complex modes of motion exist?Each fixed point of the period advance Poincaré map for the period corresponds to a periodic solution of the system (1).We have performed a scan of the phase space of the Poincaré map within with steps It turns out that for sufficiently small values of t 2 = t 4 and strong asymmetry of control (t 1 ≪ t 3 ) the Poincaré map has more than one fixed point.
Using the parameter continuation method for the detected fixed points, we constructed bifurcation diagrams, examples of which are given in Figure 3a-c.These diagrams show the dependence of the coordinate v 1 of fixed points (main p 1 and additional p 2 , p 3 ) on the parameter t 1 for various fixed values of t 2 .For coordinates v 2 and ω, bifurcation diagrams have a similar form and are not presented here.From Figure 3a-c it is clear that additional fixed points p 2 and p 3 arise in pairs due to tangent bifurcation.The bifucarcation value t 1 increases as t 2 decreases and reaches its highest value at t 2 = 0.
Figure 3d-f show the dependence of the fixed point multipliers on the parameter t 1 for t 2 = 0. We see that the main fixed point p 1 is always stable (Figure 3d), the additional fixed point p 3 is always unstable (Figure 3f), and the fixed point p 2 , which is asymptotically stable at the moment of birth, loses stability as the parameter t 1 decreases (Figure 3e).
Stable fixed points correspond to the motion of the robot near a circle.Moreover, on the limit cycle corresponding to the fixed point p 1 , the robot moves in a "natural" way (see Figure 3g).And for a fixed point p 2 , sideways motion that is almost a circle takes place (see Figure 3h).
Figures 4a-c and 5a-c show the corresponding dependencies v 1 (t), v 2 (t), ω(t) at t 1 = 0.055 for both stable fixed points.Figures 4d and 5d show fragments of the trajectory of the robot.We see that the trajectories are piecewise smooth due to the jumps in the velocity component v 2 .(a-c) Evolition of variables v 1 (t), v 2 (t), ω(t) and (d) a fragment of the trajectory of the robot corresponding to the fixed point p 2 of Poincaré map with the parameter control value t 1 = 0.055.Due to impulsive nature of selected control, the curve v 1 (t) and trajectory are continuous and piecewise smooth, and curves v 2 (t) and ω(t) have jumps.

Neimark-Sacker Bifurcation: Evolution of Invariant Curves of the Poincaré Map
In the bifurcation diagram (Figure 3c) we saw that the fixed point p 2 loses stability as the parameter t 1 decreases.The loss of stability occurs according to the scenario of the Neimark-Sacker bifurcation [41,42].In this case, an attracting invariant curve appears in the neighborhood of a fixed point when it becomes unstable.
To make sure that a Neimark-Sacker bifurcation does occur, we need to verify according to [43] a number of conditions for the moment of bifurcation:

1.
No strong resonances To check this condition, we calculate arg µ 1,2 before bifurcation at t 1 = 0.04985 and after the moment of bifurcation at t 1 = 0.0498 we have Assuming that arg µ 1,2 changes fairly smoothly as t 1 changes, we see that there are no strong resonances.

2.
Strictly monotonic increase in |µ 1,2 | in the vicinity of the bifurcation value This condition can be easily verified in Figure 3c.We see that at the moment of bifurcation |µ 1,2 | intersects the line |µ 1,2 | = 1 transversally.Therefore, condition (13) is satisfied.

3.
An additional condition is associated with the construction of the normal form (which is problematic in numerical calculations) and allows us to establish the type of bifurcation.However, the type of bifurcation can be determined using the numerical solution.Since a stable invariant curve arises, the bifurcation is supercritical.
Next, we will study the process of evolution of invariant curves of the Poincaré map using the parameter continuation method.In particular, we will change the values of t 1 from 0.0498 down to values close to zero.Next, we will perform a similar procedure, increasing the parameter t 1 from almost zero values to the value 0.0498.Examples of invariant curves are shown in Figure 6.Additionally, we will estimate Lyapunov exponents (see Figures 7 and 8).To estimate Lyapunov exponents we used Benettin's algorithm [44,45].
The invariant curve resulting from the Neimark-Sacker bifurcation has a fairly simple shape (see Figure 6e,f).As t 1 decreases, this curve grows and disappears at t 1 ≈ 0.03989.At this moment, its Lyapunov exponents take values close to zero (see the red lines in Figures 7 and 8).
At t 1 ≈ 0.04233373687 an additional invariant curve of a rather complex shape appears in the phase space of the Poincaré map (see the black lines in Figure 6d).As the parameter t 1 increases, this dynamical mode disappears (see Figures 7 and 8).As the parameter t 1 decreases, the shape of the curve changes insignificantly (see Figure 6b,c).In addition, instead of this invariant curve, a strange attractor may arise (see Figure 6a) in very small intervals t 1 (see the jumps in the largest Lyapunov exponent in Figures 7 and 8).
6.The red dots correspond to the modes obtained by continuation with respect to the decreasing parameter t 1 .The black dots correspond to the modes obtained by continuation with respect to the increasing parameter t 1 .

Quasiperiodic Nature of the Robot Motion in the Case of an Invariant Curve
The invariant curves shown in the previous subsection correspond to invariant tori in the phase space of the system.The motion on the invariant torus is quasiperiodic and can be formally represented as the following multiple series (in the resonant case Ω 1 /Ω 2 ∈ Q the motion becomes periodic): The rotation angle of the robot can be obtained with direct integration (14): where ⟨ω⟩ = c ω 0, 0 is the mean value of the angular velocity, and ω(t) is the so-called quasiperiodic antiderivative.
Let us represent the equation for the trajectory of the robot in complex form In the last expression, the sum is a quasiperiodic function with frequencies Ω 1 and Ω 2 .Therefore, it can be represented as Let us perform term-by-term integration (17), assuming that (Ω, k) + ⟨ω⟩ ̸ = 0, ∀k ∈ Z 2 , i.e., there are no resonances.Integration yields the following expression for the trajectory: We see that the function Z(t) is three-periodic and, due to the absence of resonances, describes a compact curve on the plane (X, Y).Remark 2. Under the resonance condition (Ω, k) + ⟨ω⟩ = 0, the mean motion of the robot will be directed with the velocity c z k .The search for such modes of motion can be performed numerically but is not within the scope of this paper.
Let us construct a period advance map (under period T) for the function Z(t): Recall that T is the period of the control action.
If the frequencies Ω 1 and Ω 2 were Diophantine, then the series on the right side of (19) would represent a periodic function.In this case, the trajectory of the map Z n will consist of repetitive fragments rotated relative to each other by an angle ∆αn * , where n * satisfies the relation In the case of non-Diophantine frequencies Ω 1 and Ω 2 , relation (20) will be approximate: Ω 2 Tn * ≈ 2πm * .This means that the trajectory of the map Z n will consist of fragments similar in shape and rotated relative to each other by an angle ∆αn * .
Let us illustrate this using numerical experiments.Figure 9d shows an invariant curve of simple shape that exists at t 1 = 0.0399196451225613.The frequency Ω 2 can be found using the Fourier transform of the dependencies v n 1 , v n 2 , and ω n .The frequency ⟨ω⟩ will be defined as In the numerical experiments we have used the value N = 10 6 .Analysis of the numerical solution shows that We see that the frequencies Ω 1 = 2π, Ω 2 , ⟨ω⟩ are non-Diophantine and that the trajectory Z n is quasiperiodic (see Figure 9f) and consists of similar fragments (see Figure 9e).
A similar situation is observed for invariant curves of more complex shape, see Figure 10.In addition to the base frequency Ω 1 = 2π, this mode of motion is characterized by the following additional frequencies: Thus, we see that all frequencies are non-Diophantine.The trajectory Z n (see Figure 10f) consists of similar fragments of a more complex shape (see Figure 10e).
The behavior of the trajectory of the map Z n described above is not valid for strange map attractors.We see that the behavior of the variables v n 1 , v n 2 , and ω n is irregular (see Figure 11a-c).The trajectory Z n is drifting in nature (see Figure 11f) and consists of fragments that have no similarities (see Figure 11e).

Conclusions
In this work, we have presented a numerical study of the equations constructed in [22] to describe the motion of an aquatic robot with an internal rotor.We have studied the impulsive control law of the rotor when it instantly changes the rotation velocity.This case is interesting for a number of reasons.Firstly, for such a control the fastest locomotion of the robot is implemented [37].Secondly, within the framework of the considered model under a strong asymmetry of the control law (the time of rotation of the rotor clockwise is much longer/shorter than the time of counterclockwise rotation of the rotor), asymptotically stable periodic and quasiperiodic and chaotic modes of motion coexist.In fact, multistability and complex dynamics occur that have not been found previously in the dynamics of the drop-shaped robot.
An asymptotically stable mode of motion has been found, which is not the robot's self-propulsion in forward direction, but sideways motion in a circle.Of significant interest is the study of various designs of drop-shaped robots for which such motion modes are implemented in practice.In addition, this study shows that chaotic modes of motion can occur, which should be avoided in practice when developing this kind of robot.
It is important to keep in mind that the model used in this work is extremely simplified.Although it allows one to analyze the dependence of the system dynamics on the parameters, all the results obtained must be verified experimentally.We also note that it would also be possible to use a model of the motion of the robot in a viscous fluid using the Navier-Stokes equations, but this approach only allows for simulation of certain motions and provides no insight into the dependence on parameters.
where u 1 and u 2 are the projections of the fluid velocity vector onto the curvilinear axes, p is the pressure, ρ is the fluid density, ν is the kinematic viscosity, and w 1 = v 1 − ωx 2 (ξ, η), For Equation (A4), kinematic boundary conditions were specified at the foil boundary: Figure A1 shows a fragment of the computational mesh.
Figure A1.Computational mesh constructed using the complex variable boundary element method [47].
Given distributions u 1 , u 2 , and p, the forces F 1 , F 2 and the torque G acting on the body from the fluid are determined with the following integrals along the foil contour: (A6) Thus, after solving the auxiliary hydrodynamic problem, the sequence of values (A2) is supplemented with the following values: (A7) We see that the required coefficients λ (k) ij and c k are included into the expressions (A3) linearly and all other quantities are known at the current stage.Thus, λ (k) ij and c k can be found using the classical least squares method.The found values are referred to one meter, since the problem of plane-parallel motion is being solved.The required values of the coefficients are obtained as follows: where 0.0335 is the height in meters of the wetted part of the robot in the experiment.

Figure 1 .
Figure 1.(a) A fish-like aquatic robot.(b) OXY is a fixed coordinate system, Cx 1 x 2 is a moving coordinate system.

Figure 2 .
Figure 2. (a) Time dependence (4) of the angular velocity Ω(t) of the rotor; (b) time dependence of control k(t) corresponding to the dependence (4).

Figure 3 .
Figure 3. (a-c) Bifurcation diagrams of fixed points of the Poincaré map.The solid lines correspond to asymptotically stable motions, and the dashed lines correspond to unstable ones.(d-f) Multipliers of the fixed point at t 2 = 0. (g) Trajectory and position of the robot corresponding to the fixed point p 1 at t 1 = 0.055.(h) Trajectory and position of the robot corresponding to the fixed point p 2 at t 1 = 0.055.In panels (g,h) the coordinates x, y are measured in meters, and the robot is shown not to scale.Arrows show the direction of motion.

Figure 4 .
Figure 4. (a-c) Evolition of variables v 1 (t), v 2 (t), ω(t) and (d) a fragment of the trajectory of the robot corresponding to the fixed point p 1 of Poincaré map with the parameter control value t 1 = 0.055.Due to impulsive nature of selected control, the curve v 1 (t) and trajectory are continuous and piecewise smooth, and curves v 2 (t) and ω(t) have jumps.

Figure 5 .
Figure 5.(a-c) Evolition of variables v 1 (t), v 2 (t), ω(t) and (d) a fragment of the trajectory of the robot corresponding to the fixed point p 2 of Poincaré map with the parameter control value t 1 = 0.055.Due to impulsive nature of selected control, the curve v 1 (t) and trajectory are continuous and piecewise smooth, and curves v 2 (t) and ω(t) have jumps.

Figure 7 .
Figure 7. (a) The largest Lyapunov exponent λ 1 .(b) The middle Lyapunov exponent 2 .(c) The smallest Lyapunov exponent λ 3 .The red curves correspond to the modes obtained by decreasing the parameter t 1 .The black curves correspond to the modes obtained by increasing the parameter t 1 .

Figure 8 .
Figure 8. Fragment of the diagram shown in Figure 7. Strong peaks in the panel (a) correspond to appearance of strange attractor.Jumps of exponents λ 2 , λ 3 near value 0.0415 (see panels (b,c)) correspond to appearance additional invariant curve.

Figure 11 .
Figure 11.The variable of the map: (a) v n 1 (b) v n 2 , (c) ω n .(d) Portrait of the strange attractor.(e) Fragments of the trajectory Z n corresponding to the intervals between peaks v n 1 are depicted by different colors.(f) The trajectory of the map Z n at t 1 = 0.0399120977758897.Confidence interval of the highest Lyapunov exponents ⟨λ 1 ⟩ ≈ 0.00669 ± 0.00016 with reliability 0.999.
ωx 1 (ξ, η) are the components of the transfer velocity.The Lamé coefficient D and the terms β 1 , β 2 arising due to the curvature of the grid lines have the form: