2.1. Governing equations

Consider the plane motion of a circular cylindrical wing of radius $a$, performing simultaneously translational and rotational harmonic oscillations (heaving and pitching) with the same frequency $\omega$ and amplitudes $A$, $\varTheta$, respectively, (see figure 1) in a viscous incompressible fluid. To describe the fluid flow, let us choose a polar coordinate system $(r,\theta )$ moving translationally with the centre of the cylinder,

(2.1a,b)\begin{equation} x = x_O(t) +r \cos(\theta), \quad y = y_O(t) +r \sin(\theta), \end{equation}

where $(x, y)$ are coordinates in a fixed inertial Cartesian coordinate system and functions $x_O(t) , y_O(t)$ determine the law of translational motion of the centre of the cylinder. We also normalize the spatial variables, time and velocity by $a$, $\omega ^{-1}$, $U = A \omega$, respectively. Using as the main flow characteristics the vorticity $w$ and the stream function $\psi$ associated with the velocities $v_r$, $v_ \theta$ of the fluid motion by the relations

(2.2a–c)\begin{equation} v_r=\frac{1}{r} \frac{\partial \psi}{\partial \theta},\quad v_\theta={-}\frac{\partial \psi}{\partial r},\quad w ={-}\frac{1}{r}\left(\frac{\partial (r v_\theta)}{\partial r}-\frac{\partial v_r}{\partial \theta}\right), \end{equation}

we write the governing equations in the form

(2.3)$$\begin{gather} \frac{\partial w}{\partial t} -\frac{\gamma}{r} \frac{\partial(\psi, w)}{\partial(r,\theta)}=\beta^{{-}1} \varDelta w, \quad w =\varDelta \psi , \end{gather}$$

(2.4)$$\begin{gather}\frac{\partial(\psi, w)}{\partial(r,\theta)}=\frac{\partial \psi}{\partial r} \frac{\partial w}{\partial \theta}-\frac{\partial \psi}{\partial \theta} \frac{\partial w}{\partial r}, \end{gather}$$

where $\varDelta$ is the Laplace operator, $\beta$ and $\gamma$ are the dimensionless frequency and the amplitude of translational oscillations, respectively. These dimensionless complexes are defined as follows:

(2.5a,b)\begin{equation} \beta=\frac{a^2\omega}{\nu}, \quad \gamma=A/a. \end{equation}

Hereinafter, $\gamma$ will be considered as a small parameter.

Figure 1.

Scheme of synchronous rotational–translational oscillations of the cylinder. The dashed line shows the rightmost position of the cylinder and the solid line indicates neutral position. Here $A$ and $\varTheta$ are the amplitudes of translational and rotational oscillations.

We supplement the equations of fluid motion with the following boundary conditions:

(2.6)$$\begin{gather} r = 1: \quad \psi = 0, \quad \frac{{\partial \psi }}{{\partial r}} = \xi \cos{t}, \end{gather}$$

(2.7)$$\begin{gather}r \to \infty : \quad \psi \sim{-}r\sin{\theta} \cos{t} + u^*_{\infty} r \cos{(\theta+\alpha)}. \end{gather}$$

On the cylinder surface ($r=1$) the no-slip conditions (2.6) are set. Here $\xi$ is the third dimensionless parameter of the problem. It defines the ratio of the amplitudes of rotational ($\varTheta$) and translational ($A$) oscillations,

(2.8)\begin{equation} \xi=a\varTheta/A=\varTheta/\gamma. \end{equation}

The condition (2.7) at infinity consists of two terms. The first one specifies a homogeneous oscillatory flow in the direction of cylinder oscillations (parallel to the $Ox$ axis). The second term specifies a stationary flow directed at an angle of $\alpha$ to the $Oy$ axis, it determines the average velocity of the cylinder in the fluid. Further, we will consider the case of a cruising motion. Under conditions of small oscillation amplitudes, the order of cruising speed, obviously, cannot be higher than $O(\xi \gamma )$, therefore, we will assume that $u^* _{\infty } = u_ {\infty } \xi \gamma$. The values of $u_{\infty }$ and $\alpha$ will be considered as unknown and determined from the condition that the time-average forces acting on the cylinder are equal to zero.

2.2. Force characteristics

Let us define the hydrodynamic forces $f_x, f_y$ and the moment $m$ acting per unit length of the cylinder. Normalizing the forces to $f_{sc} = 2 {\rm \pi}\rho a^2U \omega$ (where $\rho$ is the density of the fluid), the moment to $a f_{sc}$, we represent them in terms of vorticity as (e.g. Dynnikova & Andronov Reference Dynnikova and Andronov2018)

(2.9) \begin{equation} \left. \begin{array}{c@{}} \displaystyle f_y-if_x =\dfrac{i}{2}\sin{t}+\dfrac{1}{2{\rm \pi}\beta}\int_0^{2{\rm \pi}}\left(\dfrac{\partial w}{\partial r}-w\right)_{r=1}\, {\rm e}^{{\rm i}\theta} \,{\rm d} \theta,\\ \displaystyle m =\dfrac{2 \xi}{\beta}\cos{t}-\dfrac{1}{2{\rm \pi}\beta}\int_0^{2{\rm \pi}}(w)_{r=1} \,{\rm d} \theta. \end{array}\right\} \end{equation}

The first term in (2.9a) is the Krylov–Froude force, which appears as a result of the transition to a moving coordinate system. The terms proportional to $w$ and $\partial {w} / \partial r$ represent the contributions of the shear stresses and pressure, respectively.

The forces and moments are used to determine the time-average energy dissipation. Normalized to $D{sc} = U f_{sc}$, it can be written as

(2.10)\begin{equation} \left. \begin{array}{c@{}} d =d_x+\xi^2 d_\theta,\\ d_x={-}\langle f_x(t) \cos{t}\rangle, \quad d_\theta=\xi^{{-}1}\langle m(t) \cos{t}\rangle. \end{array}\right\} \end{equation}

Angle brackets hereinafter denote averaging over the oscillation period. The sign in the expression for $d_\theta$ is determined by the fact that the angular velocity of rotation of the cylinder is $-\xi \cos {t}$.

3.1. Expansion in the small parameter

Without making any assumptions about the parameter $\beta$, let us assume that $\gamma \ll 1$ and $\xi \sim 1$. We will construct a solution of the problem in the form of a series in powers of $\gamma$, limiting ourselves to the first two terms,

(3.1a,b)\begin{equation} \psi = {\psi ^{(0)}} + \gamma {\psi ^{(1)}} + \cdots ,\quad w = {w^{(0)}} + \gamma {w^{(1)}} + \cdots . \end{equation}

In the leading term, we obtain the Stokes approximation of the original problem:

(3.2) \begin{equation} \left. \begin{array}{c@{}} \dfrac{\partial }{{\partial t}}{w^{(0)}} - {\beta ^{ - 1}}\varDelta {w^{(0)}} = 0,\quad {w^{(0)}} = \varDelta {\psi ^{(0)}};\\ r = 1:\quad \psi^{(0)} = 0, \quad \dfrac{{\partial \psi^{(0)} }}{{\partial r}} = \xi \cos{t}, \quad r \to \infty : \quad \psi^{(0)} \sim{-}\cos{t}\cdot r\sin{\theta}. \end{array}\right\} \end{equation}

The Stokes solution of the problem describes a purely periodic motion of a fluid.

For the next term we have

(3.3) \begin{equation} \left. \begin{array}{c@{}} \dfrac{\partial }{{\partial t}}{w^{(1)}} - {\beta ^{ - 1}}\Delta {w^{(1)}} = {G},\quad {w^{(1)}} = \Delta {\psi ^{(1)}},\quad {G} = \dfrac{1}{r} \dfrac{\partial(\psi^{(0)}, w^{(0)})}{\partial(r,\theta)}; \\ r = 1:\quad \psi^{(1)} =\dfrac{{\partial \psi^{(1)} }}{{\partial r}} = 0, \quad\quad r \to \infty : \quad \psi^{(1)} \sim \xi u_\infty r \cos{(\theta+\alpha)}. \end{array} \right\} \end{equation}

The solution of the second subproblem contains both periodic and stationary components. The latter is often called secondary flow or steady streaming. As will be shown below, steady streaming, in the considered case of rotational–translational oscillations, provides the directed motion of the cylinder in the fluid.

3.2. Stokes approximation

In the Stokes approximation, the linear problem (3.2) splits into two. The first of them is the problem of an oscillatory fluid flow around a fixed cylinder, the second is the problem of rotational oscillations of a cylinder in an undisturbed fluid. Accordingly, the total flow is represented as the sum of these two flows

(3.4)\begin{equation} \left. \begin{array}{c@{}} \psi^{(0)}= \sin{\theta} \mathrm{Re} ( \varPsi_{1}(r) \, {\rm e}^{{\rm i}t}) +\xi \mathrm{Re} ( \varPsi_{2}(r) \, {\rm e}^{{\rm i}t}),\\ w^{(0)}=\sin{\theta}\, \mathrm{Re} ( W_{1}(r) \, {\rm e}^{{\rm i}t})+\xi \mathrm{Re} ( W_{2}(r) \, {\rm e}^{{\rm i}t}). \end{array}\right\} \end{equation}

The solutions of the corresponding one-dimensional problems are

(3.5)\begin{equation} \left. \begin{array}{c@{}} {\varPsi _1} ={-}r +\dfrac{1}{r}+\dfrac{S_x}{r} + {z^{ - 2}}{W_1},\quad {W_1} ={-}2z \dfrac{{{K_1}({z r})}}{{{K_0}({z})}}, \quad S_x(z) = \dfrac{{2{K_1}({z})}}{{{z}{K_0}({z})}},\\ {\varPsi _2} = S_\theta + { z^{ {-}2}}{W_2},\quad {W_2} ={-}z \dfrac{{{K_0}({z r})}}{{{K_1}(z)}},\quad S_\theta(z) = \dfrac{{{K_0}({z })}}{{{z}{K_1}({z})}}, \quad z =\sqrt{i \beta}. \end{array}\right\} \end{equation}

Here $K_j$ are modified Bessel functions. As can be seen, outside the Stokes boundary layer with a thickness of the order of $\beta ^{-1/2}$, the vorticity (its $W_1$ and $W_2$ components) decays exponentially and the flow becomes potential. Additional viscous terms $S_x r^{- 1}$ and $S_\theta$, which disappear in an effectively inviscid fluid (at $\beta = \infty$), correct the non-circulation flow of an ideal fluid outside the boundary layer.

Based on the obtained solution, it is possible to determine the hydrodynamic forces and moment acting on the cylinder and the time-average energy dissipation (in the first approximation):

(3.6)\begin{equation} \left. \begin{array}{c@{}} f^{(0)}_x= \text{Im} \left[{\rm e}^{{\rm i}t}\left(\dfrac{1}{2}+S_x\right)\right], \quad m^{(0)}=\dfrac{2 \xi}{\beta}\cos{t}-\xi \text{Im} [{\rm e}^{{\rm i}t}(S_\theta)],\\ \quad d^{(0)}_x={-}\dfrac{1}{2} \text{Im} [S_x], \quad d^{(0)}_\theta=\dfrac{1}{\beta}-\dfrac{1}{2} \text{Im} [S_\theta]. \end{array} \right\} \end{equation}

For $\beta = \infty$ (effectively inviscid fluid) functions, $S_x = S_\theta =0$ vanish, therefore the moment and force acting on the cylinder are purely inertial. They are proportional to the angular and linear acceleration of the cylinder with the added mass coefficients equal to $1/2$. For $\beta <\infty$, additional force components appear. They are completely determined by the functions $S_x(z)$ and $S_\theta (z)$.

In the practically important case of large $\beta$, the behaviour of $S_x(z)$ and $S_\theta (z)$ is described by the asymptotics

(3.7)\begin{equation} z \to \infty: \begin{array}{c} S_\theta(z) \sim \dfrac{1}{z}-\dfrac{1}{2z^2}+\dfrac{3}{8z^3}-\dfrac{3}{8z^4}+O(z^{{-}5}),\\ S_x(z) \sim \dfrac{2}{z}+\dfrac{1}{z^2}-\dfrac{1}{4z^3}+\dfrac{1}{4z^4}+O(z^{{-}5}). \end{array} \end{equation}

Substituting (3.7) into (3.6) we obtain the following expansion of functions $d_x$, $d_\theta$ for large values of $\beta$:

(3.8)\begin{equation} \beta \to \infty: \begin{array}{c} d^{(0)}_{\theta }=\dfrac{\beta^{{-}1/2}}{2 \sqrt{2} }+\dfrac{3 \beta^{{-}1}}{4}+\dfrac{3\beta^{{-}3/2}}{16 \sqrt{2} }+O(\beta^{{-}5/2}),\\ d^{(0)}_x=\dfrac{\beta^{{-}1/2}}{\sqrt{2}}+\dfrac{\beta^{{-}1}}{2 } -\dfrac{\beta^{{-}3/2}}{8 \sqrt{2}}+O( \beta^{{-}5/2}). \end{array} \end{equation}

3.3. Secondary flow problem

Let us define the first-order terms in the asymptotic representation (3.1a,b). We will be interested only in their stationary components, which describe the steady streaming. The same components should be kept on the right-hand side $G$ of (3.3). Separating them, we get

(3.9)\begin{equation} \langle G \rangle =G_x\sin{2\theta}+\xi G_\theta \cos\theta, \end{equation}

where

(3.10a,b)\begin{equation} G_x= \frac{1}{4r} \mathrm{Re} ( \bar{W}_1\varPsi'_1-\bar{\varPsi}_1W'_1), \quad G_\theta=\frac{1}{2r} \mathrm{Re} (\bar{W}_1\varPsi'_2 -\bar{\varPsi}_1W'_2). \end{equation}

Here and below, the overbar denotes the complex conjugation, and the prime mark denotes the derivative with respect to $r$. Carrying out direct calculations and using the notation

(3.11a,b)\begin{equation} \varPhi_k(r;\beta)=\frac{K_k(\sqrt{i\beta}r)}{K_k(\sqrt{i\beta})}, \quad S(\beta)=\frac{K_2(\sqrt{i\beta})}{K_0(\sqrt{i\beta})}, \end{equation}

we get

(3.12)$$\begin{gather} G_x(r;\beta) = \frac{\beta}{2} {\text{Im}}\left(\frac{\bar{S}\varPhi_0}{r^2}+2S\varPhi_2 \overline{\varPhi}_0-S\varPhi_2 \right), \end{gather}$$

(3.13)$$\begin{gather}G_\theta (r;\beta) = \frac{\beta}{r}{\text{Im}} (S)|\varPhi_1|^2 -\frac{\beta}{2}{\text{Im}} ( ( 1- \bar{S}r^{{-}2}) \varPhi_1). \end{gather}$$

Finding the solution of (3.3) in the form

(3.14)\begin{equation} {\psi^{(1)}} = \varPsi_x(r)\sin{2\theta}+\xi\varPsi_\theta(r)\cos\theta, \end{equation}

we obtain the following problems for $\varPsi _x$ and $\varPsi _\theta$:

(3.15a–c)$$\begin{gather} {\mathcal{L}}_2^2{\varPsi_x}={-}\beta {G_x}, \quad {\mathcal{L}}_1^2{\varPsi_\theta}={-}\beta {G_\theta}, \quad {\mathcal{L}}_k= \frac{1}{r}\frac{{\rm d}}{{\rm d} r}r\frac{{\rm d}}{{\rm d} r} - \frac{k^2}{r^2}, \end{gather}$$

(3.16a,b)$$\begin{gather}r = 1:\ {\varPsi_x} ={\varPsi'_x}= {\varPsi_\theta} ={\varPsi'_\theta}=0, \quad r \to \infty :\ \varPsi'_x=0 , \quad \varPsi'_\theta = u_\infty. \end{gather}$$

The right-hand sides of (3.15a,b) decrease exponentially with increasing $r$. As a consequence, there are unique solutions $\tilde {\varPsi }_x$, $\tilde {\varPsi }_\theta$ of the inhomogeneous (3.15a,b) that exponentially decrease as r goes to infinity. Using the method of variation of an arbitrary constant (details are given in Appendix A), they can be represented in the following integral form:

(3.17a,b)\begin{align} \tilde{\varPsi}_x={-}\frac{\beta r^4}{48}\int_1^\infty G_x(rx)( x^2-1 )^3 \frac{{\rm d}\kern0.06em x}{x}, \quad \tilde{\varPsi}_\theta={-}\frac{\beta r^4}{16}\int_1^\infty G_\theta(rx)( x^4-4x^2\ln{x}-1)\,{{\rm d}\kern0.06em x}. \end{align}

Adding to $\tilde {\varPsi }_x$, $\tilde {\varPsi }_\theta$ general solutions of homogeneous equations, we get

(3.18a,b)\begin{align} \varPsi_x= \tilde{\varPsi}_x-\frac{B_x}{r^2}+C_x\left( 1-\frac{1}{r^2}\right),\quad \varPsi_\theta= \tilde{\varPsi}_\theta-\frac{B_\theta}{r}+C_\theta\left( r-\frac{1}{r}\right) +D\left( r\ln{r}-\frac{r}{2} + \frac{1}{2r}\right). \end{align}

The boundary conditions for $r = 1$ are satisfied by choosing arbitrary constants $B_{x, \theta }$, $C_{x, \theta }$:

(3.19a–c)\begin{equation} B_{x,\theta}=\tilde{\varPsi}_{x,\theta}(1), \quad C_\theta={-}\frac{\tilde{\varPsi}'_\theta(1)+\tilde{\varPsi}_\theta (1)}{2},\quad C_x={-}\frac{\tilde{\varPsi}'_x(1)+2\tilde{\varPsi}_x(1)}{2}. \end{equation}

The function $\varPsi _x$ goes to the constant $C_x$ as $r \to \infty$. This does not contradict the boundary conditions at infinity. The function ${\varPsi }_\theta$ has the classical logarithmic singularity at infinity. It is impossible to determine the constant $D$ and satisfy the condition for the stream function at $r \to \infty$ for an arbitrary velocity of a stationary flow without taking into account the convective terms (Stokes’ paradox). Only with a special choice of $u_\infty$ and $\alpha$,

(3.20a,b)\begin{equation} u_\infty=C_\theta(\beta),\quad \alpha=0, \end{equation}

is the condition at infinity satisfied, and the constant $D$ vanishes. It turns out that this choice exactly corresponds to the cruising motion condition, that the force averaged over the period, acting on the body from the fluid, is equal to zero.

This can be shown by calculating, according to (2.9), the contributions to the force $f_y$ from the four terms of ${\varPsi }_\theta$ in (3.18b). The second and third do not generate vorticity, and therefore do not make any contribution. The contribution from the first term, up to a factor, is determined by the integral (Appendix A)

(3.21)\begin{equation} \int_1^\infty G_\theta (x;\beta) x^2 {{\rm d}\kern0.06em x} =0, \end{equation}

the explicit calculation of which also gives zero. The contribution of the last term is $D$. To fulfil the cruise flight condition we must set $D = 0$. Therefore, the last term in (3.18b) should be discarded. In this case, the cylindrical wing will move along the $Oy$ axis with a cruising speed $\xi \gamma C_\theta$.

The change of the cruising speed coefficient $u_\infty = C_\theta$ depending on dimensionless frequency $\beta$ is shown in figure 2. The obtained results show the presence of two types of cruising motion of the cylinder: forward ($u_\infty >0$) and reverse ($u_\infty <0$). The first of them is realized for $\beta > \beta _\theta = 13.5355$ and is characterized by a monotonic growth of the cruising speed with increasing frequency $\beta$. The maximum speed of such a motion is $3/_4\gamma \xi U=3/_4\varTheta U$, which is comparable to the amplitude of the velocity $U$ of translational oscillations of the cylinder. The reverse motion occurs when $\beta <\beta _\theta$. The cruising speed coefficient reaches the smallest value $u_\infty = -0.0534$ at $\beta = 2.18$. This value is more than an order of magnitude less than the maximum value of the cruising speed coefficient, $u_\infty = 0.75$, of direct motion.

Figure 2.

Dependence of the main integral characteristics on the dimensionless frequency ${\beta }$. Dashed lines represent asymptotics of different orders.

For the case $\beta \to \infty$, $u_\infty (\beta )$ is determined by the following asymptotic relationship:

(3.22)\begin{equation} u_\infty = C_\theta = \frac{3}{4}-7\sqrt{2}\beta^{{-}1/2}+\frac{1117}{16}\beta^{{-}1}- \frac{3183}{8\sqrt{2}}\beta^{{-}3/2}+O( \beta^{{-}2}). \end{equation}

The asymptotic behaviour of other main integral characteristics of the steady streaming at high oscillation frequencies is determined by the formulae

(3.23)$$\begin{gather} C_x = \frac{3}{4}-\frac{603}{32\sqrt{2}}\beta^{{-}1/2}+\frac{3681}{32} \beta^{{-}1}-\frac{405159}{512\sqrt{2}}\beta^{{-}3/2}+O( \beta^{{-}2}), \end{gather}$$

(3.24)$$\begin{gather}B_x = \frac{13}{2\sqrt{2}}\beta^{{-}1/2}-\frac{137}{4}\beta^{{-}1} +\frac{3963}{16\sqrt{2}}\beta^{{-}3/2}+O( \beta^{{-}2}), \end{gather}$$

(3.25)$$\begin{gather}B_\theta = \frac{11}{2\sqrt{2}}\beta^{{-}1/2}-\frac{123}{4}\beta^{{-}1} +\frac{2957}{16\sqrt{2}}\beta^{{-}3/2}+O( \beta^{{-}2}). \end{gather}$$

The asymptotics for $\beta \to \infty$ are shown by dashed lines in figure 2 together with the exact solution (solid lines). Table 1 also shows the boundaries of their range of applicability, which were estimated as the values of $\beta$ at which the relative approximation error reaches 5 %.

Table 1.

The limits of applicability of asymptotics of various order of $n$ with respect to the parameter $\beta$. The values (of $\beta$) in the table are given for the case when the error of the asymptotics reaches 5 %.

5.1. The limits of applicability of the asymptotic theory

To investigate the limits of applicability of the asymptotic theory, we carried out a numerical simulation of the flow past a flapping cylindrical wing. A numerical solution was constructed for the full problem (2.3)–(2.7) using the finite volume method. For this purpose, the numerical scheme developed in Nuriev, Kamalutdinov & Egorov (Reference Nuriev, Kamalutdinov and Egorov2019); Nuriev et al. (Reference Nuriev, Egorov and Kamalutdinov2021) was used. The details of the simulations are described in Appendix D.

In figure 5 a graph of the asymptotic dependence of cruising speed coefficient on frequency and the numerical simulation data obtained in the range of parameters $150 <\beta <6500$, $0.031 <\gamma <0.8$, $0.4 \leq \xi \leq 15$ are presented. As one can see, the values of the velocity coefficient obtained numerically for the same frequency $\beta$ change depending on the values of $\gamma$ and $\varTheta = \xi \gamma$. Recall here that in the framework of the asymptotic expansion, only the main stationary term was found, which determines the cruising speed as a linear function of $\varTheta$. The next terms of the asymptotic expansion define corrections proportional to higher powers of $\varTheta$ and $\gamma$ (and their combinations), which start to play a significant role as the amplitudes of rotational and translational oscillations increase (for this reason, for relatively large values of $\varTheta$ and $\gamma$, there is no single universal similarity parameter). The deviation of the numerical data from the linear (in relation to $\varTheta$) asymptotic approximation just indicates the magnitude of nonlinear corrections for different values of $\varTheta$ and $\gamma$.

Figure 5.

Dependence of the cruising speed coefficient $u_\infty$ on the dimensionless frequency $\beta$. The asymptotic solution is shown by the solid line, the numerical results are indicated by markers.

Nevertheless, the dependence of cruising speed coefficient on frequency predicted by the main stationary term of asymptotic model is well traced in the entire designated range of parameters. At moderate frequencies $150 \leq \beta \leq 470$, most of the numerical results have minimal deviations ($\pm 12\,\%$) from the theoretical dependence. In the high-frequency region $\beta > 470$, the asymptotics determines the upper values of the velocity coefficient (which are achieved at low values of $\gamma$).

The deviations $err = |1-u^{num}_\infty /u_\infty | \cdot 100\,\%$ of numerical values of the speed coefficient $u^{num} _ \infty$ from the theoretical curve $u_\infty (\beta )$ for different values of $\gamma$ and $\varTheta$ are shown in the bubble plots in figure 6 (each marker in figure 5 corresponds to a bubble in figure 6). The size of the bubble indicates the magnitude of deviation value (a large size corresponds to a large deviation), the coordinates show the values of $\gamma$ and $\varTheta$, the colour indicates the value of $\beta$. This representation of data allows one to identify parametric regions with minimal deviations between analytical and numerical results, which can be considered as the areas of applicability of the theory.

Figure 6.

Deviations $err = |1-u^{num}_\infty / u_\infty |\times 100$ of numerical values of the cruising speed from the theoretical ones depending on the parameters $\gamma, \varTheta, \beta$. The size of the marker (bubble) indicates the magnitude of deviation value, the colour indicates the value of $\beta$ (see the colour map).

In the diagrams (figure 6), blue rectangles indicate areas where the deviations $err$ did not exceed 12 % for all values of $\beta$ in the indicated ranges. As you can see, for moderate-frequency oscillations, the asymptotic model has a very wide range of applicability: $0<\gamma \leq 0.8$, $0.3 \leq \varTheta < 1.7$. Here the theory is valid for cases when the amplitude of translational oscillations does not exceed $4/5$ of the radius, and the angle of rotation is not more than $98^{\circ }$ and not less than $17^{\circ }$. The lower constraints on the parameter $\varTheta$ are obviously related to the specificity of the interaction between translational and rotational oscillations at small angles of rotation. The theory can be applied in the range $\varTheta < 0.3$ for relatively large values of $\xi$, however, detailed calculations to identify the exact dependences in this narrow range have not been performed. In the region of high values of the dimensionless frequency, the limits of applicability of the theory in terms of the amplitudes $\gamma$ significantly decrease: here the deviation $err<12\,\%$ is observed for $\gamma \leq 0.3$ (see figure 6).

The accuracy of the cruising speed determination (in the areas of applicability) is undoubtedly related to the accuracy with which the asymptotic theory in the obtained approximation can describe the structure of the flow. Figure 7 shows comparisons between the flow patterns of secondary streaming obtained analytically (right-hand side of the image) and numerically (left-hand side of the image). The visualization is presented for $\beta = 157$, where the structural features of the flows are most clearly distinguishable (all RC are of a sufficiently large size). According to the asymptotic theory (see § 4), the flow structure depends only on the parameter $\xi$. At a given frequency, this parameter varied in the range $0.46 \leq \xi \leq 11$. The values of the second parameter $\gamma$ in the calculations were selected so as to be within and on the boundaries of the selected applicability area.

Figure 7.

Flow patterns of secondary streaming obtained analytically (right-hand black and white part of the image) and numerically (left-hand half of the image) for $\beta =157$ and $\gamma =0.64$, $\xi = 0.46$ (a), $\gamma =0.64$, $\xi = 0.72$ (b), $\gamma =0.64$, $\xi = 2$, (c), $\gamma =0.13$, $\xi = 4$, (d), $\gamma =0.16$, $\xi = 7.5$ (e), $\gamma =0.16$, $\xi = 11$ ( f). Arrows show the direction of propulsive motion. The values of $err$ indicate the deviation of the numerical values of the cruising speed from the theoretical ones.

In the largest part of the computational domain $0.72 \leq \xi \leq 7.5$, regime D was localized; its boundaries are in good agreement with the asymptotic theory (see the flow regime map). As you can see in figure 7(c,d), in this regime the asymptotic model well predicts the changes in the size and position of RC in the flow region depending on the parameter $\xi$. This is also reflected in the estimates of the cruising speed: at the interior points $\xi = 2, 4$, the values of $err$ (see also figure 6) are minimal and equal to $err = 5.57, 1.92$, respectively. Note that the calculation for $\xi = 2$ was performed at $\gamma = 0.64$, i.e. when a peak-to-peak amplitude of translational oscillations is 1.28 times greater than the cylinder radius.

Changes of the flow structure near the lower boundary of regime D are shown in figure 7(a,b) for $\xi =0.46, 0.72$ and $\gamma =0.64$. According to the flow regime map, these points are located in the transition zones between regimes D and E. The features of the flow transformations in the wake of the body, observed in numerical calculations, are somewhat different from the structural changes predicted by the asymptotic theory. As can be seen in figure 7(a) at $\xi = 0.46$ in regime E, four external RC should be observed in the wake, in fact (according to the results of numerical calculations) only two cells of a larger size are formed. According to the classification proposed in § 4, the observed regime is of type F, not E. When passing from the case in figure 7(b) to the case figure 7(c), $err$ values increase from 11.7 to 16.9, thus, we reach the lower limit of applicability of the theory for the parameter $\varTheta$.

In figure 7(e, f) the changes in the flow occurring at large values of $\xi = 7.5, 11$ are shown. According to the theory, in this parametric range, the flow is transformed from regime D to regime A: a pair of RC located in the front part of the cylinder for smaller $\xi$ enclose the entire cylinder. As can be seen, the numerical results confirm such structural changes. However, the sizes of the attached cells predicted by the numerical and analytical models are different. This is directly reflected in the values of $err$, which begin to increase in regime A. This is how the upper bound of the region of applicability of the theory in the parameter $\varTheta$ is reached.

Let us consider the changes in the flow structure when approaching the upper limit of applicability in terms of the oscillation amplitude $\gamma$, for this we analyse the vorticity field. In figure 8(a) the instantaneous flow pattern in regime D is shown (the dynamics of the flow is presented in the supplementary movies which are available at https://doi.org/10.1017/jfm.2022.287). As you can see, the vortex wake formed behind the body dissipates very quickly. Away from the body, (which is confirmed by the analytical solution) an ideal fluid flow is observed. In the vicinity of the boundary of the applicability range (figure 8b,c), vortex streets are formed behind the body, they decay more slowly than the vortex wake in regime D. Depending on $\beta$, the cruising speed in these regimes differs by 10 %–20 % from the theoretical one. In the case in figure 8(b) the speed is lower than the theoretical one, and in the case in figure 8(c) it is higher. A vortex street in the form in figure 8(c) was observed only in the range of moderate values of $\beta$. With a further increase of the amplitude $\gamma$, the flow loses symmetry about the $Oy$ axis (this type of stability loss of a symmetric flow behind a flapping wing was discussed in detail by Dynnikova et al. (Reference Dynnikova, Dynnikov, Guvernyuk and Malakhova2021)). As a result, an additional stationary component of hydrodynamic force arises, acting parallel to the $Ox$ axis. The cruising speed in this mode in the direction of $Oy$ is no more than $2/3$ of the theoretical.

Figure 8.

Vorticity field. Instantaneous flow patterns for (a) $\beta =785$, $\varTheta =1.75, \gamma =0.16$; (b) $\beta =785$, $\varTheta =1.9, \gamma =0.48$; (c) $\beta =393$, $\varTheta =1.59, \gamma =0.8$; (d) $\beta =785$, $\varTheta =1.9, \gamma =0.64$. Arrows show the direction of propulsive motion.

5.2. Propulsive efficiency

As already shown above, the considered propulsive mechanism makes it possible to achieve a relatively high cruising speed; at the same time, the question of the energy efficiency of such a motion remains unclear. We define the efficiency $\varepsilon$ as the ratio of the power required for the uniform motion of a body in a fluid (under steady towing conditions) at a constant speed $V$ to the power required by the propulsor to move the body with the same cruising speed $V$. The power, in these cases, is exactly equal to the energy dissipation by the fluid.

When the cylinder moves at a constant speed $V$, the energy dissipation can be calculated as

(5.1)\begin{equation} D_0=C_{D}(Re)\rho a V^3.\end{equation}

Here $C_D$ is the drag coefficient depending on the Reynolds number

(5.2)\begin{equation} Re =\frac{2aV}{\nu}. \end{equation}

For moderate Reynolds numbers, $C_D$ for a round cylinder is of the order of unity (e.g. Wieselsberger Reference Wieselsberger1921).

For the considered flapping cylindrical wing, the (dimensional) energy dissipation can be defined as

(5.3)\begin{equation} D=2{\rm \pi}\rho U^2a^2\omega (d_x+\xi^2 d_\theta)=2{\rm \pi} \rho U^2\nu \beta \gamma^{{-}2} (\gamma^2d_x+\varTheta^2 d_\theta). \end{equation}

Taking into account that $V=U\varTheta u_\infty$ and $Re =2\gamma \varTheta \beta u_\infty$, it can be represented as

(5.4)\begin{equation} D=2{\rm \pi}\rho a( \varTheta^2 d_\theta +\gamma^{2}d_x) \left(\frac{V}{\gamma\varTheta u_\infty}\right)^3. \end{equation}

Using (5.1) and (5.4), the efficiency can be calculated by the formula

(5.5)\begin{equation} \varepsilon=\frac{D_0}{D}= \frac{C_D}{2{\rm \pi}} \frac{\gamma^3\varTheta^3 u_\infty^3}{\varTheta^2 d_\theta +\gamma^{2}d_x}. \end{equation}

For the first quantitative estimates, we can use the previously obtained asymptotic dependences (3.8), (3.22) for $d_\theta, d_x, u_\infty$. In this case, the dissipation and speed coefficients depend only on $\beta$. This, in turn, means that at a fixed frequency, the highest efficiency must be reached at the maximum allowable values of $\varTheta$ and $\gamma$.

As shown earlier, in the range of moderate frequencies, the upper limit of applicability is $\gamma = 0.8$, $\varTheta = 1.7$, i.e. according to the formula (5.5), an efficiency exceeding 20 % can be achieved in this case for $\beta > 314$. At high frequencies, the maximum values of the amplitudes are much lower $\gamma = 0.3, \varTheta = 1.7$; accordingly, comparable theoretical efficiency can be achieved for $\beta > 6000$.

To refine these estimates, we use the numerical results presented in the previous section of the work. Data for some typical parameter values are shown in table 3, where for each simulated case, the numerical values of the damping coefficients $d_x, d_\theta$, the numerical and theoretical values of the cruising speed coefficient $u_\infty ^{num}, u_\infty (\beta )$, the numerical and theoretical estimates of the efficiency $\varepsilon _ {num}, \varepsilon _{a}$, are presented.

Table 3.

Analytical and numerical estimates of the key characteristics of the propulsive motion. The column with the parameters of the most efficient motion is marked in bold.

Since the efficiency of movement significantly depends on the cruising speed (as $u_\infty ^3$), the deviation of this value from the theoretical one most significantly affects the estimates of $\varepsilon$. This, in particular, leads to the fact that at high oscillation frequencies, where the cruising speed coefficient (obtained numerically) $u_\infty ^{num} < u_\infty (\beta )$ rapidly decreases near the upper boundary of applicability, the theory significantly overestimates the efficiency (see cases $\beta = 3142, 6283$). In the range of small amplitudes in the vicinity of the upper limit of applicability, the cruising speed is higher than the theoretical one $u_\infty ^{num} \geq u_\infty (\beta )$, which should lead to some increase in efficiency. However, the underestimation of the speed is compensated by the underestimation of the value of $d_x$ in the denominator, which is caused by the low accuracy of the approximation of $d_x \approx d^{(0)}_x$ at large amplitudes. Therefore, here the theoretical and numerical estimates of $\varepsilon _ {num}, \varepsilon _ {a}$ coincide quite well.

The maximum efficiency of motion $\varepsilon \approx 26\,\%$, achieved at $Re \approx 420$ (the corresponding flow regime is shown in figure 8c), is a relatively high score for this range of Reynolds numbers. This result can be compared with the efficiency estimates (calculated in a similar way) of fish-like movement (Maertens et al. Reference Maertens, Triantafyllou and Yue2015) and with the efficiency of various low-Reynolds-number swimmers (Hancock & Newman Reference Hancock and Newman1953; Purcell Reference Purcell1977; Becker, Koehler & Stone Reference Becker, Koehler and Stone2003). The maximum efficiency of the carangiform swimmer at $2500 < Re <5000$ is approximately $40\,\%$, and of the anguilliform swimmer is approximately $30\,\%$. Low-Reynolds-number swimmers have significantly lower efficiency: the Purcell's swimmer has a peak efficiency of $0.77\,\%$, the undulating rod is approximately $7.4\,\%$ and the helical swimmer is approximately $8.6\,\%$. Thus, in terms of relative energy consumption, a cylindrical flapping wing can be attributed to the most efficient propulsors in the range of $Re \sim 10^2\text{--} 10^3$.