The propulsive effect generated by fast travelling-wave vibrations

A prototype problem for the propulsion arising from a spatially-dependent source is analysed asymptotically. Such driving is a phenomenon that, while fairly common in biological settings, is less well understood in mechanical and engineering contexts. If two parallel plates bound a fluid and the plates are free to translate relative to each other, recent calculations have shown that vibrations applied to one plate in the form of travelling waves may induce motion of the other. The degree of the propulsion obtained is a function of both the wavelength and the phase speed of the imposed vibrations and there seems to an optimal wavelength for which propulsion is greatest. In order to shed light on these observations, here we undertake a formal asymptotic analysis of the process in the circumstance when the vibrations move quickly. It is shown how the flow structure takes on a form consisting of a bulk core layer supplemented by suitable edge layers attached to the two plates. The efficiency of the mechanism relies on the generation of a mean flow component whose properties are considered. It is found that there is excellent agreement between our analytical predictions and the previous simulations; in particular we are able to find the critical wavelength for maximum propulsion and explain the modifications that arise in the short and long wavelength limits.

The flow configuration. Travelling vibrations on the lower plate induce a velocity U top of the upper plate effect, which consisted of a pair of parallel plates that can translate independently of each other. If the gap between the plates is filled with incompressible fluid then the motion of one surface can induce movement in the other. It is the peristaltic process that underpins the mechanism by which distributed propulsion systems induce motion.
Our aim in this study is to shed further light on this process by making particular study of the case when the imposed vibrations travel quickly. To order to set the scene, depicted in Fig. 1 is the configuration investigated in HF22. We suppose there are two flat surfaces that are taken to be very long in the X -direction; furthermore they are parallel and separated by a distance 2h. The lower plate is vibrated so that the dimensionless slot geometry (based on h) is given by Y L (X, t) = −1 + A cos[α(X − ct)] and Y U (X, t) = 1; here the subscripts L and U refer to the lower and upper plates respectively. Furthermore, t denotes the time while A is the amplitude of the vibrations; these are taken to be of wavenumber α and to move at a wavespeed c. This motion of the lower wall induces a movement of the upper wall, which we suppose translates at a speed U top .
If we write the governing continuity and Navier-Stokes equations in terms of the co-ordinates x = X − ct and y = Y that move with the wave phase speed what remains is the system where u=(u, v) is the velocity vector with components in the (X, Y ) directions. If ν denotes the kinematic viscosity, these velocities have been scaled on ν/ h. If the fluid density is ρ, we scale the pressure p on ρν 2 / h 2 and the time t is based on h 2 /ν. The shear stress that acts on the upper plate pulls it until it attains the velocity at which the mean shear stress vanishes. Accordingly, we require that v = 0, ∂u ∂ y mean = 0 on y = 1 ( 2 ) while at the lower surface we demand that Moreover, since our interest is with the relative movement of the plates as opposed to any effect that may arise owing to a mean pressure gradient, we ensure that the pressure field is periodic in the x-direction. We note that the definitions adopted here mean that the Reynolds number appears in the boundary conditions rather than within the equations themselves. Finally, we emphasise that the speed of the top plate U top is not specified and is part of the problem to be tackled. The task then is to solve the system (1) subject to the boundary conditions (2,3) with the speed of the top plate u = U top emerging as part of the solution. Haq & Floryan [3] computed the form of U top as a function of the wavenumber α for a selection of wavespeeds c. Their method involved solving the governing equations to spectral accuracy and combined this with the idea of the so-called method of immersed boundary conditions. Extended accounts of this technique may be obtained from [4] or [5] but the key idea of the approach is the recognition that the rather intricate geometry of the lower boundary does not lend itself to a straightforward application of spectral techniques. To circumvent this problem, the computational area is embedded within a larger, more conveniently shaped domain. The system of equations are then solved in this extended region, and constraints introduced as alternatives to the boundary conditions on the wavy boundary. The range of applicability of this method may be enhanced using an over-constrained formulation [6].
Some results typical of those found in [3] are shown in Fig. 2 which illustrates the form of the quantity U top /A 2 as a function of α for various wavespeeds c. We note that for (a given) relatively small vibration amplitude, up to perhaps a few percent of the slot width, the speed of the top plate is very closely proportional to A 2 ; furthermore in practical applications the imposed vibrations are not large. The introduction of the vibrations yields values of U top /A 2 that assume a rather complicated structure when regarded as functions of α. It appears that a noticeable propulsive effect occurs at moderate to large values of c. It also seems that long wavelength vibrations are inefficient in this regard and, intriguingly, the propulsion does not increase monotonically with α, as might have been expected. Instead, there is some critical value α * at which the propulsive effect seems to be greatest and this phenomenon actually weakens as the vibration wavelength is reduced. It is only once α exceeds about two that the propulsion strengthens again resumes its generally upwards trend at a rate roughly proportional to α 2 .
The main objective here is to analyse the mechanisms that give rise to this general behaviour. We pursue a formal large-c asymptotic analysis of the problem in an effort to understand the subsequent size of the propulsive effect. The focus on large wavespeeds is not devoid of real interest since carefully positioned piezoelectric elements put on the side of a fluid slot can be used as devices to provoke oscillations in the flow. Such elements typically operate at high frequencies which then naturally lead to fast waves.
In order to explain the form of the structure that arises at large wavespeeds we organise the remainder of this article as follows. We commence in Sect. 2 by rewriting the problem relative to a more convenient co-ordinate system and then, in Sects. 3-5, undertake a systematic analysis of the governing system. We pull together the key components of our findings and present a few final remarks in Sect. 6.

The transformed equations
Problems posed in a domain with an irregular boundary often present a formidable analytical challenge. It is frequently beneficial to transform variables in a convenient way so as to regularise the geometry although, except in very rare cases, this will normally lead to a more involved set of governing equations and/or boundary conditions. Nevertheless, this is often preferable to having to deal with the original domain and is a price well worth paying. To this end we define and, in terms of these new co-ordinates, the field equations now become within these the coefficients F 1 -F 8 are defined by Written in terms of these new variables the relevant boundary conditions become It is the system (5) subject to the boundary conditions (6) together with the requirement that there is no mean pressure gradient along the slot that we shall investigate in the large-c limit. To that end, let us define the small parameter ε so that we also set the amplitude of the corrugation

The flow structure
We shall find that the flow adopts a three-layer structure comprising a core region together with a pair of wall layers at the bounding plates. A schematic of the configuration is shown in Fig. 3. It is remarked that across the central portion of the slot, indicated as region I , the ξ -dependent parts of the flow field are essentially driven by the motion of the lower plate. The solution in this core zone needs to be supplemented by wall layers adjacent to the two edges η = ±1 in order for the boundary conditions there to be satisfied. In the following sections we set out the details of the problem in the three separate zones and crossreference them to show how the various solutions link. Whenever this type of analysis is undertaken there is an unavoidable need to switch between the various regimes for the flow structure in one region provides information as to how to proceed elsewhere. This can lead to a very confusing account in which the focus seems to continuously change and the reader can easily be overwhelmed by a plethora of notation. In an attempt to minimise this problem we shall discuss the solutions in the three regions in separate Sects. 3-5 and indicate the connection points as we proceed. We trust that isolating the details of the various zones in self-contained parts is less confusing than might be the case lest we attempt to follow a narrative that inevitably would require extensive swapping between zones. Across the majority of the slot away from the edges η = ±1 we consider a flow field given by It is important to explain some of our notation adopted here. In these expansions those quantities with an overbar are mean flow parts; they are functions of η alone. The other terms do depend on ξ ; their size is motivated by the fact that the velocity component v = ε −1Ā α sin ξ on η = −1.

Zeroth order
Leading order balances in the governing equations (5) show that we can identify a solution of the system in which the function V 00 satisfies V 00 − α 2 V 00 = 0 where the dash denotes a derivative with respect to η. Given that we expect that V 00 → αĀ as η → −1 and V 00 → 0 as η → 1 we can deduce the solution This then fixes the other components so that We will have need to know the form of the pressure as η → ±1. If we define then V 0 (η) → Q sin ξ as η → −1 so that in this limit Similarly, in the limit η → 1 then (15) and this shows that that

First order
At first order we can write the vertical velocity component in the general form with similar expansions for U 1 and P 1 . It turns out that V 11 (η) and V 12 (η) satisfy the equations so that progress requires knowledge of the mean-flow correction quantityū 0 . The streamwise momentum equation simply gives that so it is a just a linear function of η. The condition emanating from the bottom wall layer (that is the form of (33) as Z → ∞) requires thatū 0 → 0 as η → −1. On the other hand, at this stage we do not know the behaviour ofū 0 as the upper plate η = 1 is approached. Its value can only be determined once we have found the mean stress on the upper wall and this requires us to undertake a careful analysis of the upper edge region I I I . Accordingly, for the moment we shall just suppose thatū 0 → B 0 as η → 1 whereupon we conclude that the required value of the constant B 0 will be tied down in due course. Given this result we can proceed. We notice that in order to match with the solutions we have developed in the following sections, it follows that these requirements follow from (35b) and (41) respectively. We can then solve for the necessary elements of the problem at this order. It follows that These quantities also contain second harmonic terms proportional to sin 2ξ and cos 2ξ but these are not required for the subsequent calculations so are not written down. Of more significance is the finding which is the solution of V 11 − α 2 V 11 = 0 subject to the boundary conditions as η → ±1 deduced from the requirements (19). Moreover α P 11 = −V 11 so that we conclude that as η → 1 so

Second order
We next turn our attention to the second order problem and, in particular, to the form of the mean flow correction termū 1 . It transpires that we shall not need to determine explicitly other terms at this order although we record that where the . . . indicate terms that are redundant for our purposes. The governing equations imply that the function V 21 (η) satisfies V 21 − α 2 V 21 = 0 while The mean terms give d 2ū and we need to solve this subject to the requirement that the solution match to the form of U ‡ 1M in the bottom layer as Z → ∞ (see Eq. 37). Again, in exactly the same manner as we remarked in the process of determiningū 0 , we cannot know the behaviour ofū 1 as η → 1 prior to a closer examination of the upper edge layer. Accordingly we require that where the appropriate value of B 1 will be determined. Then with the constants C 1 and C 2 fixed by imposing the boundary conditions. Routine manipulations show that moreover, for matching purposes it is helpful to note that Having determined the important flow quantities across the bulk of the channel, we now need to switch attention to the wall layers in order to determine the entire flow structure.

The lower layer I I: η = −1 + O(ε)
We need to supplement the main region by thin wall layers at η = ±1. The extent of these layers can be motivated by reference to the governing system of equations (5). In the bulk of the channel it is the inertial terms in the equations that drive the flow and the viscous terms are subdominant. Since we have the wavespeed c = ε −2 it is an elementary task to verify that if the wall layer is of depth then the viscous terms balance the streamwise c(∂u/∂ξ ) term when ∼ ε.
At the bottom wall we therefore anticipate a thin layer of width O(ε). We define the O(1) co-ordinate Z using and the solutions then take the forms

The zeroth order problem
The leading-order governing system becomes Direct matching with the core layer behaviour (14) gives which enables us to solve for U ‡ 0 so that The form of V ‡ 0 then can be deduced by integrating the continuity equation subject to the requirement that V ‡ 0 vanishes on Z = 0. Then we have We remark on the behaviour of these solutions in the limit Z → ∞ which corresponds to moving back into the core. It may be seen that in this limit.

The first order problem
For the purposes of matching with the main layer solutions we only require the mean part U ‡ 1M of U ‡ 1 . The streamwise momentum equation at this order is whose mean components can be integrated to determine We conclude that At the top of the channel we put η = 1 − εz; we adopt this form as the fluid then occupies the region where z > 0. Here the flow variables take the forms Now the leading order balances give that while the pressure term P 0 is independent of z; by matching with (16) we deduce that The streamwise momentum equation is solved to give where the constant B 0 arises from the need to match with the core solutionū 0 (18) as η → 1. It should be remembered that this constant is not yet known. The form of V 0 follows from an integration of the continuity equation subject to the requirement that V 0 vanishes on z = 0. Then 5.1 The first order problem

The mean flow
The continuity equation at first order tells us that If we denote the mean part of V 1 by V 1M , the ξ -independent component of (42) yields which vanishes exponentially as z → ∞. The streamwise momentum balance at the same order becomes from which we can determine the mean component U 1M . It follows that whereupon we conclude

The ξ -dependent component
Now we can switch our focus to the remainder of the top layer at this order. There are components of the solution proportional to cos 2ξ and sin 2ξ but these higher harmonics are not be required for the subsequent work. The cross-slot momentum equation gives the simple result that ∂ P 1 /∂z = 0 so that P 1 = P 1 (ξ ). Its value can be inferred by matching with the main layer solution P 11 so that where the value ofP 11 is as stated in result (22). Given this result, we can now solve the Eq. (44) for the components of which are given by here we have defined the scaled co-ordinate z ≡ γ z. We can deduce the equivalent parts within the definition of on integration of the continuity Eq. (42) we determine

The second-order problem
The ξ -independent part of the streamwise velocity U 2 is rather complicated. The momentum equation gives an expression for d 2 U 2M /dz 2 which can be integrated once. This leads to the conclusion that We have now assembled all the results required to determine the speed of the upper wall induced by the travelling waves. We next describe how we can use the requirement for zero mean stress on the upper wall to derive the size of the propulsion effect. wavespeed c for a selection of vibration wavenumbers. It is seen that there is a pleasingly good correlation between the two, especially at the higher values of c where the theory ought to be more accurate. Despite our expectation that good agreement ought to be seen at large c, it is noteworthy that the results are still seem to be of some validity even down to quite modest values. An alternative view of the results is provided in Fig. 6. Here we consider the propulsion effect as a function of α for three values of c used in the earlier Fig. 2. Now we notice that the two-term result (55) is in excellent agreement with the computated results over much of the wavenumber range. In particular the locations and sizes of the local maximum and minimum in the propulsion effect are especially well captured.
We observe that the impressive performance of the asymptotic predictions deteriorate as α → 0 but this can be explained by consideration of changes to the solution structure that emerge in this limit. Then the key balances in the bulk region are modified; it can be verified that the convective streamwise term no longer drives the cross-slot viscous contributions once α ∼ O(ε −2 ). At this point the edge zones have grown to such an extent that they fill the whole slot and are indistinguishable from the bulk flow. This limit corresponds to the situation when we have 'sloshing' waves which are characterised by fluid moving backwards and forwards with a small net forward motion as evinced by the form of the mean flow term (18) when the constant B 0 defined by (52) is small. A typical sloshing motion is shown in Fig. 7. In the opposite case when α 1 it is evident from solution (11) that the ξ -dependent part of the flow becomes squashed into a O(α −1 ) zone located towards the bottom of the bulk region. Above this all that remains is a mean flow; to that extent it can be concluded that when α is large the wave-like part of the motion, though restricted to a small region close to the corrugated wall, is sufficiently powerful so as to modify the mean flow across the entire slot. This is typical Fig. 6 The size of U top /A 2 for three of the wavespeeds used in Fig. 2. The numerical results are shown as solid lines whilst the two-term asymptotic prediction (55) is designated by the dashed curves of a 'moving-wall' type of response in which the fluid moves forward everywhere except close to the lower plate and is powerful enough to significantly alter the propulsion effect felt by the upper plate; an example of this type of motion is presented in Fig. 8. Once α ∼ O(ε −1 ) the wave activity resides completely within the wall layer. The wall layer equations then necessitate separate numerical solution and analytical progress is somewhat limited.
In this work we have developed an expression for the speed the upper wall. While this is an interesting quantity there are other characteristics of the problem that could equally well be considered such as the forms of the velocity fields or the flux through the channel. In the interests of brevity we do not write out the various expressions here, but they can be deduced routinely by using the solutions outlined above. The system we have explored here is a model of an alternative propulsion system in which part of the energy inputted is used to mitigate the resistance rather than trying to overcome it. Here we have concentrated on using small-scale vibrations when the phase speed of the vibrations is relatively large and have managed to explain many of the phenonmena predicted by the computations summarised in [3]. Finally, we point out that related calculations described in [7] have looked at problems in which vibrations are imposed on both walls of a channel. The methods described in this work provide the basis for the discussion of the large phase-speed properties that apply to that more complicated configuration. Fig. 8 A "moving-wall' motion corresponding to the parameter choices c = 1600, α = 10 and A = 0.01. Shown is the streamwise velocity profiles at four equally spaced location at x = nλ/4 where λ = 2π/α is the wavelength of the vibrations and n = 0, 1, 2 and 3. We remark on how the ξ -dependent motion is concentrated next to the lower plate Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

Conflict of interest
The authors have no competing interests to declare that are relevant to the content of this article. The referees are thanked for their insightful comments which led to several improvements in the paper.