On wave-driven propulsion

Abstract A theory is presented for wave-driven propulsion of floating bodies driven into oscillation at the fluid interface. By coupling the equations of motion of the body to a quasipotential flow model of the fluid, we derive expressions for the drift speed and propulsive thrust of the body which in turn are shown to be consistent with global momentum conservation. We explore the efficacy of our model in describing the motion of SurferBot (Rhee et al., Bioinspir. Biomim., vol. 17, issue 5, 2022), demonstrating close agreement with the experimentally determined drift speed and oscillatory dynamics. The efficiency of wave-driven propulsion is then computed as a function of driving oscillation frequency and the forcing location, revealing optimal values for both of these parameters which await confirmation in experiments. A comparison with other modes of locomotion and applications of our model with competitive water sports is discussed in conclusion.


Introduction
From flagellar beating of spermatozoa to water-walking insects, from the wave-like lateral flexions of fish to flying, nature has evolved myriad strategies for locomotion and propulsion in fluidic environments (Hu et al. 2003;Bush & Hu 2006;Bühler 2007;Gaffney et al. 2011;Eloy 2012;Smits 2019).A perhaps lesser-known propulsion mechanism is enacted by honeybees (Apis mellifera) trapped on the surface of water: in oscillating its wings, the bee generates a fore-aft asymmetric wavefield that contributes in part to forward motion (Roh & Gharib 2019).Inspired by the striken honeybee, Rhee et al. (2022) designed SurferBot: a centimeter-scale interfacial robot that achieves wave-driven locomotion speeds of around ∼1 cm/s by inducing asymmetric vibrations with a small eccentric motor.Wave-driven propulsion may also be realized by designing floating bodies with mass-and/or geometric-asymmetries that generate thrust and torques by exchanging momentum with waves generated by the object as it oscillates on the free surface of a globally vibrating fluid bath (Ho et al. 2021;Barotta et al. 2023).At much larger scales, Benham et al. (2022) achieved wave-driven locomotion speeds of around ∼1 m/s by forcing a canoe into oscillations by jumping up and down on its gunwales, a technique known as gunwale bobbing.In this case a simple waveequation model was derived by treating the canoe as an oscillating pressure source with a prescribed motion, ignoring the effects of dispersion and considering an idealised Gaussian canoe shape.Although outside the scope of the present work, other studies have investigated alternative routes to wave-driven propulsion including symmetry-breaking due to the presence of boundaries (Tarr et al. 2023) and using acoustic radiation forces (Sabrina et al. 2018;Roux et al. 2022;Martischang et al. 2023).
Figure 1: Thrust scaling F T (1.1) plotted against drag scaling F D (1.2) for a variety of different bodies oscillating at the water surface (forces are taken as per unit width).Parameter values for each case are given in Table 1 in Appendix A. We observe that thrust balances drag with a ratio of 0.001 − 0.1.
Based on classical results (Longuet-Higgins & Stewart 1964;Longuet-Higgins 1977), the prevailing opinion is that the thrust necessary to drive locomotion originates from an asymmetric radiation of wave momentum (Ho et al. 2021;Rhee et al. 2022).A scaling argument for two-dimensional wave-driven propulsion (i.e.restricted to the vertical plane) is as follows: Consider a body of length L oscillating at the interface of a fluid with density ρ and kinematic viscosity ν in the presence of gravity g.If the oscillations have amplitude A, frequency ω, and wavenumber k, the thrust due to oscillating waves (per unit width) scales † like where v = Aω is the speed associated with the oscillations (Longuet-Higgins & Stewart 1964).Forward motion at speed U results in an inertial drag force that scales like where C D is a drag coefficient.By considering a number of known examples of wavedriven propulsion (Fig. 1), we demonstrate that thrust (1.1) balances drag (1.2) with a ratio of 0.001 − 0.1, indicating the range of prefactor values associated with (1.1).
Whilst the foregoing scaling argument performs well in capturing the thrust-drag balance, it fails to provide further insight into the problem.For example, there is no way of predicting the prefactors of the scalings or the related efficiency of propulsion.Further, without a model of the wavefield, momentum conservation cannot be verified.We herein present a theoretical model for wave-driven propulsion in the case of a raft floating at the fluid interface undergoing small-amplitude oscillations due to an external force.The equations of motion for the raft are coupled to a quasi-potential flow model of the fluid dynamics (Dias et al. 2008).Propulsion is demonstrated in the form of a constant drift speed that results from a time-averaged thrust force due to the oscillations.We compare our predictions for the wavefield and raft motion with the experimental data of SurferBot (Rhee et al. 2022), and then use our model to optimize the efficiency of propulsion by varying the frequency of vibrational forcing and the position at which Figure 2: (a) A raft of length L oscillating on the surface z = η(x, t) of a fluid of depth H self-propels with a drift-velocity U thanks to a self-generated wavefield.(b) Schematic diagram of the raft dynamics (in the moving frame).The raft is subject to an oscillatory external force F A = (F A,x (t), F A,z (t)) applied at position x A (in the frame of the raft).
The applied force gives rise to small horizontal and vertical oscillations ξ(t) and ζ(t), in addition to planar rotations with angle θ.
the forcing is applied.Some possible extensions and applications of our model to insect locomotion and competitive water sports are discussed in conclusion.

Theoretical model for an oscillating raft
Whilst scaling arguments for wave-driven propulsion show some success in predicting drag and thrust (see Fig. 1), there is no existing theory to model wave-driven dynamics that in turn can be used quantitatively to predict and maximise the efficiency of propulsion.In the following section we build such a model, first by formulating the equations of motion for a floating raft in two dimensions, then by accounting for the coupled fluid flow problem using quasi-potential flow theory.

Equations of motion for the raft
Let us consider a two-dimensional body of fluid with the horizontal and vertical coordinates given by x and z, respectively.A floating raft of length L resting on the fluid surface is subject to an applied periodic force F A = (F A,x (t), F A,z (t)) imposed at position x A = (x A , z A ) measured from the raft centre of mass (Figs. 2(a) and (b)).In response to F A it is assumed that the centre of mass of the raft X undergoes small oscillations in the x-and z-directions, whilst translating horizontally with an a priori unknown constant drift speed U .The action of F A also induces a torque that gives rise to oscillations in the x − z plane, the rotation vector given by θ = θ(t)j where the unit vector j points out of the page.
The equations of motion for the position and orientation of the centre of mass of the raft are (2.1a) where m is the mass (per unit width) of the raft and I = mL 2 /12 is moment of inertia of the raft (per unit width), taken as the value for a rod rotating about its centre of mass.
The second term on the right-hand-side of Equation (2.1a) is the force arising from fluid pressure, p, minus its atmospheric value p a , where the vector n is the unit normal to the raft surface, S, pointing in the positive z direction.The drag force is approximated using an empirical drag coefficient such that 2) The torques on the right-hand-side of Equation (2.1b) are those arising from their corresponding forces in Equation (2.1a).Drag on rotational motion is assumed to be negligible.We note the omission of a gravitational term in (2.1a) which is balanced by the pressure integral term in the steady state.While we do not go into details, the steady contribution to the pressure is a combination of surface tension (acting at the edges of the raft) and Archimedes buoyancy.We assume that the mass of the raft is such that this balance is sustained and hence we ignore the gravitational term from the momentum Equation (2.1a).
Consistent with experiments (Rhee et al. 2022), the applied force is taken to be periodic and can be written as where FA,x , FA,z are small complex constants (incorporating both amplitude and phase).
In response, the position of the centre of mass of the raft is where (2.5) represent small horizontal and vertical oscillations, respectively (Fig. 2(b)).The profile of the raft is thus described by where (2.7) Likewise, the fluid pressure reads p = p a + pe iωt , (2.8) where p = p(x, z) is a complex function.
Before going further, we pause to discuss the steady drift speed, U .The hypothesis of this paper is that steady drift is a result of a balance between the propulsion due to waves and the resultant inertial drag.We propose that the propulsion due to waves is given by the time-averaged x-component of the pressure integral term in (2.1a) and that this thrust force is given by where the overhead-bar notation denotes time-averaging and i is a unit vector in the horizontal.As shown in Appendix B, the time-average of the drag term (2.2) is wellapproximated by the component due to steady drift U , while drag due to oscillations is comparatively small.Hence, the drag in the x-direction † approximates to From the above analysis, it is clear that the time-averaged thrust (2.9) and drag forces (2.10) are both second order.In other words Consequently, by taking a square root of the thrust-drag balance (F T = F D ), we see that U = O(ǫ).
To fully determine U by equating (2.9) and (2.10) requires a fluid dynamics model for the pressure p in Equation (2.9).Likewise, the fluid pressure determines the firstorder dynamics of the raft.Specifically, after substituting (2.3), (2.5), (2.7) and (2.8)  into Equations (2.1) and linearizing, the first-order oscillating dynamics are described by Equations ( 2.11) provide a relationship between the constants FA,x , FA,z , x A , and the constants ξ, ζ, θ, given p calculated by formulating the corresponding fluid flow problem.

Quasi-potential flow
The body of fluid is assumed to be irrotational and infinitely wide with a finite depth H and a free surface at z = η(x, t), such that −∞ < x < ∞ and −H z η.Following the approach of Dias et al. (2008), the velocity for a weakly damped † flow can still be modelled as the gradient of a potential function φ(x, z, t) to good approximation.Due to incompressibility, φ satisfies (2.12) The linearised kinematic condition applies on the water surface, such that On the raft the linearised kinematic condition takes the form In addition, the unsteady version of Bernoulli's equation applies within the fluid which, upon linearisation, becomes (2.15) After applying the dynamic boundary condition p − p a = −γη xx on the free surface, (2.15) yields where Laplace pressure introduces the surface tension γ.Finally, the bottom surface is impermeable, such that (2.17) In the same manner as §2.1, we seek solutions of the form φ = φe iωt , η = ηe iωt , where φ = φ(x, z), η = η(x) are complex functions.
Henceforth, we formulate the fluid flow problem in the moving reference frame, x → x + U t.However, we restrict our attention to the case where the drift speed is much smaller than the typical wave speed, such that U ≪ √ gL.In this way, we focus on the flow resulting from the oscillations, rather than the drift velocity.
(2.18e)For the purpose of numerical simulation, the domain is rendered finite of length ℓ in the x direction.Radiative boundary conditions (Equation (2.18d)) are chosen at the left-and right-hand boundaries to avoid reflection, where k satisfies the dispersion relation The final step in linking together the fluid-raft interaction is calculating the thrust force F T (2.9) in terms of the pressure and normal vector, noting that to leading order n • i ≈ −η x .From (2.15) the pressure on the raft is given by p = −ρ(iω φ + g η + 2ν φzz )| z=0 . (2.20) Meanwhile, on the free surface the wavefield is given by the solution to the kinematic condition (2.13), such that which is accompanied by radiative boundary conditions (similar to (2.18d)) at x = ∓ℓ and continuity boundary conditions on the raft, namely η = ζ ± θL/2 on x = ±L/2.

Summary
Let's now summarise the governing system of equations.Given constants FA,x , FA,z and x A , the dynamics of the raft are determined by Equations (2.11), coupled to the fluid flow problem given by Equations (2.18) via the pressure (2.20).The drift speed U is then calculated by equating thrust (2.9) and drag (2.10), where F T follows from the solution via p and η.The problem is solved numerically in MATLAB using a combination of finite differences † to solve the PDE system (2.18) and Newton's method to resolve the coupling with (2.11).Example code is provided in the Supplemental Materials.

Model results, validation and application to SurferBot
We begin this section by exploring the numerical results of our model, demonstrating that it is satisfies a conservation of momentum balance.We then consider a specific application of our model to SurferBot (Rhee et al. 2022) comparing our results directly with experimental data.The velocity potential is plotted in Fig. 3 in the case of pure pitching (a) and pure heaving (b).As anticipated, F T = 0 for both these cases since both heaving and pitching are required to generate thrust ‡.
Next, we validate our numerical scheme using a conservation of momentum balance.As shown in Appendix C and first discussed by Longuet-Higgins & Stewart (1964), an expression for the thrust force can be derived via integration of the Euler equations.This results in an equivalence between F T and the momentum flux of the form The term on the right-hand side is the difference in mean momentum-flux between backward and forward directions (including the pressure difference).Hence, when larger waves propagate backwards than forwards, this results in a positive thrust.In Appendix C we show how to calculate this balance in the case where surface tension and viscosity are neglected.Calculating the left-and right-hand sides of Equation (3.1) numerically for a range of ω yields very close agreement between the two (Fig. 4), giving us confidence that the numerical scheme is accurate and consistent with global momentum conservation.
It should be noted that whilst our model accurately conserves momentum, it neglects any possible contribution to thrust from vortex shedding since the fluid is considered irrotational.It is well-known that some forms of insect locomotion rely on both vortex shedding and capillary wave propagation (Bush & Hu 2006;Bühler 2007).To incorporate thrust from any vortices generated by the motion would require a model for rotational fluid flow, which we leave for a future study.However, by comparing our model with SurferBot in the next section, we will see that our predictions perform well despite ‡ Although both heaving and pitching are required to generate thrust in the case of vanishingly small drift velocity, it is not clear if this remains true at large drift velocities in general (Benham et al. 2022).neglecting vorticity, suggesting that in the present case vortices are not the leading-order driver of propulsion (Barotta et al. 2023).

Comparison with SurferBot
To demonstrate a specific application of our model, we calculate the dynamics of SurferBot (Rhee et al. 2022), for which dimensional parameter values are given in Table 1 in Appendix A. In this case, the force is applied via an eccentric motor positioned at x A = −3 mm from the centre of the raft.The forces F A,x , F A,z , are assumed to be out of phase, and are estimated using a scaling law relating the mass (per unit width) of the raft m = 2.6 g/3 cm and the oscillation amplitude A = 150 µm, such that ( FA,x , FA,z ) = mω 2 A (i, 1). (3.2) This results in a drift speed of U ≈ 2 cm/s, compared with an experimentally measured speed of ∼ 1.8 cm/s.A comparison between the experimental and theoretical selfgenerated wavefield is plotted in Figures 5(a Overall, the results of our theoretical model exhibit good qualitative agreement with experiments and yields an excellent prediction of the drift speed, U .However, there is disagreement in the amplitude of the wavefield and hence the oscillation amplitude of SurferBot, in addition to the decay length-scale of the waves.However, these discrepancies are likely attributed to some obvious differences between experiment and theory, including the precise form of the oscillatory forcing, three-dimensional effects of both the SurferBot geometry and wavefield, or indeed the flexure of the SurferBot body contributing to wave damping.

SurferBot
Present theory

Efficiency and optimization
Furnished with a theoretical model describing the SurferBot dynamics, we can interrogate the efficiency of propulsion and optimize the parameters of SurferBot thereafter.The efficiency of propulsion is given in terms of the applied and useful (propulsive) power.The total applied power is given by the sum of linear and angular contributions †.The propulsive power is taken as P T = F T U and the efficiency is defined as χ = P T /P A .Hence, in the case of SurferBot (before optimization), we calculate an efficiency value of χ = 1.8% from our mathematical model.
The two parameters we attempt to optimize are the driving frequency ω/2π and the horizontal position of the forcing, x A .We assume that the horizontal force FA,x = 0 (as it doesn't contribute to forward thrust anyway) and that the total power (4.1) is fixed ‡ at the value calculated for SurferBot, P A = 7.2 × 10 −4 W/m.Efficiency is plotted as a function of ω/2π and x A in Figure 6 showing an optimum frequency of 16 Hz and an † There is also a contribution FA,x ξ to P A which vanishes because FA,x and ξ are out of phase due to (2.11a).
‡ The applied force FA,z must change to maintain constant power at variable frequency.optimum motor position of 5 mm behind the raft centre, both of which could be readily tested in a future experiment.
In a similar manner to the discussion in §1, it is also possible to derive a very rough estimate of efficiency using scaling arguments.To do so we may assume that the total power P A (4.1) can be split into the contribution from propulsive power, which scales like P T ∼ F T U , and the power radiated away by the waves, which scales like P wave ∼ F T c g , where c g = dω/dk is the group velocity of the waves.The efficiency then scales like where Ma = U/c g is Mach number.In the case of SurferBot we find χ ∼ 3.5%, which in broad agreement with the model calculation χ = 1.8%.However, a scaling argument obviously provides no means for optimization.Nevertheless, it indicates that may increase with Mach number, which is probably due to a decrease in forward-propagating waves as the motion approaches the supersonic limit (which we discuss further shortly).We also note that due to the scaling (4.2), it is estimated that subsonic flows Ma < 1 have a maximum possible efficiency of χ ∼ 1/2, which we indicate with dashed lines in Figure 6.

Discussion and perspectives
Using linear quasi-potential flow theory, coupled with the equations of motion for an oscillating raft, we have demonstrated how the raft propels itself forward using its own waves.The close agreement between our linear model and the experimental data of SurferBot (Rhee et al. 2022) indicates that we have captured the key physics at play for this self-propelling centimetre-scale robot.Nevertheless, there are several further improvements to the model that are worth discussing.
The assumption of a vanishingly small drift velocity is worth revisiting.In a future study a finite drift-speed could be incorporated into the fluid-flow problem.This would require transforming Equations (2.18) to a moving reference frame, thereby introducing new derivative terms of the potential.Consequently, the dispersion relation (2.19) would need to be modified and may result in multiple-wavenumber solutions for forward and backward travelling waves.An initial resolution to this problem could involve studying a small but finite drift speed using the method of asymptotic expansions.
A finite drift-speed incorporated into the model would result in a Doppler-shifted wavefield, as discussed by Benham et al. (2022).In this case, a key dimensionless parameter is the Mach number, Ma, which is defined as the ratio between the drift speed U and the group velocity of the waves c g .The Mach number is closely related to whether or not forward (in addition to backward) waves are propagated.Specifically, if the drift speed is larger than the wave speed for a particular wavenumber, this indicates that such a wave can only propagate behind the object and not in front of it.Since the momentum balance (3.1) indicates that larger forward waves reduce the forward thrust, we would expect more efficient motion at higher Mach numbers (which is consistent with the scaling The value of the Mach number and hence the of the wavefield varies considerably for different cases of wave-driven propulsion.example, the wavefield of SurferBot (Ma = 0.04) is markedly different to that of gunwale bobbing on paddleboard (Ma = 1.7) †.The former produces significant waves both forwards and backwards, whereas the latter has purely backward waves (e.g.compare Figure 3(a) of (Rhee et al. 2022) with Figure 1(a) of (Benham et al. 2022)).
By considering a finite drift-speed, one can then study how symmetry is broken when oscillations commence (i.e. the start-up problem for wave-driven propulsion).One could investigate, for example, how perturbations grow or decay from a stationary state when subjected to oscillations.This would then help determine the necessary conditions for initiating forward propulsion under different operating conditions.Therein lie numerous optimization problems, such as the optimal control of a single oscillator (or multiple oscillators) with variable positions or masses, not to mention the optimization of the shape of the raft itself.One could also consider different optimization objectives, such as maximising speed as well as efficiency.
A rotational fluid model could be developed to account for propulsion due to both vortex shedding and wave propagation.Whilst the present irrotational model performs well in the case of SurferBot, it is well-known that vortices play a significant role in the locomotion of water striders which are smaller and faster than SurferBot (Bühler 2007).At larger scales (e.g. in the case of gunwale bobbing) it is not known what role vorticity plays.However, by developing a rotational fluid flow model, one could establish the regimes in which vorticity can be neglected as a means of propulsion.
One could extend the model to account for a three-dimensional flow field rather than the two spatial dimensions studied here.In such a model a more realistic representation of the SurferBot geometry (or other examples) could be rendered.This could then capture the anisotropic wavefield emitted by the rectangular-shaped SurferBot and what effect this has on the thrust calculations studied here.Further, the model described herein assumes that the raft is subject to an external, local force, motivated in part by the physical design of the SurferBot and the eccentric mass motor that drives propulsion.An interesting extension to the present work would be to consider another form of symmetrybreaking wherein a raft with variable mass density along its length self-propels due to the global vibration of the fluid bath (Ho et al. 2021;Barotta et al. 2023), a scenario akin to walking droplets (Couder et al. 2005;Bush 2015).
In closing, it is worth discussing some of the potential applications of this work.
Firstly, as indicated by Rhee et al. (2022), studies on wave-driven propulsion may serve as a way to better understand how insects move when floating on the water surface.Secondly, boat oscillations during rowing and canoe races (due to the strokes of athletes) may affect performance due to wave-body interactions.For example, Dode et al. (2022) studied how fluctuations in the horizontal speed may affect rowing efficiency.However, in the case of the vertical component of these oscillations, the impact on efficiency is uncertain et al. 2022).Hence incorporating a finite drift-speed into the current model could be leveraged to optimize stroke styles to improve race times.Another extension of this work could include the interactions between multiple bodies.This may find application in the context of ducklings benefiting from their mother's wavefield (Yuan et al. 2021), or in rowing/canoe races where boats may gain an advantage by interacting with the wavefields of their opponents.

Appendix C. Momentum balance
We start by neglecting viscosity and writing down the two-dimensional (nonlinear) Euler equations for incompressible fluid flow in the presence of gravity, which are By adding together (C 1) (multiplied by u) and (C 2) we get This can be integrated vertically to give where we have used the impermeability and kinematic boundary conditions This equates to Equationn (2.9) upon linearisation of the normal vector.
Next, we apply the foregoing calculation to a potential flow example in the absence of viscosity and surface tension, γ = ν = 0, since the calculation is much simpler.In the case of an oscillating raft, as described in the main text, the pressure within the fluid is given by the (linearised) Bernoulli's equation p = p a − ρ(gz + φ t ).
(C 12) However, on the free surface z = η (outside the raft region), the pressure is set by the dynamic boundary condition p = p a .Hence, the force (C 11) becomes The term on the left hand side of (C 10) approximates to We confirm this balance by numerical calculation in Figure 4 of the main text, which shows extremely close agreement between the right and left hand sides of (C 15).These calculations are for a fixed raft amplitude A/L = 1 (which is achieved by setting ζ/L = −1/2, θ = 1), whilst varying the frequency.All force values are normalised by the scaling ρgL 2 for convenience.We also compare the numerically calculated thrust values with the scaling (1.1), incorporating the difference in amplitude between aft and fore waves.This provides the force scaling where we use a prefactor of 3/4 to be consistent with other authors (Ho et al. 2021).Since this scaling is attributed to Longuet-Higgins, we compare it to the numerically calculated values in Figure 4 of the main text, labelled as the L-H scaling.The L-H scaling captures the thrust qualitatively but not quantitatively, further indicating the need for the model derived in this study. REFERENCES .10) Since the drag term (2.10) is approximately the same as it would be for steady flow past the raft, we use Blasius' theory of flow past a flat plate to estimate C D = 1.33 • Re −1/2 in terms of the Reynolds number Re = U L/ν (Schlichting & Gersten 2016).

Figure 4 :
Figure 4: Comparison between the numerically calculated thrust F T and momentum flux across the domain (Equation (3.1)) showing close agreement between the two.Comparison is also made with the scaling proposed by Longuet-Higgins (see Appendix C) exhibiting only qualitative agreement.Surface tension and viscosity are neglected for the purpose of these calculations.
) and (b), whilst the time-varying position of the back and front of SurferBot are shown in Figures5(c,e) (experimental) and 5(d,f) (theoretical).

Figure 5 :
Figure 5: Comparison between the experimental wave amplitude from (a) SurferBot and (b) the mathematical model.In (b), the theoretical wave amplitude (red) is scaled by 2L/πx to be consistent with the far-field behaviour of Bessel functions (since in practice the waves extend radially to the far-field).The unscaled wave amplitude is shown as a dotted black line.(c,d,e,f) Comparison between the experimental and theoretically predicted position of the aft (c,d) and fore (e,f).In the theoretical predictions (d) and (f), x aft = X − (1/2, θ/2)L and x fore = X + (1/2, θ/2)L.Experimental data reproduced with permission.

Figure 6 :
Figure 6: (a,b) Optimization of SurferBot efficiency χ by varying the frequency and the motor position whilst maintaining constant total applied power P A .The optimum values (illustrated with black dots) are the outcome of a two-parameter optimisation, and hence the plotted curves are slices of the two-dimensional efficiency surface at its maximum.