Dynamic Analysis and Design Optimization of a Drag-Based Vibratory Swimmer

Many organisms achieve locomotion via reciprocal motions. This paper presents the dynamic analysis and design optimization of a vibratory swimmer with asymmetric drag forces and fluid added mass. The swimmer consists of a floating body with an oscillatory mass inside. One-dimensional oscillations of the mass cause the body to oscillate with the same frequency as the mass. An asymmetric rigid fin attached to the bottom of the body generates asymmetric hydrodynamic forces, which drive the swimmer either backward or forward on average, depending on the orientation of the fin. The equation of motion of the system is a time-periodic, piecewise-smooth differential equation. We use simulations to determine the hydrodynamic forces acting on the fin and averaging techniques to determine the dynamic response of the swimmer. The analytical results are found to be in good agreement with vibratory swimmer prototype experiments. We found that the average unidirectional speed of the swimmer is optimized if the ratio of the forward and backward drag coefficients is minimized. The analysis presented here can aid in the design and optimization of bio-inspired and biomimetic robotic swimmers. A magnetically controlled microscale vibratory swimmer like the one described here could have applications in targeted drug delivery.


Introduction
Vibrations are a characteristic of many physical systems. Many biological and robotic systems are vibration-driven, i.e., they use periodic motions of their limbs, bodies, or links to generate forces for locomotion. Some examples are flying insects and birds, fish, snakes and worms, swimming micro-organisms, swimming robots, and flapping-wing biomimetic air vehicles. Vibrations may also be used for open-loop stabilization or to change the transient response properties of a system, which is called vibrational control [1][2][3]. The vibrational mechanics of smooth dynamical systems is a well-developed research area. But many vibrational systems are only piecewise-smooth. For example, mechanical systems with dry (Coulomb) friction, mechanical systems with aero-or hydrodynamic forces, electrical circuits with switches, mechanical systems with colliding bodies, and biolocomotion systems [4].
Most often, due to the presence of time-periodic inputs, vibrational systems are time-periodic, and hence time-varying systems. In general, dynamic analysis of time-varying systems is more complicated than for time-invariant systems [5]. However, for classes of time-periodic systems one may use the averaging theorem to transform the original time-periodic system into a time-invariant system which represents the averaged dynamics [6][7][8]. According to the averaging theorem, every stable equilibrium point of the averaged dynamics corresponds to a stable periodic orbit of the original, time-periodic system with the same stability properties of the equilibrium point, and which "hovers" in a small neighborhood around that equilibrium point. Averaging techniques are widely used for dynamical analysis of smooth, time-periodic systems. However in [9][10][11], the averaging theorem is also used for vibration analysis of piecewise-smooth, time-periodic systems. Although more complicated, higher order averaging techniques are also available for the analysis of time-periodic systems [12,13].
One application of vibrational mechanics and averaging is for the dynamic analysis of flying and swimming biolocomotion systems and for the analysis, control, and optimization of biomimetic robotic systems. Most of the available robotic swimmers and fliers use actuated tails, fins, wings, or flagella for motion [14][15][16][17]. Design optimization and optimization of kinematic gaits and inputs of these types of swimmers and fliers have been active research areas [18][19][20][21]. For a review of robotic fliers and macro-and micro-size swimmers see [22][23][24][25]. Dynamics, control, and optimization of mobile robots, such as robotic swimmers and fliers, with no actuated members outside of their bodies have also been subjects of vast research recently. This type of robots use internal rotors, internal oscillatory masses, or deformations of their bodies to generate asymmetric fluid or friction forces for motion [26][27][28][29][30][31][32][33].
In this paper we consider the vibrational analysis and design optimization of a drag-based vibratory swimmer. Unlike conventional robotic swimmers, the proposed swimmer does not have any actuated or moving parts on the exterior of its body, which causes its design, fabrication, and maintenance to be simple. The swimmer with mass M consists of a symmetric body which floats on the surface of a fluid, e.g., water. A mass m oscillates inside the body with constant amplitude and frequency. A rigid, thin, inclined plate (fin) is attached to the bottom of the body and is submerged in the fluid. Hydrodynamic forces act on the fin when it moves relative to the fluid. According to the principle of conservation of linear momentum, the swimmer body oscillates in response to the oscillation the mass m with the same frequency, but a different amplitude. Without the inclined fin, the average velocity of the swimmer is zero. By adding the inclined fin and thus breaking the symmetry of the external swimmer geometry, an asymmetric drag force is applied to the system. Due to the asymmetric drag force, the average velocity of the swimmer will not be zero anymore and the swimmer moves in one direction on average with a certain mean velocity. The fin also changes the total mass of the system by adding an asymmetric fluid added mass to the mass of the system.
We performed finite element simulations with COMSOL Multiphysics R and determined the drag coefficient of an inclined rigid fin attached to a wall and immersed in a fluid, for different values of the fin angle, and established an approximate relationship between the drag coefficient and the fin angle. Using the averaged dynamics of the system and the relationship we found between the drag coefficient and the fin angle, we determine the optimum fin angle for the swimmer to move at maximum average speed. We present a general theorem for design optimization of a class of vibrational systems with asymmetric forces acting on the system. Through experiments, we verified the dynamic model, optimization technique, and finite element simulation results used for analysis. The results show adequate agreement between analytical and experimental design optimization. The idea presented in this paper may be used for the dynamic analysis and design optimization of macro-and microsized robotic swimmers. e.g., for microrobots used for targeted drug delivery. Though the swimmer presented in this paper floats on the fluid surface, the general class of systems analyzed also includes vibratory swimmers which are immersed in a fluid. For microscale swimmers, one may use micro-cantilever beams as the vibratory mass.
This paper is organized as follows. In Section 2.1 we discuss the averaging of a class of mechanical systems with high-frequency inputs. We use the averaging technique described in Section 2.1 to determine the average dynamics of the vibratory swimmer in Section 2.2, after discussing the equations of motion of the swimmer. Section 2.3 presents the design optimization of the vibratory swimmer to maximize its mean velocity. In Section 2.4 we discuss the results from numerical simulations to determine the hydrodynamic forces acting on a rigid, inclined plate submerged in a fluid. Section 3.1 discusses the averaging and design optimization of vibratory boat with both asymmetric drag and asymmetric fluid added mass forces. The swimmer used for experiments and experimental results are presented in Section 3.2. The main results of the research are summarized in Section 4.

Averaging of Piecewise-Smooth Systems with High Frequency Inputs
Consider the first order systemẋ where x is the n × 1 state vector, f (x) is a vector field containing smooth or piecewise-smooth functions, and u(t) = u 1 (t), . . . , u n (t) T is the input vector. Consider high-frequency, high-amplitude inputs of the form where v i (t) are zero-mean, T-periodic functions, and ω is the "high" frequency. Substituting inputs of the form (2) into system (1), the dynamic equations of the system arė where v(τ)dτ = ( v 1 (τ)dτ, . . . , v n (τ)dτ) T . Taking the derivative with respect to τ, (5) becomes Comparing Equations (4) and (6), one concludes that where Since the inputs v i (τ) are zero-mean, T-periodic functions, g(y, τ) is also T-periodic in its second argument. The initial conditions y 0 are determined using transformation (5) for τ = 0. System (7) is in the standard averaging form and its averaged dynamics are represented by [6,7] Replacing the fast time scale τ with the slow time scale t, the averaged dynamics (8) can be written in the form˙ȳ =ḡ(ȳ),ȳ(0) = y 0 .
According to the averaging theorem, if the averaged dynamics (9) possesses a hyperbolically stable equilibrium point, then the original time-periodic system (7) possesses a hyperbolically stable periodic orbit in an O( ) neighborhood of that equilibrium point. From transformation (5) it is evident that the average values of x and y are equal, i.e.,x =ȳ.

The Vibratory Swimmer: Dynamic Analysis and Averaging
In this section we derive the equation of motion of the proposed vibratory swimmer and determine its averaged dynamics. The vibratory swimmer considered in this paper consists of a symmetric main body floating on the surface of a fluid. A smaller body oscillates back and forth inside the main body. A rigid, inclined fin is attached under the swimmer and is submerged inside the fluid as shown in Figure 1. Here, we ignore the symmetric hydrodynamic forces acting on the main body and only consider the drag force acting on the fin. We also ignore the pitch and plunge motions of the main body and only consider its one-dimensional motion in the horizontal x-direction. Suppose that the oscillatory mass performs a harmonic oscillation of the form where y is the position vector of the oscillatory mass relative to the main body and Y ω and ω are the amplitude and frequency of its vibrations. The x-component of the linear momentum of the system is where M is the total mass of the system including the oscillating mass and the added mass of the fluid, m is the mass of the oscillating body, and x is the position vector of the body. Using Newton's second law, the equation of motion of the system is where F D is the drag force acting on the system. Replacingÿ using the expression above, the equation In this paper we neglect the symmetric drag force acting on the main body and only consider the drag force of the fin in the form where ρ is the fluid density, A is the fin's total surface area, v =ẋ is the velocity of the main body, α is the fin angle (see Figure 1), and C D (α) is the drag coefficient of the fin. It is evident that due to the asymmetry of the fin, the drag coefficient C D and the added mass of the fluid are different in forward (v > 0) and backward (v < 0) motion of the swimmer. In this section we consider a symmetric added mass, hence a constant M, and an asymmetric drag coefficient. We discuss the effects of asymmetric added mass in Section 3.1.
Consider the drag coefficient to be For simplicity we use c(α) = ρAC D (α) 2M andȲ = mY M , and write Equation (13) where using the drag coefficients C D p (α) and The coefficient c(α) can then be written as Replacing c(α) in Equation (16) and simplifying, the equation of motion iṡ Equation (19) is in the form of system (3) with n = 1 and a piecewise-smooth function. Using the fast time scale τ = ωt and the transformation v = u +Ȳ sin τ (20) and following the procedure discussed in Section 2.1, the averaged dynamics of system (19) are determined to be˙ū where φ = sin −1ū Y and the initial condition u 0 = v 0 . Figure 2 shows time histories of the original, time-periodic system (19) and its averaged dynamics (21). The physical parameters are M = 0.4 kg, m = 0.2 kg, A = 0.04 m 2 , ρ = 1000 kg/m 3 , Y = 0.2 m/s, ω = 20 rad/s, C D p = 0.1, and C D n = 0.2. The initial conditions are x 0 = 0 and v 0 = 0. By linearizing the averaged dynamics about their equilibrium pointū e = 0.0135 m/s, one can show that this equilibrium point is stable. Therefore the original time-periodic system possesses a stable periodic orbit in a small neighborhood of that equilibrium point.

Optimization of a Class of Vibratory Systems
In this section we consider optimization of a class of vibratory systems, including the vibratory swimmer presented in Section 2.2. For the case of the vibratory swimmer, the goal is to determine the fin angle α (see Figure 1) to maximize the average velocityv of the main body for a certain set of physical parameters of the system. Theorem 1. Consider the first order time-periodic systeṁ where f (x) is a smooth or piecewise-smooth function, u(t) is a T-periodic, zero-mean function, p ∈ D ⊆ R is a parameter, and where c p (p) and c n (p) are functions of class C r , r ≥ 1. Suppose that for any value of p ∈ D, the average dynamics of time-periodic system (22) possesses a hyperbolically stable equilibrium pointx e =x e (p). Then the extrema ofx e coincide with the extrema of c p (p) c n (p) .

Proof. See Appendix A.
The dynamics (16) of the vibratory swimmer is in the form of (22). Therefore, according to Theorem 1, the average velocityv of the main body is maximum if the ratio c n (α) is minimum, where α is the fin angle. In this case, while the main body is oscillating back and forth, it also moves forward with maximum speed on average. Therefore the maximum average speed of the swimmer depends on the ratio c p (α) c n (α) of the damping coefficients. Remark 1. From system (22) one concludes that if the drag force acting on the rigid fin is proportional to any positive power of the velocity, i.e., F D ∝ v r , r ∈ R + , the average velocity of the swimmer is still maximum if the ratio c p (α) c n (α) is minimum.

Numerical Simulation of Hydrodynamic Forces Acting on a Rigid, Inclined Fin
In this section we discuss the numerical computation of the drag coefficient of a rigid, rectangular fin immersed in a viscous fluid using COMSOL Multiphysics R software (The COMSOL Group, Stockholm, Sweden).The fin is attached to an infinite large, fixed surface parallel to the flow, such that the fluid does not flow over the top of the plate. The fin has an angle α with the infinite surface and the flow direction, as shown in Figure 3. The CFD module was used for the computations. No-slip boundary conditions were used at the walls, a constant velocity condition was used at the inlet of the domain, and a pressure boundary condition was used at the exit. Domain size and resolution studies were performed to ensure that the domain size of 14 m × 4 m and CFL number of 1.0 were sufficient to predict a converged value for the drag coefficient. We used two-dimensional numerical simulations to determine the drag coefficient C D of the fin for different values of α ∈ [0, π] rad. To determine the drag coefficient, we divided the total drag force found in the simulations by the dynamic pressure and the fin frontal area. The simulation results were used for design optimization of the vibratory swimmer. Figure 4 presents a snapshot of the simulations for α = 50 • after the simulation has reached steady state. The fluid inlet velocity is U ∞ = 0.1 m/s and the Reynolds number is Re = ρvL µ = 22,500, where ρ and µ are the density and dynamic viscosity of the fluid, respectively (assumed to be water at standard pressure and temperature here), and L is the length of the fin, which is 0.3 m in the simulations, with a thickness of 5 mm. Figure 5 shows the drag coefficient C D versus the fin angle α. For optimization purposes, we approximate the drag coefficient with the following 12 th degree polynomial in terms of α, where α is in radians C D (α) ≈ 0.018α 12 − 0.33α 11 + 2.583α 10 − 11.48α 9 + 31.74α 8 − 56.38α 7 +63.98α 6 − 44.5α 5 + 17.31α 4 − 3.766α 3 + 1.44α 2 + 0.132α  Using Theorem 1 and the approximate function of C D (α), one determines the optimum fin angle α m ≈ 9.4 • which generates the maximum average velocity of the swimmerv m ≈ 8 mm/s, determined using Equation (21). Figure 6 presents the displacement of the swimmer for the three values of fin angle α 1 = 5 • , α 2 = α m , and α 3 = 15 • . Figure 7 shows the average displacement of the swimmer versus time for the parameter values listed in the caption. The rest of the physical parameters are the same as used in simulations to generate Figure 2. It can be seen that though during the transient motion at the beginning of motion the optimum angle α m may not generate the maximum velocity, during steady state motion the swimmer with the optimum angle moves with the maximum speed. Therefore, after a "long enough" time, the swimmer with the optimum fin angle travels farther than the swimmer with any other angle.

The Vibratory Swimmer with an Asymmetric Added Mass
In the vibratory swimmer introduced in Section 2.2 and depicted in Figure 1, the flow of the fluid around the asymmetric inclined fin is not identical during forward (v > 0) and backward (v < 0) motions. The asymmetric flow of the fluid in the forward and backward motion of the swimmer generates both asymmetric drag force and asymmetric added mass of the fluid. In this section we discuss the averaging and design optimization of the swimmer when considering realistic asymmetric drag forces and added masses. This is the form of the model that we will compare with experimental results in Section 3.2.
Consider the equation of motion (13) of the system with the total mass of the system M(α) = M 0 + M f (α), where M 0 is the mass of the system without considering the fluid added mass, hence it is constant, and M f (α) is the fluid added mass which depends on the fin angle α and the direction of motion of the swimmer. Therefore, the total mass of the system is the total mass can also be written as M(α) =M(α) +M(α) sgn(v). Defining a(α) = 1 2 ρAC D (α), the equation of motion of the system iṡ ω cos ωt (25) whereā(α) andã(α) are determined using a p (α) = 1 2 ρAC D p (α) and a n (α) = 1 2 ρAC D n (α) similar to (18). Using the fast-time scale τ = ωt, Equation (25) can be written as whereȲ p (α) = mY M p (α) ,Ȳ n (α) = mY M n (α) , and r(α) = M p (α) M n (α) =Ȳ n (α) Y p (α) . By taking derivative with respect to τ from (27) and comparing the result with Equation (26), one gets M p (α) and b n (α) = a n (α) M n (α) . The time periodic Equations (28) and (29) are in the standard averaging form. Using the averaging theorem, the averaged dynamics of the system is dx dτ where t 1 = π + sin −1ū Y p (α) and t 2 = 2π − sin −1ū Y p (α) . After some math and replacing the fast-time scale τ with the slow-time scale t, the averaged dynamics can be written as where φ = sin −1ū Y p (α) . Knowing the fluid added mass function in terms of the fin angle α, i.e., M f = M f (α), one can use the averaged dynamics for design optimization of the swimmer, e.g., determine the angle α for maximum speed of the swimmer on average (v). For this goal, one determines α for: (1)ū e to be the equilibrium point of the second equation of the averaged dynamics (31), and (2)ū e maximizesv. Therefore ifū e is an equilibrium point of the averaged dynamics, from Equations (31) one can writē v e = F(ū e , α) where F(ū, α) and G(ū, α) are the right hand side of the first and second equations of (31), respectively. Then the problem reduces to a constrained optimization problem; maximizev in (32) under constraint (33). Using the method of Lagrange multipliers, we use the function The equilibrium pointū e that maximizes or minimizesv e , and the corresponding angle α are then determined from the following set of equations Since the determination of the fluid added mass is complicated and outside the scope of this paper, in order to show the effectiveness of the averaging and optimization results in this section, we use a hypothetical added mass and write the total mass of the system in the form where α is in radians. Using M 0 = 0.4 kg, the total mass M of the system in terms of α is presented in Figure 8. With the rest of the physical parameters being the same as in Section 2.4, using Equations (35) and (31), one determines the optimum angle α m ≈ 8.7 • and the maximum average velocity of the swimmerv m ≈ 0.67 cm/s, respectively. The simulation results using three fin angles α 1 = 5 • , α 2 = α m , and α 3 = 13 • and zero initial conditions are shown in Figure 9.

Experiments
To validate the theoretical and numerical results presented in previous sections, we performed experiments with a prototype vibratory swimmer. Here, we describe the device we used in the experiments and present the experimental results. The device used for the experiments consists of a small swimmer floating in a water tank, as shown in Figure 10. A rigid rectangular fin with an adjustable angle is attached to the bottom of the swimmer and is submerged in water. A small servomotor in the swimmer moves (vibrates) a mass back and forth, causing the swimmer to move back and forth with the same frequency. A battery is used to power the servomotor. For the main body we used a rectangular piece of polystyrene. The mass of the swimmer is M = 104 g, the vibratory mass is m = 43 g, and the rigid fin is 11.5 cm × 4.5 cm. The amplitude and frequency of oscillations are 3 cm and 2π rad/s, respectively. Interested readers can watch Video S1 in the Supplementary Materials to see the prototype swimmer in operation. To perform an experiment, we chose a test fin angle, started the system with zero initial speed, and measured the time required for the swimmer to travel a certain rectilinear distance. This determined the swimmer's average velocity corresponding to that fin angle. To minimize the effects of disturbances and uncertainties, we repeated the experiment at least five times for each fin angle. We also chose the maximum distance from the starting point possible as a stopping point for the experiments. Figure 11 shows the experimental average swimming versus fin angle compared with the average swimming speed predicted by the theoretical model from Section 3.1. Though the theoretical results predict the maximum average velocity happens for an optimal fin angle of α m ≈ 10 • , the experiments show that the true optimal angle is α m ≈ 20 • . Since the real flow is three-dimensional while the modeled flow used in numerical simulations for determining the drag coefficient is two-dimensional, and considering the other effects neglected in the modeling, such as the hydrodynamic forces acting on the main body, the discrepancy displayed between the theoretical and experimental results is not unexpected. Although the trends are the same in both cases, including the predicted maximum average swimming speed, which differs by only about 7% between the experiments and theory, the experiments reveal a higher optimal fin angle than that predicted by the two-dimensional theoretical analysis. From Equation (12) it is evident that the inertial force generated by the acceleration of the sliding mass, i.e., mÿ, causes the back and forth motion of the swimmer. For a swimmer with symmetric drag and added mass, if acceleration of the sliding mass is also symmetric during one cycle, the swimmer does not move on average. However, if the acceleration of the sliding mass, and hence the inertial forces, are not symmetric during one cycle, the swimmer moves in one direction on average. Note that the acceleration can be a zero-mean, periodic function of time, however not symmetric over one cycle. The slider-crank mechanism used in the experiments generates periodic, zero-mean, asymmetric acceleration for the slider. Therefore, even though for the fin angles α = 0 and α = 90 • the drag acting on the system and the added mass are symmetric, the swimmer still possesses non-zero velocity on average. This phenomenon can be clearly seen in Figure 10 for α = 0 and α = 90 • .

Discussion
In this paper we introduced a mathematical model for a vibratory surface swimmer consisting of a symmetrical body containing an additional vibratory mass inside the body, and an asymmetric rigid fin attached to the bottom of the body and immersed in fluid. Due to the symmetric inertial forces generated by the oscillations of the vibratory mass and the asymmetric hydrodynamic forces acting on the fin, the swimmer moves in one direction on average while performing asymmetric oscillations. Using the averaged dynamics of a class of vibratory systems and constrained optimization, we determined the optimum design for the rigid fin to generate the maximum average speed for the swimmer. The optimum speed predicted by the model and found in the associated physical experiments we performed with a prototype vibratory swimmer different by no more than 7%, though the theory predicted an optimal fin angle of around 10 • while the actual optimal fin angle was found experimentally to be close to 20 • .
Generalizing our results beyond the specific parameters and geometries used here, we showed that for a class of vibratory systems with asymmetric damping or drag forces in their forward and backward motions, to maximize or minimize the equilibrium state of the averaged system, one should design the system to maximize or minimize the ratio of the damping or drag coefficients during the forward and backward motion. The class of systems considered in this paper includes both vibratory swimmers moving on a fluid surface (vibratory surface vessels) and vibratory swimmers that are completely submerged in a fluid (vibratory underwater vehicles). The results do not depend on assumptions of swimmer scale or the details of mass oscillation actuation. Therefore the results of this paper can be used for the design optimization of a wide range of underwater vehicles, biomimetic robotic swimmers, and swimming microrobots, including those that could be used for drug delivery.