1 Introduction
Numerous undulatory swimming animals produce thrust with their caudal fins, while the anterior part of their body remains almost rigid (Lauder & Tytell Reference Lauder and Tytell2005). Motivated by applications to artificial propellers, these caudal fins have inspired studies on thrust generation by oscillating rigid foils (e.g. Freymouth Reference Freymouth1988; Koochesfahani Reference Koochesfahani1989; Schouveiler, Hover & Triantafyllou Reference Schouveiler, Hover and Triantafyllou2005; Godoy-Diana, Aider & Wesfreid Reference Godoy-Diana, Aider and Wesfreid2008) and oscillating flexible foils (e.g. Alben et al. Reference Alben, Witt, Baker, Anderson and Lauder2012; Moored et al. Reference Moored, Dewey, Smits and Haj-Hariri2012; Dewey et al. Reference Dewey, Boschitsch, Moored, Stone and Smits2013; Quinn, Lauder & Smits Reference Quinn, Lauder and Smits2014; Paraz, Schouveiler & Eloy Reference Paraz, Schouveiler and Eloy2016). In these studies, the thrust production and the propulsion efficiency are often linked to the characteristics of the wake, with the idea that most of the dynamical information is contained in the fish ‘footprint’ (Müller et al. Reference Müller, Van Den Heuvel, Stamhuis and Videler1997).
The question of whether the oscillating frequency of a foil, or in the context of fish swimming, the tailbeat frequency, is somewhat related to the wake instability frequency has first been addressed by Triantafyllou, Triantafyllou & Grosenbaugh (Reference Triantafyllou, Triantafyllou and Grosenbaugh1993). In their seminal paper, they considered the stability of the average experimental flow behind a pitching foil, as reported by Koochesfahani (Reference Koochesfahani1989). The principal wake parameter is the Strouhal number,
$\mathit{St}=fA/U$
, with
$f$
the oscillating frequency,
$U$
the flow velocity and
$A$
the width of the wake, taken to be the peak-to-peak amplitude of the foil trailing edge. From linear stability analyses, they found that these jet wake profiles are convectively unstable and are likely to behave as noise amplifiers when excited close to the resonance frequency with a maximum amplification for
$0.25<\mathit{St}<0.35$
(see also Huerre & Monkewitz Reference Huerre and Monkewitz1990).
Studying experimentally rigid foils activated in pitch and heave, Triantafyllou et al. (Reference Triantafyllou, Triantafyllou and Grosenbaugh1993), and later Schouveiler et al. (Reference Schouveiler, Hover and Triantafyllou2005), showed that the Froude efficiency (i.e. the ratio between thrust power and input power) reaches a maximum within the same range of Strouhal number. They argued that swimming performance is intimately linked to the characteristics of the wake instability. Similar correlations between the flapping frequency and the frequency associated with the maximum amplification of the jet wake have also been reported in numerical simulations (e.g. Lewin & Haj-Hariri Reference Lewin and Haj-Hariri2003). It has also been claimed that swimming animals benefit from this efficiency peak by swimming within the same range of Strouhal number:
$0.25<\mathit{St}<0.35$
(Triantafyllou et al.
Reference Triantafyllou, Triantafyllou and Grosenbaugh1993; Taylor, Nudds & Thomas Reference Taylor, Nudds and Thomas2003), although this argument is still debated today (Eloy Reference Eloy2012; van Leeuwen, Voesenek & Müller Reference van Leeuwen, Voesenek and Müller2015). More recently this so-called ‘wake resonance theory’ has been re-examined by Moored et al. (Reference Moored, Dewey, Smits and Haj-Hariri2012, Reference Moored, Dewey, Boschitsch, Smits and Haj-Hariri2014), who performed stability analyses on averaged experimental profiles using a locally parallel flow assumption. They concluded that multiple local maximal efficiencies can be reached, each corresponding to a frequency of a different unstable mode of the averaged wake profile. However, despite correlations between these two frequencies, the causal link between the stability properties of the averaged wake profile and the swimming performances has so far proved difficult to establish.
The wake resonance theory is based on a simplifying assumption: both the stability analyses and the experimental studies only consider the wake generated by a thrust-producing foil and neglect the influence of the rest of the body. Yet, when a self-propelled body swims at constant speed, thrust and drag balance on average, and the wake behind such a body is momentumless.
Considering that the features of wakes behind swimming bodies can be reduced to a combination of spatially localised forces acting on the fluid, a family of momentumless wakes has been proposed by Afanasyev (Reference Afanasyev2004). In this work, the flow induced by non-translating localised forces (described, for instance, by Cantwell Reference Cantwell1986) is generalised to the case of forces translating at constant velocity. When a doublet of such translating forces is considered, a family of momentumless wakes parametrised by the intensity of the force doublet, the swimming velocity and the Reynolds number is obtained (figure 1). The objective of the present work is to address the stability of such momentumless wakes.
Illustration of the wake produced by a force doublet
$Q$
placed in a uniform flow velocity
$U$
.
The paper is organised as follows. In § 2, the analytical description and the dimensional analysis of the wake profile is given. The locally parallel stability characteristics of the wake profiles are then computed in § 3. The profiles are shown to undergo a transition from a convective to absolute instability in the Reynolds number–force doublet intensity space. The non-parallelism of the wake base flow and the nonlinear behaviour of the unstable modes are addressed in § 4 by performing a global stability analysis through time stepping of the Navier–Stokes system. In § 5, we discuss whether the present approach may apply to swimming animals by considering empirical data on animals with different lengths. We also discuss the connection between the momentumless wake profiles and jet profiles. Finally, some tentative conclusions are given in § 6.
2 Family of momentumless wakes
A self-propelled body moving at constant velocity experiences a drag equal and opposite to the thrust produced. Both thrust and drag trigger a momentum transfer between the self-propelled body and the fluid. In the fluid this momentum is evidenced by the emission of vortex dipoles in two dimensions and vortex rings in three dimensions. In the far field, these two opposite forces can be reduced to a translating force doublet in a fixed coordinate framework as being described by Afanasyev (Reference Afanasyev2004).
In the Stokes approximation, the two-dimensional streamfunction of a single impulsive force located at the origin of the coordinate system
$(x^{\ast },y^{\ast })$
(noting dimensional quantities with
$^{\ast }$
) is given by (Cantwell Reference Cantwell1986; Afanasyev Reference Afanasyev2004)
$$\begin{eqnarray}\unicode[STIX]{x1D713}_{I^{\ast }}^{\ast }(I^{\ast },x^{\ast },y^{\ast },t^{\ast })=\frac{I^{\ast }y^{\ast }}{2\unicode[STIX]{x03C0}({x^{\ast }}^{2}+{y^{\ast }}^{2})}(1-\exp (-\unicode[STIX]{x1D709}^{2})),\end{eqnarray}$$
with
$$\begin{eqnarray}\unicode[STIX]{x1D709}=\sqrt{\frac{{x^{\ast }}^{2}+{y^{\ast }}^{2}}{4\,\unicode[STIX]{x1D708}^{\ast }t^{\ast }}},\end{eqnarray}$$
$I^{\ast }$
the impulsive force intensity
$([I^{\ast }]=L^{3}T^{-1})$
and
$\unicode[STIX]{x1D708}^{\ast }$
the kinematic viscosity. The streamfunction of two opposite forces separated by distance
$\unicode[STIX]{x1D716}^{\ast }$
is
$$\begin{eqnarray}\unicode[STIX]{x1D713}^{\ast }(x^{\ast },y^{\ast },t^{\ast })=\unicode[STIX]{x1D713}_{I^{\ast }}^{\ast }(x^{\ast }+\unicode[STIX]{x1D716}^{\ast }/2,y^{\ast },t^{\ast })-\unicode[STIX]{x1D713}_{I^{\ast }}^{\ast }(x^{\ast }-\unicode[STIX]{x1D716}^{\ast }/2,y^{\ast },t^{\ast }),\end{eqnarray}$$
which, in the limit
$\unicode[STIX]{x1D716}^{\ast }\rightarrow 0$
, reduces to
$$\begin{eqnarray}\unicode[STIX]{x1D713}_{M^{\ast }}^{\ast }(M^{\ast },x^{\ast },y^{\ast },t^{\ast })=\unicode[STIX]{x1D716}^{\ast }\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}_{I^{\ast }}^{\ast }}{\unicode[STIX]{x2202}x^{\ast }}=\frac{M^{\ast }x^{\ast }y^{\ast }}{\unicode[STIX]{x03C0}({x^{\ast }}^{2}+{y^{\ast }}^{2})^{2}}(1-(1+\unicode[STIX]{x1D709}^{2})\exp (-\unicode[STIX]{x1D709}^{2})),\end{eqnarray}$$
with
$M^{\ast }([M^{\ast }]=L^{4}T^{-1})$
the doublet intensity
$$\begin{eqnarray}M^{\ast }=\lim _{\unicode[STIX]{x1D716}^{\ast }\rightarrow 0,I\rightarrow \infty }I^{\ast }\unicode[STIX]{x1D716}^{\ast }.\end{eqnarray}$$
Let a force, either a single force or a force doublet, move with a constant speed
$U^{\ast }$
in the negative
$x^{\ast }$
direction. By solving the diffusion–advection equation corresponding to the Oseen approximation, we can obtain the streamfunction of a moving force with intensity
$J^{\ast }$
(
$[J^{\ast }]=L^{3}T^{-2}$
) for a single force or
$Q^{\ast }$
(
$[Q^{\ast }]=L^{4}T^{-2}$
) for a force doublet
$$\begin{eqnarray}\unicode[STIX]{x1D713}_{J^{\ast }\,\text{or}\;Q^{\ast }}^{\ast }(J^{\ast }\text{ or }Q^{\ast },x^{\ast },y^{\ast },t^{\ast })=\int _{0}^{t^{\ast }}\unicode[STIX]{x1D713}_{I^{\ast }\,or\;M^{\ast }}^{\ast }(J^{\ast }\text{ or }Q^{\ast },x^{\ast }-U^{\ast }(t^{\ast }-\unicode[STIX]{x1D70F}^{\ast }),y^{\ast },t^{\ast }-\unicode[STIX]{x1D70F}^{\ast })\,\text{d}\unicode[STIX]{x1D70F}^{\ast }.\end{eqnarray}$$
Taking the derivative of (2.6) with respect to
$y^{\ast }$
, we obtain the expression for the steady streamwise velocity
$u^{\ast }(x^{\ast },y^{\ast },t^{\ast })$
.
To make this problem dimensionless, we choose the constant swimming speed
$U^{\ast }$
as the reference velocity, and
$$\begin{eqnarray}\unicode[STIX]{x1D6FF}^{\ast }=\sqrt{\frac{\unicode[STIX]{x1D708}^{\ast }x^{\ast }}{U^{\ast }}},\end{eqnarray}$$
as the reference length, which is analogous to a boundary-layer reference length. Noting dimensionless variables without stars, we have
$$\begin{eqnarray}x=\frac{x^{\ast }}{\unicode[STIX]{x1D6FF}^{\ast }}=\mathit{Re}=\frac{U^{\ast }\unicode[STIX]{x1D6FF}^{\ast }}{\unicode[STIX]{x1D708}^{\ast }}=\sqrt{\frac{U^{\ast }x^{\ast }}{\unicode[STIX]{x1D708}^{\ast }}},\end{eqnarray}$$
$\mathit{Re}$
being the Reynolds number.
According to the expression of the streamfunction (2.6), for sufficiently large integration time, we obtain a steady-state solution with dimensionless streamwise velocity profiles for the single force and force doublet given by (in the framework attached to the translating body)
$$\begin{eqnarray}\displaystyle & \displaystyle u_{J}(\mathit{Re},y)=1+\frac{J}{2\unicode[STIX]{x03C0}\mathit{Re}}\int _{0}^{\infty }\unicode[STIX]{x1D6F7}_{J}[\mathit{Re}-(t-\unicode[STIX]{x1D70F}),y,t-\unicode[STIX]{x1D70F}]\,\text{d}\unicode[STIX]{x1D70F}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle u_{Q}(\mathit{Re},y)=1+\frac{Q}{\unicode[STIX]{x03C0}\mathit{Re}^{2}}\int _{0}^{\infty }\unicode[STIX]{x1D6F7}_{Q}[\mathit{Re}-(t-\unicode[STIX]{x1D70F}),y,t-\unicode[STIX]{x1D70F}]\,\text{d}\unicode[STIX]{x1D70F}, & \displaystyle\end{eqnarray}$$
with
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D6F7}_{J}(\mathit{Re},y,t)=\frac{\text{e}^{-\unicode[STIX]{x1D709}^{2}}}{(\mathit{Re}^{2}+y^{2})^{2}}((\text{e}^{\unicode[STIX]{x1D709}^{2}}-1)(\mathit{Re}^{2}+y^{2})+2y^{2}(1+\unicode[STIX]{x1D709}^{2}-\text{e}^{\unicode[STIX]{x1D709}^{2}})), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D6F7}_{Q}(\mathit{Re},y,t)=\frac{\mathit{Re}\,\text{e}^{-\unicode[STIX]{x1D709}^{2}}}{(\mathit{Re}^{2}+y^{2})^{3}}((\mathit{Re}^{2}-3y^{2})\text{e}^{\unicode[STIX]{x1D709}^{2}}-(\mathit{Re}^{2}-3y^{2})(1+\unicode[STIX]{x1D709}^{2})+2\unicode[STIX]{x1D709}^{4}y^{2}), & \displaystyle\end{eqnarray}$$
and
$$\begin{eqnarray}\unicode[STIX]{x1D709}=\frac{1}{2}\sqrt{\mathit{Re}}\,\sqrt{\frac{\mathit{Re }^{2}+y^{2}}{t}}.\end{eqnarray}$$
The dimensionless intensities of the single force and doublet are given by
$$\begin{eqnarray}J=\frac{J^{\ast }}{U^{\ast }\unicode[STIX]{x1D708}^{\ast }}\quad \text{and}\quad Q=\frac{Q^{\ast }}{{\unicode[STIX]{x1D708}^{\ast }}^{2}},\end{eqnarray}$$
these two quantities being connected via the doublet size
$\unicode[STIX]{x1D716}^{\ast }$
through
$$\begin{eqnarray}Q=\frac{\unicode[STIX]{x1D716}^{\ast }J^{\ast }}{{\unicode[STIX]{x1D708}^{\ast }}^{2}}=\unicode[STIX]{x1D716}\,\mathit{Re}\,J,\end{eqnarray}$$
where
$\unicode[STIX]{x1D716}=\unicode[STIX]{x1D716}^{\ast }/\unicode[STIX]{x1D6FF}^{\ast }$
.
The expression (2.10) describes a family of momentumless wakes, parametrised by the dimensionless force doublet intensity
$Q$
and the Reynolds number
$\mathit{Re}$
(
$\mathit{Re}$
can also be viewed as a dimensionless streamwise distance). In the same manner, equation (2.9) is a family of jet wakes parametrised by the force intensity
$J$
and
$\mathit{Re}$
. In § 3, we will study the linear stability of the momentumless wakes with a locally parallel flow assumption. In § 4, we will consider a non-parallel case in which we also need the lateral velocity component which is given by
$$\begin{eqnarray}\displaystyle & \displaystyle v_{Q}(\mathit{Re},y)=-\frac{Q}{\unicode[STIX]{x03C0}\mathit{Re}^{2}}\int _{0}^{\infty }\unicode[STIX]{x1D711}_{Q}[\mathit{Re}-(t-\unicode[STIX]{x1D70F}),y,t-\unicode[STIX]{x1D70F}]\,\text{d}\unicode[STIX]{x1D70F}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D711}_{Q}(\mathit{Re},y,t)=\frac{y\,\text{e}^{-\unicode[STIX]{x1D709}^{2}}}{(\mathit{Re}^{2}+y^{2})^{3}}(-(3\mathit{Re}^{2}-y^{2})\text{e}^{\unicode[STIX]{x1D709}^{2}}+(3\mathit{Re}^{2}-y^{2})(1+\unicode[STIX]{x1D709}^{2})+2\unicode[STIX]{x1D709}^{4}\mathit{Re}^{2}).\qquad & \displaystyle\end{eqnarray}$$
The linear stability of the jet–wake profile family (2.9) will be considered in § 5, in connection with the momentumless wake instability results.
3 Linear stability analysis
A locally parallel linear stability analysis is performed by numerically solving the linearised Navier–Stokes equations around the wake profile given by (2.10). Perturbations are expressed as normal modes, that is
$$\begin{eqnarray}\left(\begin{array}{@{}c@{}}u\\ v\\ p\end{array}\right)=\left(\begin{array}{@{}c@{}}\hat{u} (y)\\ \hat{v}(y)\\ \hat{p}(y)\end{array}\right)\text{e}^{\text{i}(\unicode[STIX]{x1D6FC}x-\unicode[STIX]{x1D714}t)},\end{eqnarray}$$
with
$\unicode[STIX]{x1D714}\in \mathbb{C}$
their frequency and
$\unicode[STIX]{x1D6FC}\in \mathbb{C}$
their wavenumber. For a given force doublet intensity
$Q$
and Reynolds number
$\mathit{Re}$
, there is an infinite number of frequencies solution of a Orr–Sommerfeld-type eigenvalue problem, each corresponding to a different eigenmode.
To solve this eigenvalue problem, a Chebyshev collocation discretisation has been used in the
$y$
-direction. Since we want to account for the infinite domain in the
$y$
-direction, and because the base flow has large variations in a narrow region around
$y=0$
(rendering the problem relatively stiff), we use the following mapping (Peyret Reference Peyret2002)
$$\begin{eqnarray}y=\frac{a\unicode[STIX]{x1D702}}{\sqrt{1+(a/H)^{2}-\unicode[STIX]{x1D702}^{2}}},\quad -H\leqslant y\leqslant H,\end{eqnarray}$$
with
$\unicode[STIX]{x1D702}\in [-1,1]$
the Chebyshev collocation points. The typical parameter value used for
$a$
is
$a=1.5$
and
$H$
has been chosen large enough for the perturbation to vanish at
$y=\pm H$
. We found that for
$H$
larger than
$30$
the unstable modes could be captured accurately and in order to guarantee the absence of finite height effects we used
$H=100$
. Given the stiffness of the problem, a thorough convergence study has been performed and
$1000$
discretisation points have been considered, ensuring that the eigenvalues
$\unicode[STIX]{x1D714}$
are converged almost to machine accuracy.
For each unstable eigenvalue (when
$\text{Im}(\unicode[STIX]{x1D714})>0$
), the question naturally arises whether the corresponding eigenmode is absolutely unstable or not. When the flow is absolutely unstable, the linear instability dynamics of the flow is of an oscillator type, in contrast to convectively unstable flows where initially triggered perturbation are convected downstream (Huerre & Monkewitz Reference Huerre and Monkewitz1990). To assess the convective/absolute character of the instability, we can use the cusp map method described in Schmid & Henningson (Reference Schmid and Henningson2000). This method consists in mapping the complex wavenumber plane (
$\unicode[STIX]{x1D6FC}$
-plane) into the complex frequency plane (
$\unicode[STIX]{x1D714}$
-plane), by computing the most unstable complex frequency for a given complex wavenumber (figure 2). A cusp is obtained when
$$\begin{eqnarray}\frac{\text{d}\unicode[STIX]{x1D714}_{0}}{\text{d}\unicode[STIX]{x1D6FC}}=0,\end{eqnarray}$$
with
$\unicode[STIX]{x1D6FC}\in \mathbb{C}$
. Absolute instability occurs when the imaginary part
$\unicode[STIX]{x1D714}_{0,i}>0$
and convective instability occurs when
$\unicode[STIX]{x1D714}_{0,i}<0$
. From a practical point of view, the computation is started by looking for the most unstable complex frequencies
$\unicode[STIX]{x1D714}$
for a real wavenumbers
$\unicode[STIX]{x1D6FC}$
. Then, the imaginary part of the complex wavenumber is gradually decreased, and for each value of this negative imaginary part
$\unicode[STIX]{x1D6FC}_{i}$
, we look for the most unstable frequencies, until we find a cusp as shown in figure 2 (corresponding here to an absolutely unstable case because the cusp appears for
$\unicode[STIX]{x1D714}_{0,i}>0$
). (Note that in the appendix A we compare and validate our procedure by considering a wake model proposed in the literature.)
Mapping from (a) the
$\unicode[STIX]{x1D6FC}$
-plane to (b) the
$\unicode[STIX]{x1D714}$
-plane for a wake profile with
$\mathit{Re}=1$
and
$Q=1000$
(each line corresponds to a different value of
$\unicode[STIX]{x1D6FC}_{i}$
). This mapping shows a cusp (marked by the circle) on the upper half of the
$\unicode[STIX]{x1D714}$
-plane, thus indicating an absolute instability.
To find the transition between absolute and convective instability, we look for wake parameters
$\mathit{Re}$
and
$Q$
, such that cusp appears for real frequency, i.e.
$\unicode[STIX]{x1D714}_{0,i}=0$
. By tracing the neutral curve and the absolute–convective transition curve together in the
$(\mathit{Re},Q)$
-plane, we obtain the instability map shown in figure 3(a). The
$(\mathit{Re},Q)$
-plane is divided into three regions: a stable region, a convectively unstable region and an absolutely unstable region. As an example, two profiles corresponding respectively to an absolutely unstable and a convectively unstable wake are shown in figure 3(c).
(a) Instability map in the
$\mathit{Re}$
–
$Q$
plane for a momentumless wake. AU, CU and S stand respectively for absolutely unstable, convectively unstable and stable. (b) Dimensionless frequency
$\unicode[STIX]{x1D714}_{0,r}$
at the convective–absolute transition. (c) Absolutely unstable profile for
$Q=10^{5}$
,
$\mathit{Re}=10$
(solid line), and convectively unstable profile for
$Q=10^{5}$
,
$\mathit{Re}=30$
(dashed line). Parameters of the two profiles are marked by coloured rectangles in (a).
For larger Reynolds numbers
$\mathit{Re}$
(i.e.
$\mathit{Re}\gtrsim 20$
), it is found that the transition from absolute to convective instability occurs for
$Q\propto \mathit{Re}^{3}$
, while the neutral curve satisfies
$Q\propto \mathit{Re}^{2}$
. In this framework of the locally parallel flow assumption, we also obtain the dimensionless frequency at the convective–absolute transition, which seems to converge asymptotically to
$0.93$
as both the doublet intensity and the Reynolds number increase as shown in figure 3(b).
4 Global wake dynamics
4.1 Global linear stability
Given the non-parallelism of the wakes described by (2.6), a question naturally arises: are the stability predictions based on the locally parallel flow assumption reliable? To address this question, we perform a global linear stability analysis of the non-parallel base flow.
For some inlet position
$x_{0}^{\ast }$
, we define
$\mathit{Re}_{0}$
as the Reynolds number formed with the reference length
$\unicode[STIX]{x1D6FF}_{0}^{\ast }=\sqrt{\unicode[STIX]{x1D708}^{\ast }x_{0}^{\ast }/U^{\ast }}$
at inflow. We define
$X$
as the dimensionless (using
$\unicode[STIX]{x1D6FF}_{0}^{\ast }$
) distance from inflow and the corresponding local Reynolds number may be written as
$$\begin{eqnarray}\mathit{Re}=\mathit{Re}_{0}\sqrt{1+\frac{X}{\mathit{Re}_{0}}}.\end{eqnarray}$$
The non-parallel evolution of the base flow may be taken into account by simply using (4.1) in the base-flow (2.10) formula. The transverse coordinate
$Y$
made dimensionless with
$\unicode[STIX]{x1D6FF}_{0}^{\ast }$
is then
$$\begin{eqnarray}Y=y\sqrt{1+\frac{X}{\mathit{Re}_{0}}}.\end{eqnarray}$$
The Navier–Stokes equations linearised around the non-parallel base flow
$\boldsymbol{U}(X,Y)$
become
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}t}\boldsymbol{u}=-(\boldsymbol{U}\boldsymbol{\cdot }\unicode[STIX]{x1D735})\boldsymbol{u}-(\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D735})\boldsymbol{U}-\unicode[STIX]{x1D735}p+\frac{1}{\mathit{Re}_{0}}\unicode[STIX]{x1D6FB}^{2}\boldsymbol{u}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}=0, & \displaystyle\end{eqnarray}$$
where the streamwise and transverse components of
$\boldsymbol{U}$
are given by (2.10) and (2.16) respectively.
The domain we consider is shown by a straight dashed line in figure 4(a) where I, II, III (corresponding to
$\mathit{Re}_{0}=20,21.83,22.4$
respectively) are the 3 different upstream inflow boundaries which have been used to solve the linearised Navier–Stokes system, for the force doublet intensity
$Q_{0}=1.9\times 10^{5}$
. The inflow position
$\text{I}$
has been chosen to be inside the absolutely unstable parameter region, whereas
$\text{III}$
is inside the convectively unstable domain,
$\text{II}$
being approximately on the absolute–convective transition boundary. In all the computations performed the outflow boundary corresponds to the point
$\text{IV}$
on the dashed line. The streamwise velocity profiles at the positions
$\text{I}$
,
$\text{III}$
and
$\text{IV}$
are shown in figure 4(b), which illustrates the important non-parallelism. Also, the profiles exhibit strong variations in
$Y$
which makes necessary to use a high discretisation when aiming at solving the linearised Navier–Stokes system (4.3) and (4.4).
(a) Set of non-parallel profiles at
$Q_{0}=1.9\times 10^{5}$
on the
$Q$
–
$\mathit{Re}$
-plane that illustrates our system. I is the inlet at
$\mathit{Re}_{0}=20$
while IV is the outlet at
$\mathit{Re}=36.87$
. (b) Streamwise variation of the base flow used in the direct numerical simulation at three different locations marked on (a).
In the context of large-scale global stability problems, time stepping approaches have become increasingly popular during the last decade. A comprehensive presentation of those matrix-free methods for open-flow stability problems can be found in Bagheri et al. (Reference Bagheri, Åkervik, Brandt and Henningson2009). Writing formally the linearised Navier–Stokes system as
$\unicode[STIX]{x2202}\boldsymbol{q}/\unicode[STIX]{x2202}t=\boldsymbol{A}\boldsymbol{q}$
, the solution at time
$T$
for any initial condition
$\boldsymbol{q}_{0}$
can formally be written as
$\boldsymbol{q}(T)=\text{e}^{\boldsymbol{A}T}\boldsymbol{q}_{0}$
, whatever method of time integration for the linearised Navier–Stokes system is used. Here, we are interested in the most amplified (or least damped) eigenvalue and we apply a Rayleigh iteration procedure by iterative time stepping, computing
$\boldsymbol{q}^{(k+1)}=\text{e}^{\boldsymbol{A}T}\boldsymbol{q}^{(k)}$
, by successively integrating over the time interval
$T$
the linearised Navier–Stokes system.
For this purpose, the code documented in Marquillie & Ehrenstein (Reference Marquillie and Ehrenstein2002), which is based on Chebyshev collocation method in the transverse
$Y$
-direction and finite difference in the streamwise direction has been adapted to our problem. In order to obtain a fully resolved flow field, the transverse
$Y$
-direction is discretised with 600 points in the range
$-40\leqslant Y\leqslant 40$
, which extends sufficiently far from the region with significant variations of the base-flow profiles (cf. figure 4
b). Note that in the direct numerical simulation the coordinate system is now normalised by
$\unicode[STIX]{x1D6FF}_{0}^{\ast }$
. The streamwise direction is therefore discretised using
$\unicode[STIX]{x0394}X=\unicode[STIX]{x0394}x/\unicode[STIX]{x1D6FF}_{0}^{\ast }=0.02$
and a fourth-order finite-difference scheme is used. Since we used 2400 discretisation points in the streamwise direction for the largest domain from
$\text{I}$
to
$\text{IV}$
, the distance from I to IV in the direct numerical simulation coordinate system is
$X=48$
.
The initial (divergence free) flow field
$\boldsymbol{q}^{(0)}$
considered is a Gaussian function
$$\begin{eqnarray}\left(\begin{array}{@{}c@{}}u\\ v\end{array}\right)=A\left(\begin{array}{@{}c@{}}-(Y-Y_{0})\\ (X-X_{0})\unicode[STIX]{x1D70E}_{Y}^{2}/\unicode[STIX]{x1D70E}_{X}^{2}\end{array}\right)\text{exp}\left(-\frac{(X-X_{0})^{2}}{2\unicode[STIX]{x1D70E}_{X}^{2}}-\frac{(Y-Y_{0})^{2}}{2\unicode[STIX]{x1D70E}_{Y}^{2}}\right),\end{eqnarray}$$
centred at
$Y_{0}=0$
and located relatively close to the inlet (with
$X_{0}=8$
for the inlet at
$\text{I}$
) the other parameters being
$\unicode[STIX]{x1D70E}_{X}=0.5$
and
$\unicode[STIX]{x1D70E}_{Y}=1$
. At the outlet, the advective boundary condition
$$\begin{eqnarray}\frac{\unicode[STIX]{x2202}\boldsymbol{u}}{\unicode[STIX]{x2202}t}+U_{ad}\frac{\unicode[STIX]{x2202}\boldsymbol{u}}{\unicode[STIX]{x2202}X}=0,\end{eqnarray}$$
has been considered. The value of
$U_{ad}=0.12$
proved to be appropriate to let the perturbation leave the domain without reflection. Also, the computed frequency of oscillation of the global structure appeared to be fairly insensitive to the exact choice of
$U_{ad}$
(the values 0.12 or 0.06 for
$U_{ad}$
giving the same global eigenvalue results). A zero Dirichlet boundary condition proved to be convenient to avoid reflections at inflow for the inlet positions II and III, the corresponding velocity profiles being at the margin of absolute instability (position II) or convectively unstable (position III) and no upstream propagating perturbations were encountered. For the inlet I, inside the absolutely unstable region, upstream propagating perturbations are expected and indeed the use of a Dirichlet condition led to spurious reflections at inflow and ultimately divergence was encountered. An advective boundary condition (4.6) has therefore also been applied at inflow in this case, with a negative advective velocity. Given the weak absolute instability, a small (in absolute value) velocity proved to be suitable and
$U_{ad}=-10^{-3}$
has been chosen as about the smallest value such that no perturbation wave reflections were encountered at inflow. The time interval
$T$
for the successive flow snapshots has to be appropriately chosen such that it satisfies the Nyquist criterion (cf. Bagheri et al.
Reference Bagheri, Åkervik, Brandt and Henningson2009), i.e. there should be at least two sampling points in one period of oscillation. The value
$T=1.1$
has been considered, which is small enough given that the period at the absolute–convective transition is roughly 7.7.
In this time stepping Rayleigh procedure, the value
$$\begin{eqnarray}\unicode[STIX]{x1D706}=\text{e}^{\unicode[STIX]{x1D70E}T},\end{eqnarray}$$
with the largest module is recovered after convergence, which provides the global leading stability eigenvalue
$\unicode[STIX]{x1D70E}=(\log \unicode[STIX]{x1D706})/T$
. Note that the time exponential is written
$\text{e}^{\unicode[STIX]{x1D70E}t}$
(and not
$\text{e}^{-\text{i}\unicode[STIX]{x1D714}t}$
as in the local analysis). To obtain both the instability amplification rate and the frequency, a variant of the simple Rayleigh iteration has been considered, by computing the coefficients
$\unicode[STIX]{x1D6FE}_{0,k}$
and
$\unicode[STIX]{x1D6FE}_{1,k}$
for every three steps in the procedure such that
$$\begin{eqnarray}\frac{1}{||\boldsymbol{q}^{(k)}||}(\unicode[STIX]{x1D6FE}_{0,k}\boldsymbol{q}_{j}^{(k)}+\unicode[STIX]{x1D6FE}_{1,k}\boldsymbol{q}_{j}^{(k+1)}+\boldsymbol{q}_{j}^{(k+2)})=0\end{eqnarray}$$
(by selecting two components
$\boldsymbol{q}_{j},j=j_{1},j_{2}$
of the vector fields to compute
$\unicode[STIX]{x1D6FE}_{0,k},\unicode[STIX]{x1D6FE}_{1,k}$
). Convergence implies that
$\unicode[STIX]{x1D6FE}_{0,k}\rightarrow \unicode[STIX]{x1D6FE}_{0}$
,
$\unicode[STIX]{x1D6FE}_{1,k}\rightarrow \unicode[STIX]{x1D6FE}_{1}$
, the complex conjugate pair of
$\unicode[STIX]{x1D706}$
being obtained by solving
$\unicode[STIX]{x1D6FE}_{0}+\unicode[STIX]{x1D6FE}_{1}\unicode[STIX]{x1D706}+\unicode[STIX]{x1D706}^{2}=0$
and the leading stability eigenvalue
$\unicode[STIX]{x1D70E}=\unicode[STIX]{x1D70E}_{r}\pm \text{i}\unicode[STIX]{x1D70E}_{i}$
is then calculated. This procedure has been applied for the different inlet locations and the global eigenvalues are given in table 1.
Most amplified global mode: real part of the streamwise velocity perturbation.
Global linear eigenvalue
$\unicode[STIX]{x1D70E}=\unicode[STIX]{x1D70E}_{r}\pm \unicode[STIX]{x1D70E}_{i}$
for three different inlets: inlet I (
$\mathit{Re}_{0}=20$
), inlet II (
$\mathit{Re}_{0}=21.83$
) and inlet III (
$\mathit{Re}_{0}=22.44$
).
For the inflow at
$\text{I}$
,
$\unicode[STIX]{x1D70E}_{r}$
is positive and a globally unstable mode is hence found. The real part of the eigenfunction’s streamwise velocity component is shown in the
$(X,Y)$
-plane in figure 5. The corresponding frequency,
$\unicode[STIX]{x1D70E}_{i}=0.63$
, is different from the frequency at the absolute–convective transition calculated by the local analysis of the previous section:
$\unicode[STIX]{x1D70E}_{i}\approx 0.8$
, which would be the expected instability frequency by the global mode if the flow were only weakly non-parallel (Chomaz Reference Chomaz2005). This illustrates the influence of the base flow’s non-parallelism in the present problem.
Interestingly, for the inlet
$\text{II}$
very near the location of the absolute–convective transition, the global mode is only weakly amplified and, for the inlet
$\text{III}$
in the convectively unstable region,
$\unicode[STIX]{x1D70E}_{r}<0$
and the leading global mode is damped. This provides evidence that the existence of a finite region of absolute instability is necessary in this highly non-parallel case for a global unstable mode to emerge.
4.2 Nonlinear disturbance evolution
To assess how the global instability evolves when considering the full nonlinear Navier–Stokes equations, the same numerical procedure as described above has been used, by adding the nonlinear term
$-(\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D735})\boldsymbol{u}$
to the right-hand side of (4.3). The same Gaussian-type initial condition (4.5) is considered, choosing an amplitude
$A=10^{-2}$
. The inlet corresponds to
$\mathit{Re}_{0}=20$
and the domain is reaching from
$\text{I}$
to
$\text{IV}$
(cf. figure 4
a). The discretisation is chosen to be the same as for the global stability analysis and a semi-implicit time marching is used.
Instantaneous field structures: (a) vorticity of the perturbation and (b) total vorticity at
$t=1140$
.
The instantaneous vorticity perturbation and the total vorticity in the fully nonlinear regime (at
$t=1140$
) is shown in figure 6. A reversed Bénard–von Kármán vortex street can be observed in the second half of the flow domain. From the time evolution of the perturbation, frequency spectra associated with the nonlinear dynamics can be calculated (figure 7). The spectra have been computed at two different locations slightly off the centreline:
$(X,Y)=(20,1.06)$
and
$(X,Y)=(40,1.06)$
. Figure 8 shows that the global nonlinear structure is tuned to a unique fundamental frequency of
$\unicode[STIX]{x1D70E}_{i}=0.646$
(the harmonics being visible as well), which is in agreement with the frequency of the linear global mode.
Evolution of vorticity perturbation in time recorded at (a)
$(X,Y)=(20,1.06)$
and (b)
$(X,Y)=(40,1.06)$
. Both signals exhibit the same fundamental frequency
$\unicode[STIX]{x1D70E}_{i}=0.646$
(figure 8).
Frequency spectra of the vorticity perturbation shown in figure 7 for the time range
$400\leqslant t\leqslant 1200$
.
5 Application to swimming animals
5.1 Absolute instability
Figure 3(a) shows that, when the force doublet intensity is large enough (
$Q\gtrsim 100$
), the near wake (corresponding to small
$\mathit{Re}$
) is always absolutely unstable. Moving further away from the doublet, there is first a transition from absolute to convective instability, and then to stability. The existence of an absolute region is important, because, in that case, the wake is expected to behave like an oscillator triggering a self-sustained instability process.
To examine the consequences of the present stability analyses on a swimming animal, we consider a self-propelled body of length
$L^{\ast }$
and moving at a constant speed
$U^{\ast }$
. In two dimensions, the skin friction drag (per unit length) exerted on this swimmer is of the order of
$$\begin{eqnarray}F^{\ast }\sim \unicode[STIX]{x1D70C}^{\ast }{U^{\ast }}^{2}L^{\ast }\mathit{Re}_{L^{\ast }}^{-1/2},\end{eqnarray}$$
which is nothing other than the Blasius boundary-layer law (Schlichting & Gersten Reference Schlichting and Gersten2003), where
$\mathit{Re}_{L^{\ast }}=U^{\ast }L^{\ast }/\unicode[STIX]{x1D708}^{\ast }$
is the Reynolds number based on body length.
For a constant swimming speed, this drag has to be balanced by an equal but opposite thrust. Assuming that the points of application of these two forces are separated by a distance of order
$L^{\ast }$
, the dimensional force doublet intensity can be estimated to be
$$\begin{eqnarray}Q^{\ast }=\frac{F^{\ast }L^{\ast }}{\unicode[STIX]{x1D70C}^{\ast }}\sim {U^{\ast }}^{3/2}{L^{\ast }}^{3/2}{\unicode[STIX]{x1D708}^{\ast }}^{1/2}.\end{eqnarray}$$
At the transition from absolute to convective instability, we found
$Q=18.4\,\mathit{Re}^{3}$
for
$\mathit{Re}>10$
. By recalling that
$\mathit{Re}=U^{\ast }\unicode[STIX]{x1D6FF}^{\ast }/\unicode[STIX]{x1D708}^{\ast }$
and by using the definition of
$Q$
in (2.14), we can infer the following relation
$$\begin{eqnarray}\unicode[STIX]{x1D6FF}_{ac}^{\ast }=\frac{0.38}{U^{\ast }}(Q^{\ast }\unicode[STIX]{x1D708}^{\ast })^{1/3},\end{eqnarray}$$
which states that the value of
$\unicode[STIX]{x1D6FF}^{\ast }$
corresponding to the absolute–convective transition (denoted by
$\unicode[STIX]{x1D6FF}_{ac}^{\ast }$
), can be estimated from
$Q^{\ast }$
,
$U^{\ast }$
and
$\unicode[STIX]{x1D708}^{\ast }$
. Following the same line of reasoning, the neutral stability curve corresponds to
$Q=35.22\,\mathit{Re}^{2}$
for
$\mathit{Re}>10$
(figure 3
a) and is associated with a critical
$\unicode[STIX]{x1D6FF}^{\ast }$
above which the wake is stable
$$\begin{eqnarray}\unicode[STIX]{x1D6FF}_{c}^{\ast }=\frac{0.17}{U^{\ast }}{Q^{\ast }}^{1/2}.\end{eqnarray}$$
Let us first consider a fish swimming in water (
$\unicode[STIX]{x1D708}^{\ast }=10^{-6}~\text{m}^{2}~\text{s}^{-1}$
) having body length
$L^{\ast }=10\,$
cm and moving with a velocity
$U^{\ast }=L^{\ast }~\text{s}^{-1}$
. Based on (5.2), we have
$Q^{\ast }\approx 10^{-6}~\text{m}^{4}~\text{s}^{-2}$
. By substituting these values into (5.3) and (5.4), we find
$\unicode[STIX]{x1D6FF}_{ac}^{\ast }=0.38\,$
mm and
$\unicode[STIX]{x1D6FF}_{c}^{\ast }=1.7\,$
mm. These values are far smaller than the width of the wake or the tailbeat amplitude, which is generally approximately one fifth of the body length (Videler Reference Videler1993). This means that, for a fish, the wake is likely to lie entirely in the stable region (figure 3
a).
Now consider the case of a swimming ascidian larva. The hydrodynamics of locomotion for these small swimmers has been addressed in McHenry, Azizi & Strother (Reference McHenry, Azizi and Strother2003), focusing in particular on the contribution of viscous and inertial force to the production of thrust and drag during steady undulatory swimming. To assess the relative importance of form drag and skin friction, the authors divide the larva into a spherical head and a rectangular tail. With these hypotheses, they find that only for small Reynolds numbers (
$Re_{L^{\ast }}<10$
) is the drag almost entirely due to skin friction, the contribution of form drag becoming increasingly important at higher Reynolds numbers. In McHenry et al. (Reference McHenry, Azizi and Strother2003), the model is validated by comparisons with measurements of thrust produced by ascidian larvae of body length
$L^{\ast }=1.9~\text{mm}$
, the value provided being approximately
$6\times 10^{-6}~\text{N}$
for a mean swimming speed of approximately
$31~\text{mm}~\text{s}^{-1}$
. The corresponding Reynolds number is
$Re_{L^{\ast }}=58.9$
and the two-dimensional skin-friction Blasius drag formula (5.1) used here would predict a drag force (per unit length) of
${\sim}2.4\times 10^{-4}~\text{N}\,\text{m}^{-1}$
. Although the extrapolation of a two-dimensional model to a three-dimensional body geometry is problematic, we may consider the analysis in Ehrenstein, Marquillie & Eloy (Reference Ehrenstein, Marquillie and Eloy2014), who addressed numerically the boundary layer of a periodically flapping plate with finite width in uniform incoming flow, assessing the longitudinal skin-friction force expression
$F_{3D}^{\ast }=C_{3D}\,\sqrt{r\langle |U_{\bot }|\rangle }\,\unicode[STIX]{x1D70C}^{\ast }{U^{\ast }}^{2}{L^{\ast }}^{2}\mathit{Re}_{L^{\ast }}^{-1/2}$
, with
$r$
the plate’s width to length ratio and
$\langle |U_{\bot }|\rangle$
the mean absolute value of the dimensionless (with the swimming speed) periodic wall-normal velocity, the coefficient
$C_{3D}$
being
${\approx}1.8$
. This expression is equivalent to (5.1), but these formulas, that aim at estimating drag for a swimming body, can only give crude values that would need to be assessed with precise measurements or numerical simulations. Assuming nevertheless this drag force law, we obtain a skin-friction drag force of
$C\,4.6\times 10^{-7}~\text{N}$
, which is to be compared to the measured thrust of
$6\times 10^{-6}~\text{N}$
, given the uncertainty of the value of the coefficient
$C$
.
Using the measured larva data in formula (5.2), we find a value
$Q^{\ast }\approx 4.5\times 10^{-10}~\text{m}^{4}~\text{s}^{-2}$
. Using (5.3), we then find
$\unicode[STIX]{x1D6FF}_{ac}^{\ast }\approx 0.094~\text{mm}$
and
$\unicode[STIX]{x1D6FF}_{c}^{\ast }=0.12~\text{mm}$
, which correspond to approximately 5 % and 6.3 % of the body length respectively. In this particular case, the wake is thus likely to be unstable and there may even be an absolute–convective transition in the near wake, provided that our model applies at least qualitatively to this swimming ascidian larva case. It is however unclear whether the instability properties of the wake will affect the swimmer performance since the instability cannot produce momentum.
5.2 Momentumless wake versus jet wake
Triantafyllou et al. (Reference Triantafyllou, Triantafyllou and Grosenbaugh1993) proposed that the swimming efficiency of a self-propelled body reaches a maximum when there is a resonance between the frequency of the wake instability and the tailbeat frequency. As it has already been noted in the introduction, the wake profiles considered by Triantafyllou et al. (Reference Triantafyllou, Triantafyllou and Grosenbaugh1993) are jets with a net positive momentum. In the present study, we consider momentumless wake profiles, but these wakes can still be decomposed into one part due to the thrust and one due to the drag. To compare our results with those of Triantafyllou et al. (Reference Triantafyllou, Triantafyllou and Grosenbaugh1993), we propose to extract the thrust part of our family of momentumless wakes and assess its stability properties.
We consider a swimming fish of body length
$L^{\ast }$
and wake half-width or tailbeat amplitude of approximately
$\unicode[STIX]{x1D6FF}^{\ast }\approx 0.1L^{\ast }$
(which is approximately the case for most undulatory swimmers). Thrust is generally produced by the caudal fin or by the posterior part of the body (Lighthill Reference Lighthill1969), while skin friction is expected to decrease along the length of the body as boundary-layer thickness increases. We can thus safely assume that the separation distance between the points of application of thrust and drag is in the interval
$0.1L^{\ast }<\unicode[STIX]{x1D716}^{\ast }<L^{\ast }$
, which means that the dimensionless doublet size is in the interval
$1<\unicode[STIX]{x1D716}<10$
. Now, using (2.15), we find that the force doublet intensity
$Q$
is connected to the thrust intensity
$J$
through the relation:
$J=Q/(\unicode[STIX]{x1D716}\mathit{Re})$
.
We now compare the stability properties of the momentumless profiles and jet profiles. To do so, we have performed a linear stability analysis of the jet wake profile of intensity
$J$
given by (2.9), using a method similar to that explained in § 3 for the momentumless wake profiles. The result of this analysis is shown in figure 9 together with the results of the momentumless wake. In these stability diagrams plotted in the
$(\mathit{Re},Q)$
or
$(\mathit{Re},J)$
-plane, for a given doublet intensity
$Q$
and Reynolds number
$\mathit{Re}$
, the corresponding jet thrust intensity (for a specific
$\unicode[STIX]{x1D716}$
) has to be chosen on the ordinate as
$J=Q/(\unicode[STIX]{x1D716}\mathit{Re})$
and the stability property may be inferred.
Stability diagram of a jet wake profile of intensity
$J$
(dashed line) and a momentumless wake profile of intensity
$Q$
(solid line).
Using again the example of a
$L^{\ast }=10\,$
cm fish swimming at constant speed
$U^{\ast }=L^{\ast }~\text{s}^{-1}$
, we have
$\mathit{Re}=U^{\ast }\unicode[STIX]{x1D6FF}^{\ast }/\unicode[STIX]{x1D708}^{\ast }\approx 10^{3}$
and
$Q=Q^{\ast }/({\unicode[STIX]{x1D708}^{\ast }}^{2})\approx 10^{6}$
. In that case, the momentumless wake is stable according to figure 9. Yet, its jet counterpart associated with the production of thrust only has an intensity
$10^{2}<J<10^{3}$
(with
$1<\unicode[STIX]{x1D716}<10$
), which corresponds to an unstable wake profile. The same holds for larger or faster fish. Hence, for most fish, the momentumless wake is stable while the jet profile due to the thrust alone is unstable.
To go further, we plot in figure 10 three jet profiles for different values of
$\unicode[STIX]{x1D716}$
together with their stability properties. As stated above, for these values of
$\unicode[STIX]{x1D716}$
, the profiles are unstable (
$\unicode[STIX]{x1D714}_{i}>0$
). Moreover, the (real) frequency
$\unicode[STIX]{x1D714}_{r}$
associated with the maximum of
$\unicode[STIX]{x1D714}_{i}$
is almost constant:
$\unicode[STIX]{x1D714}_{r}\approx 0.5$
. This frequency correspond to a Strouhal number based on the wake width
$A_{wake}$
,
$$\begin{eqnarray}\mathit{St}\approx \frac{A_{wake}\,\unicode[STIX]{x1D714}_{r}}{2\unicode[STIX]{x03C0}},\end{eqnarray}$$
where
$A_{wake}$
is estimated as the
$y$
-distance between the two inflection points of the thrust profiles. By substituting the value of
$A_{wake}$
(
$A_{wake}\approx 3.2$
for all values of
$\unicode[STIX]{x1D716}$
) and
$\unicode[STIX]{x1D714}_{r}$
into (5.5), the Strouhal number is found to be
$\mathit{St}\approx 0.25$
. This result is similar to the range of Strouhal number found by Triantafyllou et al. (Reference Triantafyllou, Triantafyllou and Grosenbaugh1993) from stability analyses of experimental jet profiles:
$0.25<\mathit{St}<0.35$
. However, as it has already been noted above, the corresponding momentumless profiles, found when both thrust and drag forces are taken into account, are stable (and therefore no Strouhal number can be defined in that case). The physical basis for this difference is however unclear. It may be due to the fact that velocity gradients are smaller in norm for the momentumless wakes since the velocity increase due to the thrust tends to be compensated by the velocity deficit due to the drag. Indeed, computing the maximum, in norm, of the velocity gradient for the three jet profiles shown in figure 10(a), we find approximately the values 0.12, 0.06 and 0.01, for respectively
$\unicode[STIX]{x1D716}=1,2$
and 10. The maximum values of the growth rates decrease monotonically with the velocity gradient norm, as seen in figure 10(b). The maximum of the velocity gradient for the corresponding stable momentumless profile, not shown here, is however much smaller (
$\max (|\text{d}u/\text{d}y|)\approx 1.5\times 10^{-4}$
).
(a) Jet wake profiles of extracted thrust for three different values of
$\unicode[STIX]{x1D716}$
and (b) their corresponding most unstable modes. The parameters used are
$\mathit{Re}=10^{3}$
,
$Q=10^{6}$
, and
$J=Q/(\unicode[STIX]{x1D716}\mathit{Re})$
.
5.3 Limit of validity
We will now examine the limit of validity of the present results. We will first discuss the validity of the Oseen approximation. We will then examine if the exchange of momentum between the swimming body and the flow can be described by a force doublet. Finally, we will discuss how three-dimensional effects may affect the results.
Evolution of the steady Oseen solution from
$t=0$
(solid line) to a new profile at
$t=20$
(dashed line) for
$Re=20.78$
.
In the present analysis, the base flow is obtained in the Stokes approximation and using the Oseen assumption. To assess the reliability of this approximation, the Oseen solution has been considered as the initial condition for the full Navier–Stokes system which has been integrated in time. The Oseen solution profile is held fixed at the inlet as a Dirichlet boundary condition and at outflow an advective boundary condition (similar to that for the flow perturbation in § 4) is considered. It is observed that inside the domain, the Oseen solution slightly evolves from
$t=0$
to
$t=20$
(see figure 11), but then the profile undergoes no significant change any more until
$t=165$
. These quasi-steady profiles slightly different from the steady Oseen solution have been considered for local stability computations, focusing in particular on position of absolute to convective instability transition. This new transition is found roughly at
$Re=20.8$
which is quite close to the value of
$Re=21.83$
(cf. figure 4
a) for the pure Oseen profile.
Pursuing the time integration of the Navier–Stokes system where the Oseen solution has been used as initial condition, the flow starts to oscillate (for
$t>165$
) and ultimately the flow dynamics reaches nonlinear saturation. The vorticity pattern is shown in figure 12 and the reverse von Kármán vortex street similar to figure 6(b) is observed. Frequency spectra of the flow are depicted in figure 13 which shows that the flow is tuned to a frequency of 0.66. This frequency is close to 0.646 found in the analysis of § 4 with the Oseen approximation as the base flow. The similarity of the results clearly lends credit that the Oseen assumption is justified, at least for the relatively low Reynolds number considered.
Instantaneous vorticity field structure at
$t=750$
.
Frequency spectra of the vorticity in time recorded at two stations off-centre: (a)
$(X,Y)=(30,-2.2)$
and (b)
$(X,Y)=(40,-2.2)$
, for the time range
$750\leqslant t\leqslant 1200$
.
The Oseen approximation has also been addressed in Gustafsson & Protas (Reference Gustafsson and Protas2012), who performed a detailed study on the solutions of the two-dimensional Oseen equations for the flow behind an obstacle for a broad range of Reynolds number (here, Reynolds number is based on the characteristic dimension of the obstacle). They compared their study with numerical simulations of Fornberg (Reference Fornberg1985) on a steady viscous flow around a two-dimensional cylinder. They concluded that the flow structures have a number of similarities with that of Oseen flow in terms of recirculation length, drag coefficient and separation angle. However, these similarities can only be observed up to a Reynolds number of approximately 100. For higher Reynolds number, the Oseen approximation is therefore debatable. In the present context, we make the assumption that the family of momentumless wakes still qualitatively describe the flow even for large Reynolds number. This assumption could however be assessed in the future by performing stability analyses of different families of momentumless wakes (e.g. profiles obtained from averaged nonlinear numerical simulations).
The family of momentumless wake profiles used in the present study is obtained by assuming that the contribution of thrust (respectively drag) forces can be reduced on average to a point force. Further, we assume that the flow field is that of a force doublet, which means that we consider the far field, i.e. distances large compared to the separation distance between the points of application of thrust and drag. Yet, as we saw above, most of the instability properties of the flow are related to the near wake, where these hypotheses do not hold. Nevertheless, the family of wake profiles used here captures qualitatively the principal features of momentumless wakes: no net momentum, profiles parametrised by their amplitude and width. As already proposed above, it would be interesting to extend the present investigation by studying the stability of average profiles obtained from numerical simulations of self-propelled bodies.
Finally, we may wonder how the present two-dimensional results may generalise to three dimensions. First, it should be noted that the flow behind swimming animals is generally ‘very three-dimensional’, or said differently, the wake in a given horizontal slice may appear to have a positive or negative momentum depending on its depth (Müller et al. Reference Müller, Van Den Heuvel, Stamhuis and Videler1997; Drucker & Lauder Reference Drucker and Lauder2000; Lauder & Drucker Reference Lauder and Drucker2002; Nauen & Lauder Reference Nauen and Lauder2002). The only exception seems to be the flow behind eels (Tytell & Lauder Reference Tytell and Lauder2004) and, of course, two-dimensional numerical simulations. An interesting avenue for future works would be to extend the stability analyses to three dimensions. But, at the present time, it is difficult to extrapolate our two-dimensional results to three dimensions. It is however safe to assume that the stability properties of a momentumless wake and its jet part will still be very different, even in three dimensions.
6 Conclusion
In this paper, we have performed a linear stability analysis of a family of momentumless wake profiles. As a base flow, we used the flow generated by a translating force doublet in the Oseen approximation, as initially proposed by Afanasyev (Reference Afanasyev2004). This flow takes into account the opposite drag and thrust that are exerted on average on a self-propelled body swimming at constant speed
$U^{\ast }$
. For this family of momentumless wakes, a transition from absolute to convective instability has been found, in contrast with the jet profiles usually considered for bio-inspired propeller wakes that are only convectively unstable. Within the locally parallel flow assumption, it was found that this transition occurs at
$$\begin{eqnarray}\unicode[STIX]{x1D6FF}_{ac}^{\ast }=\frac{0.38}{U^{\ast }}(Q^{\ast }\unicode[STIX]{x1D708}^{\ast })^{1/3},\end{eqnarray}$$
where
$\unicode[STIX]{x1D6FF}^{\ast }$
is a measure of the wake half-width and
$Q^{\ast }$
is the intensity of the force doublet. This value is too small to be meaningful for fish longer than approximately 1 cm. It may however be relevant for swimming organisms with low Reynolds number such as tadpoles or larvae. For a specific dimensionless intensity of the force doublet possibly in the range of centimetre size swimming organisms, the Oseen approximation used in the model being validated only for the corresponding low Reynolds number regimes, the linear as well as nonlinear wake dynamics has been assessed for the non-parallel base flow. The existence of a global instability behaviour tuned at a specific frequency is found, whenever the wake’s inflow boundary is chosen in the locally absolute instability region.
The thrust part may be extracted from the momentumless wake profiles. When it is intense enough, the corresponding jet profile is found to be (convectively) unstable, even though the momentumless wake from which it has been extracted may well be stable. Although further investigations on the Oseen approximation, different wake profiles, and three dimensions, may be required, the present analysis demonstrates that physical interpretations of swimming efficiency based on jet wakes have to be taken with great care and cannot be easily transposed to a whole self-propelled body. In other words, for swimming animals, the selection of Strouhal number through a ‘wake resonance’, as originally proposed by Triantafyllou et al. (Reference Triantafyllou, Triantafyllou and Grosenbaugh1993) and later revisited by Moored et al. (Reference Moored, Dewey, Boschitsch, Smits and Haj-Hariri2014), seems unlikely.
Acknowledgement
We are grateful to E. Lauga and T. Patel for their participation in the initial stage of this work.
Appendix A
In this appendix we briefly address, whether single force profiles of type (2.9) may capture conventional mean wake profile characteristics, by appropriately tuning the parameters. As an example, we consider the profile proposed by Monkewitz (Reference Monkewitz1988) and given by
$$\begin{eqnarray}\displaystyle & \displaystyle U(y)=1-\unicode[STIX]{x1D6EC}+2\unicode[STIX]{x1D6EC}F(y), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D6EC}=(U_{c}^{\ast }-U_{max}^{\ast })/(U_{c}^{\ast }+U_{max}^{\ast }), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle F(y)=(1+\text{sinh}^{2N}(y\;\text{sinh}^{-1}(1)))^{-1}. & \displaystyle\end{eqnarray}$$
For particular values of the parameters (
$\unicode[STIX]{x1D6EC}=-1.105$
and
$N=1.34$
, at
$Re=12.5$
) the stability characteristics of the Monkewitz’s profile have been recomputed, showing that the flow is absolutely unstable, the absolute frequency and wavenumber being
$$\begin{eqnarray}\unicode[STIX]{x1D714}_{0}=0.952+\text{i}0.058,\quad \unicode[STIX]{x1D6FC}_{0}=0.8-\text{i}0.505.\end{eqnarray}$$
These values are approximately equal to the values computed by Monkewitz (Reference Monkewitz1988) (see table I, page 1004), the small differences being probably due to a better convergence here.
Singlet force profile (dashed line) fitted to a generic bluff-body wake profile proposed by Monkewitz (Reference Monkewitz1988) (solid line).
A single force streamfunction given by (2.6) can be, by appropriately choosing the parameters, fitted to the generic profile proposed by Monkewitz (Reference Monkewitz1988) (figure 14), by following the same non-dimensionalisation as in this latter paper. We find that the best fit is obtained for
$x=5$
and
$J=J^{\ast }/(2\unicode[STIX]{x03C0}l^{\ast }{U^{\ast }}^{2})=-0.804$
(
$l^{\ast }$
being the reference length), when a constant velocity equal to
$2$
is added to the streamwise velocity.
We then performed the local stability analysis of this single force profile and found that the absolutely unstable mode corresponds to
$$\begin{eqnarray}\unicode[STIX]{x1D714}_{0}=0.956+\text{i}0.031,\quad \unicode[STIX]{x1D6FC}_{0}=0.8-\text{i}0.505.\end{eqnarray}$$
These values are quite close to what was found for the generic profile proposed by Monkewitz (Reference Monkewitz1988), the amplification rate being slightly lower. Treating the Reynolds number as an independent parameter, we found that the transition from absolute to convective instability is approximately at
$Re=9.75$
with
$$\begin{eqnarray}\unicode[STIX]{x1D714}_{0}=0.945,\quad \unicode[STIX]{x1D6FC}_{0}=0.8-\text{i}0.505.\end{eqnarray}$$
This value of Reynolds number for the transition from absolute to convective instability is similar to the result of Monkewitz (Reference Monkewitz1988) (see figure 2, page 1001).
From these linear stability analyses, we first conclude that our method is validated by the results of the literature. Second, these analyses show that the properties of the instability are not very sensitive to the profile family used. In other words, the single force model can qualitatively capture the instability properties of generic bluff-body wake, and similarly the doublet force model is likely to capture the instability of a momentumless wake.
