2.1 Fish shape and swimming motion

In order to capture the main parameters of BCF swimming while keeping the problem complexity manageable, we model the main body of the fish and its caudal fin but not the other fins or details on the body such as scales, finlets and other protrusions. We represent a swimming fish by a neutrally buoyant undulating body of length $L=1$ , as illustrated in figure 1. For the 2-D simulations, a NACA0012 shape is chosen at rest, whereas a danio-shaped body, shown in figure 2, is used for the 3-D simulations. The body propels itself at average speed $U_{s}$ in a fluid of kinematic viscosity $\unicode[STIX]{x1D708}$ and density $\unicode[STIX]{x1D70C}$ by oscillating its mid-line in the transverse direction $y$ . The leading edge of the body at rest is located at $x=0$ and its trailing edge at $x=1$ . In 2-D simulations, we will refer to the body as a ‘fish’ rather than a flexible foil to avoid confusion with the caudal fin, and with non-self-propelled flexible foils used as propulsors.

Figure 1. Schematic showing the fish model parameters. An elongated body of length $L$ undulates in a flow of speed $U_{s}$ with a wave travelling backward at speed $f\unicode[STIX]{x1D706}$ and amplitude $a$ at the trailing edge.

Figure 2. Three-dimensional fish geometry based on a giant danio.

We employ travelling wave kinematics that resemble those observed in fish according to either carangiform or anguilliform swimming, and include recoil. The lateral displacement, $h$ , of a point located at $x$ along the foil is given at time $t$ by the sum of body deformation $h_{0}$ , recoil $B$ and steering $y_{1}$ :

where $A(x)$ , with $A(1)=1$ , is the envelope of the prescribed backward travelling wave of wavelength $\unicode[STIX]{x1D706}$ and frequency $f$ ,

is the recoil term due to the hydrodynamic forces on the fish, and

can be used for steering (see appendix B) by adding camber to the fish, while $\unicode[STIX]{x1D6FE}$ and $\unicode[STIX]{x1D6FD}$ ensure that linear and angular momentum are conserved through the deformation. The recoil term represents periodic solid body transversal displacement and rotation. The amplitude of the translation is $a_{r}$ , while the amplitude of the rotation is $\arctan (b_{r})$ . The steering term $y_{1}$ is necessary to ensure stability but, in the steady regime, $y_{1}\ll a_{0}$ (typically $C<0.025$ ).

The parameter $a_{0}$ determines the amplitude of the deformation $h_{0}$ at the trailing edge. It is adjusted through a feedback control loop to ensure that the average net drag on the foil is 0, as described in appendix B. $h_{0}(x,t)$ can be used without the recoil and steering terms in order to fully prescribe the kinematics of the swimmer, in which case $h(x,t)=h_{0}(x,t)$ . For a freely moving body with prescribed deformation, the recoil is computed from the hydrodynamic forces on the body. In the latter case, the envelope of the actual displacement is given by $g(x)$ , with peak to peak amplitude at the trailing edge given by $a=2g(1)$ .

The prescribed kinematics of a carangiform swimmer, based on the experimental observation of steadily swimming saithe (Videler & Hess Reference Videler and Hess1984; Videler Reference Videler1993), is often modelled without recoil ( $B=0$ ) as:

This motion is for example used in Dong & Lu (Reference Dong and Lu2007), Borazjani & Sotiropoulos (Reference Borazjani and Sotiropoulos2008) and, in the rest of the paper, will be referred to as the carangiform gait. Assuming $B=0$ , experimental observations of American eels (Tytell & Lauder Reference Tytell and Lauder2004) provide that anguilliform motion can be represented by:

Figure 3(a) shows the prescribed envelope $A(x)$ for the carangiform (respectively anguilliform) swimmer defined in (2.4) (respectively (2.5)). Figure 3(b) illustrates the resulting mid-line displacement in the presence of the recoil term.

Figure 3. Carangiform and anguilliform motion for $f=1.8$ and $a_{0}=0.1$ at Reynolds number $\mathit{Re}=5000$ with recoil. (a) Prescribed amplitude envelopes, (b) mid-line displacement.

2.2 Kinematic parameters

The goal of this paper is to identify kinematic parameters that minimize the self-propelled swimming power $P_{in}$ for a given speed (Reynolds number) and body shape. In order to quantify the fitness of each motion, the quasi-propulsive efficiency $\unicode[STIX]{x1D702}_{QP}$ is used, which compares $P_{in}$ to the useful power, i.e. the resistance $R$ of the rigid–straight towed body at the same speed $U_{s}$ times the speed: $\unicode[STIX]{x1D702}_{QP}=RU_{s}/\overline{P_{in}}$ . Indeed, as discussed in Maertens, Triantafyllou & Yue (Reference Maertens, Triantafyllou and Yue2015), the Froude propulsive efficiency is not appropriate here as it is zero for a self-propelled body. Therefore, $R$ is calculated by towing a rigid body in the still flow with the prescribed velocity $U_{s}$ .

For rigid flapping foils, the parameters that characterize the motion and its performance have been extensively studied (Anderson et al. Reference Anderson, Streitlien, Barrett and Triantafyllou1998; Read, Hover & Triantafyllou Reference Read, Hover and Triantafyllou2003). The principal kinematic parameters are the Strouhal number and the maximum nominal angle of attack, and, to a lesser degree, the heave amplitude to chord ratio, and the phase angle between heave and pitch; all, typically, measured at 25 % of the chord. The Strouhal number is a wake parameter, since it characterizes the dynamics of the (unstable) wake (Triantafyllou et al. Reference Triantafyllou, Triantafyllou and Gopalkrishnan1991; Triantafyllou, Triantafyllou & Grosenbaugh Reference Triantafyllou, Triantafyllou and Grosenbaugh1993); hence the width of the wake must, in principle, be used as the characteristic length. However, the width of the wake is unavailable beforehand, so this characteristic length is approximated typically by the peak to peak motion of the trailing edge. Hence, for an undulating flexible foil, we define the Strouhal number, heave amplitude, pitch angle and nominal angle of attack at the trailing edge. These parameters and others used throughout this paper are summarized in table 1. While the motion cannot be characterized by these parameters alone, they play an important role in determining the swimming efficiency.

Changing the amplitude of motion and Strouhal number can be achieved through parameters like $a_{0}$ and $f$ (though, for a given motion and average velocity, there is a unique amplitude that ensures a steady velocity) but, in general, the pitch amplitude $\unicode[STIX]{x1D703}_{max}$ and maximum angle of attack $\unicode[STIX]{x1D6FC}_{max}$ cannot be directly controlled. Therefore, when optimizing the swimming gait, it is important to choose a parametrization that allows to adjust the pitch and angle of attack amplitudes independently of the heave amplitude and Strouhal number. This is best done by changing parameters that control the derivative of the prescribed envelope $A(x)$ at the trailing edge.

In this study, the wavelength of the travelling wave is fixed, equal to the body length, and the average swimming speed is adjusted to ensure a Reynolds number $\mathit{Re}=5000$ . Fish body shape is also fixed, while the linear and angular recoil terms are computed by integration of the hydrodynamic forces. The lateral flexing motion $h_{0}(x,t)$ (i.e. the lateral motion after the linear and angular recoil are subtracted) is characterized by four parameters: the frequency, amplitude and two parameters controlling the shape of $A(x)$ , the envelope of the unsteady bending motion. With prescribed frequency, the amplitude is automatically adjusted to satisfy self-propulsion at $\mathit{Re}=5000$ , while the two remaining parameters are varied systematically in order to identify the values that minimize the swimming power (equivalently, maximize the quasi-propulsive efficiency $\unicode[STIX]{x1D702}_{QP}$ ). In order to explore a wide range of motions, two different parametrizations are used for $A(x)$ : a quadratic envelope and a Gaussian envelope (see details in § 3.1).

2.3 Governing equations and numerical implementation

The motion of a self-propelled swimming body is determined by the coupled fluid–body dynamics. The physical parameters are non-dimensionalized by the fish body length $L$ , its intended average cruising speed $U_{s}$ and the density of water  $\unicode[STIX]{x1D70C}$ .

In order to solve the coupled fluid/body problem described above, we adapted the second-order boundary data immersion method (BDIM) presented in Maertens & Weymouth (Reference Maertens and Weymouth2015). The validation of the numerical method used in the present paper for simulating self-propelled undulating bodies is presented in appendix A. For the two-dimensional simulations, constant velocity $\boldsymbol{u}=\boldsymbol{U}_{s}$ is used on the inlet ( $x=-6$ ), periodic boundary conditions on the upper and lower boundaries ( $y=\pm 2.4$ ) and a zero gradient exit condition with global flux correction ( $x=7$ ). As shown in appendix A, the domain is large enough, such that the boundaries do not impact the results of the simulation. The Cartesian grid is uniform near the fish and uses a 2 % geometric expansion ratio for the spacing in the far field, as illustrated in figure 4. As discussed in appendix A, the grid spacing near the fish is $\text{d}x=\text{d}y=1/160$ for the optimization procedures, which allows reasonable accuracy and fast iterations. The optimized kinematics are then simulated using a finer grid with $\text{d}x=\text{d}y=1/320$ in order to guarantee an error of less than 10 %. Similarly, the fine grid is used for all the two-fish simulations.

The three-dimensional simulations are run on a $6\times 3\times 3$ domain with constant velocity $\boldsymbol{u}=\boldsymbol{U}_{s}$ on the inlet, a zero gradient exit condition with global flux correction, and periodic boundary conditions along $y$ and $z$ boundaries. The Cartesian grid is uniform near the fish with grid size $\text{d}x=\text{d}y=\text{d}z=1/100$ and uses a 4 % geometric expansion ratio for the spacing in the far field.

Figure 4. Flow configuration for the undulating NACA0012 simulations. The vorticity field for the carangiform motion with $f=1.8$ and zero mean drag is shown as an example.

The fluid and body equations are integrated over the fluid and body domains, respectively $\unicode[STIX]{x1D6FA}_{f}$ and $\unicode[STIX]{x1D6FA}_{b}$ , with a kernel of radius $\unicode[STIX]{x1D716}=2\,\text{d}x$ . The BDIM equations for the smoothed velocity field $\boldsymbol{u}_{\unicode[STIX]{x1D716}}$ are valid over the complete domain $\unicode[STIX]{x1D6FA}=\unicode[STIX]{x1D6FA}_{f}\cup \unicode[STIX]{x1D6FA}_{b}$ and enforce the no-slip boundary condition at the interface. These equations, integrated from time $t$ to time $t_{+}=t+\unicode[STIX]{x0394}t$ , are:

where $\boldsymbol{v}$ is the velocity field associated with the closest body, $\hat{n}$ the unitary normal to the closest fluid/solid boundary (pointing toward the fluid) and $d$ the signed distance to the closest boundary ( $d>0$ within the fluid, $d<0$ inside a body). $\unicode[STIX]{x1D707}_{0}^{\unicode[STIX]{x1D716}}$ and $\unicode[STIX]{x1D707}_{1}^{\unicode[STIX]{x1D716}}$ are respectively the zeroth and first central moments of the smooth delta kernel (see Maertens & Weymouth (Reference Maertens and Weymouth2015) for more details). $\boldsymbol{\unicode[STIX]{x2202}}\boldsymbol{P}_{\unicode[STIX]{x0394}t}$ , the pressure impulse, and $\boldsymbol{R}_{\unicode[STIX]{x0394}t}$ , accounting for all the non-pressure terms, are defined as:

In order to simplify the equations of motion, we consider motion within the $(x,y)$ plane, such that the translational velocity of the body centre of mass (COM), $\boldsymbol{v}_{c}$ , is a two-dimensional vector $(v_{c}^{x},v_{c}^{y})$ , and its rotation velocity is $\unicode[STIX]{x1D714}_{b}=\unicode[STIX]{x1D714}_{b}^{z}$ . We then define the generalized velocity $\boldsymbol{V}$ , location $\boldsymbol{X}$ , and force $\boldsymbol{F}$ vectors, as well as the generalized mass matrix $\unicode[STIX]{x1D648}$ :

where $\boldsymbol{F}_{h}$ is the hydrodynamic force on the body, $m$ is the mass of the body which has density $\unicode[STIX]{x1D70C}_{b}=\unicode[STIX]{x1D70C}$ and $I_{c}$ its moment of inertia with respect to the COM, computed at each time step. The motion of the body is governed by:

The coupled dynamic equations are discretized using a sequentially staggered Euler explicit integration scheme with Heun’s corrector. Sequentially staggered schemes are computationally efficient but, for large added mass, they become unconditionally unstable (Förster, Wall & Ramm Reference Förster, Wall and Ramm2007), regardless of the particular scheme used. In order to stabilize the numerical scheme, we introduce the virtual added mass matrix  $\unicode[STIX]{x1D648}_{\boldsymbol{a}}$ .

The virtual added mass, which is used in an implicit added mass scheme (Connell & Yue Reference Connell and Yue2007; Zhu & Shoele Reference Zhu and Shoele2008; Peng & Zhu Reference Peng and Zhu2009), can eliminate the instability due to large added mass, but its exact value will not affect the results. In the case of an undulating fish, the coefficients of the matrix can be estimated from the added mass of the straight fish body at zero angle of attack, which has been done here. In the present simulations, only the terms necessary to the stability have been used and it has been assessed that the exact value chosen for the added mass does not impact the results, as long as it is large enough to ensure stability but not so large as to damp the system. The virtual added mass used is a diagonal matrix with value $[0~11m~13m]$ .

We also define the total mass as:

With these new definitions, we integrate (2.9) over a time step $\unicode[STIX]{x0394}t$ in the form:

At each time step $t_{n}$ , the fluid and body velocities, $\boldsymbol{u}_{n}=\boldsymbol{u}_{\unicode[STIX]{x1D716}}(t_{n})$ and $\boldsymbol{v}_{n}=\boldsymbol{v}(t_{n})$ respectively, are calculated from the velocities and forces at the previous time steps according to (2.6) and (2.11) (for additional details, see Maertens Reference Maertens2015).

2.4 The importance of recoil

Equation (2.11) determines the recoil $B(x,t)$ resulting from the prescribed flexing $h_{0}(x,t)$ and the related hydrodynamic forces $\boldsymbol{F}$ . Due to the significant added complexity incurred by the recoil term, most of the earlier simulation studies neglected it (Dong & Lu Reference Dong and Lu2007; Borazjani & Sotiropoulos Reference Borazjani and Sotiropoulos2008). However, the amplitude of this term, and its impact on the estimated swimming power, are substantial (Reid et al. Reference Reid, Hildenbrandt, Padding and Hemelrijk2012), as illustrated below.

We consider first the carangiform motion of (2.4) with frequency $f=2.1$ . Figure 5(a) shows the dimensionless linear and angular momentum for the self-propelled fish, including the recoil, as determined by the hydrodynamic forces and adaptive amplitude $a_{0}$ . The angular and transverse momentum are larger than the longitudinal momentum, but the three amplitudes are comparable. However, the non-dimensional moment of inertia of the fish is much smaller than its mass:

where the mass and moment of inertia are non-dimensionalized by the length $L$ and density $\unicode[STIX]{x1D70C}$ . Therefore, whereas the linear momentum results in velocities smaller than 3 % of the free-stream $U_{s}$ , the rotation of the fish generates velocities at the trailing edge up to 40 % of the free stream, as shown in figure 5(b). This observation suggests that, whereas the longitudinal motion of the fish might be negligible, the transverse motion and specifically the angular recoil, are important.

Figure 5. (a) Linear and angular momentum and (b) corresponding velocities for a neutrally buoyant self-propelled NACA0012 with carangiform motion at frequency $f=1/T=2.1$ .

In order to further illustrate this result, figure 6 shows the quasi-propulsive efficiency as a function of frequency for the carangiform and anguilliform motions with and without recoil. The figure shows that, at all frequencies, the undulation with recoil requires up to 40 % more power than the undulation without recoil. Therefore, simulations that do not allow for recoil are likely to underestimate the swimming power, as discussed in Reid et al. (Reference Reid, Hildenbrandt, Padding and Hemelrijk2009). The figure also shows that the optimal frequency without recoil might differ from the optimal frequency with recoil. In the cases studied here, the optimal frequency for the carangiform undulation without recoil is around $f=1.6$ , while with recoil it is around $f=2.1$ .

Figure 6. Quasi-propulsive efficiency as a function of frequency for the carangiform and anguilliform motions with and without recoil.

In summary, we have shown here that the impact of the angular recoil of the fish on swimming performance is significant. In order to estimate meaningful values of fish swimming efficiency, it is critical to allow for recoil.

2.5 Gait optimization procedure

Toki & Yue (Reference Toki and Yue2012) and Eloy (Reference Eloy2013) combined an evolutionary algorithm with Lighthill’s potential flow slender-body model to simultaneously optimize the shape and kinematics, using, respectively, 22 and 9 parameters. In this paper, we employ viscous simulations that are far more demanding computationally, hence we use only two parameters, which allows us to find an optimum with a reduced number of evaluations, and also facilitates the visualization and interpretation of the results.

For a given kinematic parametrization and frequency, the two parameters controlling the envelope $A(x)$ are optimized using derivative-free optimization (Rios & Sahinidis Reference Rios and Sahinidis2013). We apply the BOBYQA algorithm that performs bound-constrained optimization using an iteratively constructed quadratic approximation for the objective function (Powell Reference Powell2009). For each set of parameters, the viscous simulation is run for 15 non-dimensional time units, and the average power coefficient $\overline{C_{P}}$ across the last 10 undulation periods is calculated. Note that, as can be seen on figure 36(b) of appendix A, the instantaneous power coefficient can be negative. Indeed, animals have the capability to store energy in their tendons and hence recover energy during the oscillatory cycle (Roberts & Azizi Reference Roberts and Azizi2011). Although this recovery may not be 100 %, under optimal conditions our estimate averaging the power coefficient regardless of the instantaneous sign provides a good estimate of the power consumed.

Based on the values of $\overline{C_{P}}$ , the implementation of BOBYQA provided by the NLopt free C library (Johnson Reference Johnson2013) interfaced with Matlab computes the next set of parameters. In order to avoid finding a local minimum due to numerical noise, after the algorithm has converged, it is run again using the previously found minimum as a starting point.

3.1 Gait optimization for a two-dimensional foil

For several values of undulation frequency, we optimize the deformation envelope $A(x)$ . $A(x)$ has been traditionally modelled by a quadratic function, of the form:

Unlike $c_{1}$ and $c_{2}$ , $A(0)$ and $A(1/2)$ can be restricted to a rectangle and are easy to interpret as the envelope amplitude at the leading edge and mid-chord respectively. Therefore, in the figures, each envelope is represented by $A(0)$ and $A(1/2)$ , the amplitude at the trailing edge being by definition $A(1)=1$ .

First, we fix the undulation frequency to $f=1.8$ and optimize the quadratic envelope $A(x)$ , restricting $A(0)$ to positive values. Figure 7(a) shows the efficiency as a function of $A(0)$ and $A(1/2)$ . The carangiform envelope used in previous sections is denoted by a black square, and the anguilliform gait through a diamond. It is clear that the envelopes above the dashed line, which are concave envelopes with a peak upstream from the trailing edge, have good efficiency. The efficiency decreases very quickly below the dashed line, as the envelope becomes convex with an increasing amplitude at the trailing edge. Therefore, the envelope traditionally used to model carangiform swimming is inefficient, whereas the anguilliform envelope, which is closer to a straight line, is much more efficient. Among the concave envelopes, $A(0)=0$ seems best, together with $1\leqslant A(1/2)\leqslant 1.7$ , where the efficiency reaches a value of 48 %. Since the optimal quadratic gait saturates the constraint $A(0)\geqslant 0$ , we then fix the leading edge amplitude to $A(0)=0$ and optimize the undulation frequency $f$ and the second envelope parameter $A(1/2)$ . Figure 7(b) shows the efficiency as a function of $f$ and $A(1/2)$ . Here again, around the optimal point, the efficiency is not very sensitive to the exact value of $f$ and $A(1/2)$ . The optimal quadratic envelope ( $A(0)=0$ , $A(1/2)=1$ , $A(1)=1$ ) has a maximum amplitude at $x=3/4$ and reaches an efficiency of $\unicode[STIX]{x1D702}_{QP}=49\,\%$ around $f=1.6$ .

Figure 7. $\unicode[STIX]{x1D702}_{QP}$ as a function of $A(0)$ or $f$ and $A(1/2)$ for quadratic envelopes. The black dots show the location of the points that have been used to build the thin-plate smoothing spline (tpaps function in Matlab with smoothing parameter $p=0.999$ ) represented in colour. (a) Fixed frequency $f=1.8$ . The carangiform and anguilliform motions are respectively denoted by a black square and a black diamond, and a dashed line shows the location of linear envelopes (points below this line correspond to convex envelopes, while above it the envelopes are concave). (b) Fixed leading edge value $A(0)=0$ . Note that the thin-plate smoothing is only accurate in regions with a sufficient density of points.

Traditional models for swimming motion ignore recoil and therefore identify body deformation with displacement. However, while a quadratic convex envelope can reasonably be used to describe the displacement envelope of undulating fish, the curvature envelope in saithe and mackerel has a distinctive peak around the peduncle section (Videler & Hess Reference Videler and Hess1984). The results from figure 7 also suggest that the efficiency is higher if the deformation is largest upstream of the trailing edge rather than at the trailing edge itself. Such envelopes can be better modelled by a Gaussian function of the form:

where $x_{1}$ parametrizes the location of the peak and $\unicode[STIX]{x1D6FF}$ its width, as shown in figure 8. With the Gaussian function, it is easy to change the pitch and angle of attack amplitudes at the tail by adjusting the location and width of the peak. Since the Gaussian envelope is always positive, the entire $(x_{1},\,\unicode[STIX]{x1D6FF})$ space can be used to search for an optimal gait without running into degenerate gaits.

Figure 8. Definition of the parameters for a Gaussian envelope.

Figure 9 shows the curvature amplitude for various envelopes. The curvature amplitude $C(x)$ is given by:

It can in particular be seen that while the traditional carangiform envelope leads to a curvature that is minimum at 0.25 from the leading edge and maximum at the trailing edge. This is very different from the curvature observed in saithe and mackerel (Videler & Hess Reference Videler and Hess1984). A Gaussian envelope, on the other hand, can lead to a very small curvature for $x<0.4$ with a peak around $x=0.8{-}0.9$ . This curvature profile is in much better agreement with the curvature reported in Videler & Hess (Reference Videler and Hess1984). Moreover, the curvature peak corresponds to the peak of the Gaussian, which will later be referred to as the area where the body deformation is largest.

Figure 9. (a) Deformation and (b) curvature for various quadratic and Gaussian envelopes. The carangiform envelope is the one defined in (2.4): $A(x)=1-0.825(x-1)+1.625(x^{2}-1)$ . The optim. quadratic envelop is the most efficient quadratic envelop identified in figure 7: $A(x)=1+3(x-1)-2(x^{2}-1)$ . The optim. Gauss envelopes are the most efficient Gaussian envelopes identified for $f=1.5$ and $f=2.4$ respectively, with parameters summarized in table 2.

Figure 10 shows the efficiency as a function of $x_{1}$ and $\unicode[STIX]{x1D6FF}$ in the neighbourhood of the optimal envelope (in the sense of requiring minimum input power) for five frequencies ranging from $f=1.5$ to $f=2.7$ . For all frequencies, the efficiency decreases very rapidly as $\unicode[STIX]{x1D6FF}$ is decreased below its optimal value, while the efficiency is much less sensitive to increases above this optimal value. Moreover, while for all frequencies it is possible to find a region in the $(x_{1},\,\unicode[STIX]{x1D6FF})$ space that reaches an efficiency of 50 % (see table 2 for details), the optimal envelope clearly depends on the frequency. Note that, for a given motion, the angle of attack also changes with frequency. Since we hope to find an efficient angle of attack by optimizing the envelope, it is not surprising that the optimum is frequency dependent.

Figure 10. $\unicode[STIX]{x1D702}_{QP}$ as a function of $x_{1}$ and $\unicode[STIX]{x1D6FF}$ near the optimum for Gaussian envelopes. (a) $f=1.5$ , (b) $f=1.8$ , (c) $f=2.1$ , (d) $f=2.4$ , (e) $f=2.7$ , (f) colour bar. The black dots show the location of the points that have been used to build the thin-plate smoothing spline (tpaps function in Matlab with smoothing parameter $p=0.999$ ) represented in colour. Note that the thin-plate smoothing is only accurate in regions with a sufficient density of points.

At low frequency, gaits with undulations of the entire body ( $x_{1}=0.73$ and $\unicode[STIX]{x1D6FF}=0.52$ at $f=1.5$ ) are most efficient, while at high frequency, the undulations should be restricted to a narrow region ( $\unicode[STIX]{x1D6FF}=0.21$ at $f=2.7$ ) located around 25 % of the trailing edge ( $x_{1}=0.88$ at $f=2.7$ ). However, for all frequencies, the optimized deformation envelope $A(x)$ , shown in figure 11(a), is qualitatively similar to the curvature envelope from Videler & Hess (Reference Videler and Hess1984), with a small amplitude at the leading edge, a peak of amplitude 0.1 located 10 %–30 % from the trailing edge, and a sharp decrease in amplitude at the trailing edge. The corresponding displacement envelopes $g(x)$ are shown in figure 11(b). The displacement envelopes are qualitatively similar to the carangiform displacement envelope from Videler & Hess (Reference Videler and Hess1984), with an amplitude close to $g(0)=0.02$ at the leading edge (for from 1.8 to 2.7), a minimum amplitude around $x=0.25$ and a maximum amplitude at the trailing edge.

Figure 11. Optimized (a) prescribed deformation envelopes and (b) displacement envelopes for the Gaussian parametrization. $\unicode[STIX]{x1D706}=1$ and $f=[1.5,1.8,2.1,2.4,2.7]$ .

The similarity between optimized envelopes (both deformation and displacement) and envelopes observed in fish suggest, while the motion traditionally used to model carangiform swimming is inefficient, the actual motion of carangiform swimmers could be close to optimal. It is also interesting to note that, since quadratic envelopes can only result in functions with a wide peak, they can reach the same efficiency as the Gaussian envelopes at low frequency ( $f=1.5$ ), but not at high frequency ( $f>2$ ) where a sharp peak is needed to mitigate the large pitch and angle of attack caused by fast transverse motion of the tail.

3.2 Characterization of efficient swimming gaits for a two-dimensional fish

We showed in the previous section that, by changing the location and width of the peak in a Gaussian deformation envelope, a very efficient gait can be designed for a large range of undulation frequencies.

Figure 12 shows the deformed fish for three optimized gaits. As expected from figure 11, at $f=1.5$ the entire length of the fish undergoes noticeable deformation and displacement, resulting in a swimming motion that is similar to that of an anguilliform swimmer, with a moderate curvature along the entire body: the deformation of the fish matches the large wavelength trajectory of the trailing edge and thus avoids the efficiency loss associated with a large angle of attack. At higher frequency, the front half of the fish undergoes virtually no deformation, resulting in a swimming motion very similar to a carangiform ( $f=2.1$ ) or even a thunniform ( $f=2.7$ ) swimmer. At this frequency, the region that would correspond to the peduncle deforms with a very large curvature caused by the sharp peak in the Gaussian envelope. This allows the deformation of the fish to match the small wavelength trajectory of the trailing edge and thus avoid the efficiency loss associated with a large angle of attack.

Figure 12. Superimposed body outlines over one undulation period for three frequencies.

Figure 13 shows the deformed fish and vorticity snapshots for the five optimized gaits at $t/T=0~(\text{mod }1)$ , where $T=1/f$ is the undulation period. With a Gaussian deformation envelope, a peak width specifically tailored to the undulation frequency allows for reduced angle of attack at all frequencies. Specifically, as the undulation frequency increases, the angle of attack at the trailing edge, induced by the tail displacement increases. In order to recover a small and efficient angle of attack, the (already negative) slope of the deformation envelope must be decreased, which can only be done by a sharper peak. As for thrust-producing flapping foils, a reverse Kármán vortex street forms in the wake. The width and wavelength of the reverse Kármán vortex street decreases with increasing undulation frequency, and secondary small vortices develop at low frequency.

Figure 13. Snapshots of vorticity for optimized gaits at $t/T=0~(\text{mod }1)$ . (a) $f=1.5$ , (b) $f=1.8$ , (c) $f=2.1$ , (d) $f=2.4$ , (e) $f=2.7$ , (f) colour bar.

Table 2 summarizes the parameters and properties of the five optimized gaits. The quasi-propulsive efficiency $\unicode[STIX]{x1D702}_{QP}$ of these undulatory gaits is of prime interest. The efficiency reaches 57 % for $f=2.7$ , whereas the least efficient frequency, $f=1.5$ , reaches $\unicode[STIX]{x1D702}_{QP}=49\,\%$ . Another important parameter is the Strouhal number, which is close to $St=0.35$ . The consistency of the Strouhal number for the optimized envelopes across frequencies suggests that, for a given Reynolds number, there exists an optimal Strouhal number that can be reached with a large range of frequencies. Like the Strouhal number, the maximum pitch angle $\unicode[STIX]{x1D703}_{max}$ and maximum angle of attack $\unicode[STIX]{x1D6FC}_{max}$ are almost constant across the five optimized gaits, with a value close to $\unicode[STIX]{x1D703}_{max}=31^{\circ }$ and $\unicode[STIX]{x1D6FC}_{max}=17^{\circ }$ . The corresponding phase angle between the heave and pitch of the trailing edge is $\unicode[STIX]{x1D713}=82^{\circ }$ . The results from this optimization show that, as in rigid flapping foils, the efficiency of undulating fish is primarily driven by the Strouhal number, angle of attack, heave motion (or pitch motion) and heave–pitch phase angle, all at the trailing edge. The optimal Strouhal number, pitch angle and angle of attack can be attained by tuning the envelope peak for each frequency.

Figure 14 shows the pressure field and body velocity for the optimized envelopes with frequency $f=[1.5,2.1,2.7]$ at their respective time of minimum and maximum power. For $f=1.5$ (figure 14 a,b) and $f=2.7$ (figure 14 e,f), there are three distinct sections along the upper side of the fish: high pressure near the leading edge, low pressure in the middle and high pressure near the trailing edge (and the opposite on the other side). In figures 14(b), 14(f), these sections almost exactly match transverse velocity of respectively positive, negative and positive sign, resulting in very large instantaneous swimming power. Conversely, in figure 14(a,e), the sign of the transverse velocity is reversed, resulting in a significant negative swimming power. For $f=2.1$ , the pressure changes along the fish are smaller, and the pressure is close to zero along a large portion of the fish. Moreover, unlike for $f=2.7$ , the sign changes in pressure do not match the sign changes in transverse velocity. For instance, at $t/T=0$ , the pressure along the bottom side of the fish near the trailing edge is positive (not shown here), which would result in a positive swimming power. Therefore, the minimum power is reached at a later time $t/T=0.12$ , at which point the amplitude is largest in areas where the pressure is close to zero, resulting in a very small power. Similarly, the maximum power reached at $t/T=0.34$ is not as large as for $f=2.7$ because the sections of high pressure do not exactly match the sections of large transverse velocity.

Figure 14. Snapshots of pressure field with arrows showing the body velocity. (a,b) Optimized Gaussian envelope at $f=1.5$ ; (c,d) optimized Gaussian envelope at $f=2.1$ ; (e,f) optimized Gaussian envelope at $f=2.7$ . (a,c) $t/T=0.12~(\text{mod }1)$ (minimum power for $f=1.5$ and $f=2.1$ ); (e) $t/T=0~(\text{mod }1)$ (minimum power for $f=2.7$ ); (b,d) $t/T=0.34~(\text{mod }1)$ (maximum power for $f=1.5$ and $f=2.1$ ; (f) $t/T=0.29~(\text{mod }1)$ (maximum power for $f=2.7$ ).

It must be pointed out that there are additional parameters affecting the efficiency of an undulating fish, since the efficiency ranges from $\unicode[STIX]{x1D702}_{QP}=0.49$ at $f=1.5$ to $\unicode[STIX]{x1D702}_{QP}=0.57$ at $f=2.7$ . This trend runs contrary to results from Shen et al. (Reference Shen, Zhang, Yue and Triantafyllou2003) who found that a slip ratio around $s_{r}=0.8$ ( $f=1.2$ ) is optimal for a body undergoing travelling wave motion of constant amplitude, in order to reduce separation effects and turbulence intensity. In our case, however, these results do not strictly apply because the undulations have a non-constant envelope, and especially for higher frequency are confined to a small section of the fish.

3.3 Application to a three-dimensional danio-shaped swimmer

We have so far modelled a fish by a two-dimensional fish-like flexible foil. However, fish have a highly three-dimensional geometry. In particular, most carangiform and thunniform swimmers are characterized by a region of reduced depth, around 20 % from the trailing edge, the peduncle. In order to model this region of reduced added mass with a two-dimensional geometry, it might be more appropriate to model a fish with a separate foil for the tail, as illustrated in figure 15. The fish model shown in this figure undulates with the optimized gait at frequency $f=2.4$ identified earlier, and the performance ( $\unicode[STIX]{x1D702}_{QP}=0.54$ ) is very close to that obtained with a single foil, indicating that the results are robust to changes in the geometry.

Figure 15. Snapshot of the vorticity field around a two-dimensional fish with a separate tail fin.

In the rest of this section, we consider a simplified three-dimensional fish shape, shown in figure 2, which is based on a giant danio (Devario aequipinnatu). For this geometry, we fix the undulation frequency to $f=2.4$ and optimize a Gaussian envelope for quasi-propulsive efficiency (for a fixed swimming speed $U_{s}$ , we minimize the expended power $\overline{P_{in}}$ ). In figure 16 we compare how $\unicode[STIX]{x1D702}_{QP}$ changes with the envelope parameters $x_{1}$ and $\unicode[STIX]{x1D6FF}$ for a two-dimensional fish and for the three-dimensional shape. The efficiency is generally lower with the three-dimensional shape, but the dependence on $x_{1}$ and $\unicode[STIX]{x1D6FF}$ is very similar for both geometries: the most efficient gaits are for $0.8<x_{1}<0.9$ and $0.2<\unicode[STIX]{x1D6FF}<0.3$ with a sharp decrease in efficiency for $\unicode[STIX]{x1D6FF}<0.2$ . This shows that, even though three-dimensional effects reduce the swimming efficiency, most of the conclusions about BCF swimming drawn from the two-dimensional study extend to three-dimensional shapes.

Figure 16. $\unicode[STIX]{x1D702}_{QP}$ as a function of $x_{1}$ and $\unicode[STIX]{x1D6FF}$ near the optimum for (a) 2-D and (b) 3-D geometries with $f=2.4$ . The black dots show the location of the points that have been used to build the thin-plate smoothing spline (tpaps function in Matlab with smoothing parameter $p=0.999$ ). Note that the thin-plate smoothing is only accurate in regions with a sufficient density of points.

The parameters and properties of the optimized gait for $f=2.4$ are compared to those of the carangiform gait $A_{c}(x)$ in table 3. Like in the 2-D case, the optimization decreases the power consumption by 50 % compared with the carangiform gait. As in two dimensions, the optimized gait manages to bring the phase angle between the heave and pitch motion of the trailing edge close to $90^{\circ }$ , which significantly reduces the angle of attack. As a result, the optimized gait for the three-dimensional fish shape have a pitch angle, phase angle and angle of attack very close to the optimized gait for the two-dimensional fish. However, since the 3-D effects reduce the thrust produced by the undulating motion, the Strouhal number is higher than in two dimensions, especially for the carangiform gait.

Figure 17. Prescribed deformation envelope $a_{0}A(x)$ and displacement envelope $g(x)$ for (a) carangiform gait with $f=3$ and (b) optimized gait with $f=2.4$ . Frequencies have been selected such that the displacement amplitude at the trailing edge is the same.

Figure 18. Superimposed body outlines over one undulation period for (a) the carangiform motion and (b) the optimized gait.

Figure 17 shows the deformation envelope $A(x)$ and the displacement envelope $g(x)$ for the carangiform gait at $f=3$ and for the optimized gait. Similarly to the observations made for the two-dimensional body shape, the optimal gait has maximum displacement at the trailing edge, but maximum deformation (curvature) around $x=0.85$ , which corresponds to the peduncle.

The superimposed body outlines for the optimized gait shown in figure 18(b) also look very similar to the body outlines of the optimized motions in two dimensions: the deformation of the tail follows the trajectory of the trailing edge, resulting in an efficient low angle of attack. The body outlines for the carangiform motion, on the other hand, show that the pitch of the tail is out of phase with its velocity (phase angle far from $90^{\circ }$ ), which results in a very inefficient gait, with a large angle of attack.

Finally, we show the flow structure around the 3-D fish model for both gaits in figure 19. The performance difference between the two gaits is accompanied by noticeable differences in the wake structure of the two swimmers. For both gaits, figures 19(a) and 19(b) show wakes comprised of two interconnected vortex loops per cycle, together with other smaller structures. In particular, the structure in the wake of the optimized motion is complex, with many vortex tubes interlaced with each other. Indeed, as can also be seen in the vorticity field at $z=0$ (figure 19 d), the deformation at the peduncle is quite large for the optimized gait, resulting in vortex tubes separating from the main body and then interacting with the structures shed from the tail.

Figure 19. Snapshots of the flow around a three-dimensional fish with (a,c,e,g) a carangiform and (b,d,f,h) an optimized gait. (a,b) Three-dimensional vortical structures visualized using the $\unicode[STIX]{x1D706}_{2}$ -criterion; (c,d) $z$ component of the vorticity in the $z=0$ plane; (e,f) pressure in the $z=0$ plane; (g,h) pressure in the $z=0.06$ plane.

Borazjani & Sotiropoulos (Reference Borazjani and Sotiropoulos2010) also observed in their 3-D simulations that the 2-D wake structure, dominated by a single vortex pair (or ring in three dimensions), transitions to vortex loops wrapping around each other for Strouhal number greater than $St=0.3$ . Dong, Mittal & Najjar (Reference Dong, Mittal and Najjar2006) showed that the same phenomenon happens to elliptical flapping foils of finite aspect ratio: at low aspect ratio/large Strouhal number, two vortex rings are shed each cycle. As the aspect ratio increases or the Strouhal number decreases, the tip vortices do not merge together any more and the wake consists of interconnected loops. As the Strouhal number further decreases or the aspect ratio increases, the three-dimensional effects become even weaker and the linkage between tip vortices disappears. At this point, the 3-D wake looks similar to the (reverse) Kármán vortex street observed in two dimensions.

In the carangiform example shown here, the tip vortices merge, while with the optimized gait, which has a lower Strouhal number and angle of attack, they do not. At higher Reynolds number, the Strouhal number would be smaller and a wake similar to that observed in two dimensions would probably emerge. Figure 19(c,d) shows that, near the tail, the vorticity in the $z=0$ plane looks very similar to that behind a 2-D foil. However, under the influence of the tip vortices, the vortex sheets shed by the tail do not evolve into two strong vortices as in two dimensions. As a result, whereas the pressure field around the undulating fish-like shape is very similar to the pressure around an undulating airfoil, the pressure signature in the wake shown at the plane $z=0$ is very weak (figure 19 e,f). However, the pressure signature in the plane $z=0.06$ , just above the peduncle, is much stronger (figure 19 g,h), and could potentially be used by a downstream fish to reduce its swimming energy.

Figure 20 shows a magnified view of the vortex structures generated by the carangiform motion. A red line shows the formation of a clear vortex ring at the trailing edge of the tail between figures 20(a) and 20(c). In figure 20(e), the vortex ring is fully formed and detached from the tail. Since the vortex rings are oblique, they produce a large transverse velocity, which is inefficient and results in waste of energy. We also see a spanwise narrowing of the vortex rings as they convect downstream, as also observed in the simulations of Blondeaux et al. (Reference Blondeaux, Fornarelli, Guglielmini, Triantafyllou and Verzicco2005) and Dong et al. (Reference Dong, Mittal and Najjar2006) for a respectively rectangular and elliptical pitching and heaving foil.

Figure 20. (a,c,e) Side view and (b,d,f) top view of the vortex structures at several time steps for the carangiform gait. (a,b) $t/T=0.1~(\text{mod }1)$ ; (c,d) $t/T=0.4~(\text{mod }1)$ ; (e,f) $t/T=0.7~(\text{mod }1)$ . A red line shows the formation of a vortex ring.

Figure 21 shows a magnified view of the vortex structures generated by the optimized gait. The structure of the wake is more intricate than for the carangiform motion. In particular, instead of one set of interconnected vortex tubes, there are two sets of tubes, marked in red and green in the figure. The loop marked in red is the same as observed for the carangiform gait, but at this lower Strouhal number, it never fully closes into a clearly defined vortex ring. The tubes marked in green are formed upstream of the tail and are shed from the body as a result of the large curvature at the peduncle. The resulting vortex tubes are interlaced with the vortex loops from the tail with which they have a phase difference of close to $180^{\circ }$ .

Figure 21. (a,c,e) Side view and (b,d,f) top view of the vortex structures at several time steps for the optimized gait. (a,b) $t/T=0.1~(\text{mod }1)$ ; (c,d) $t/T=0.4~(\text{mod }1)$ ; (e,f) $t/T=0.7~(\text{mod }1)$ . A red line shows a vortex shed from the tail that never fully develops into a ring, while green lines show the vortices shed from the body.

For a three-dimensional fish shape with two-dimensional undulation, as for a two-dimensional foil, the Strouhal number, pitch angle (or heave motion), nominal angle of attack and phase angle at the trailing edge are the primary parameters for efficient swimming. The key is to decrease the (already negative) slope of the deformation envelope in order to reduce the angle of attack to efficient values. These two-dimensional considerations only depend on fish geometry through its impact on recoil. Moreover, the optimization leads to a lower Strouhal number and angle of attack, which reduces the three-dimensional effects observed behind the non-optimized gait, such as formation of inefficient oblique vortex ring chains. With the optimized gait, the production of thrust is also distributed between the body and the tail, which sheds vortex structures with opposite phase. It has been shown with turbines, for instance, that distributing the thrust production (or energy capture) could significantly improve the efficiency, and it is possible that fish use the same process. Finally, while we used a simplified fish geometry with a two-dimensional undulation, fish can also rely on three-dimensional motion of their dorsal and pectoral fins to fine tune their swimming motion and save energy (Drucker & Lauder Reference Drucker and Lauder2001; Lauder & Madden Reference Lauder and Madden2007).

4.1 Flow in the wake of a self-propelled undulating fish

The reverse Kármán vortex street behind a self-propelled undulating fish is characterized by its periodic structure, with vortices moving parallel to the $y=0$ axis in stable formation. The vorticity snapshot in figure 22(a) shows that the vortices decrease in strength under the effect of diffusion, but this is a slow process and the wake is primarily characterized by its periodicity. Figure 22(b) shows that the vortices generate patches of faster flow along the $y=0$ axis and slower flow away from the centreline. The time-averaged flow (figure 22 c,d) is characterized by four shear layers of alternating sign, resulting in a jet along the centreline with strips of slowed-down flow on either side; consistent with the momentumless wake, on average, of a self-propelled body (Triantafyllou et al. Reference Triantafyllou, Triantafyllou and Grosenbaugh1993).

Figure 22. Wake behind a self-propelled undulating fish for the optimized gait with Gaussian envelope and $\unicode[STIX]{x1D706}=1$ at frequency $f=1.5$ . (a) Instantaneous vorticity field; (b) instantaneous $x$ -velocity field; (c) time-averaged vorticity field; (d) time-averaged $x$ -velocity field.

The vorticity at longitudinal distance $d$ from the trailing edge in the wake of an undulating fish can be modelled as:

where the frequency $f$ is given by the undulation frequency and the wavelength $\unicode[STIX]{x1D706}_{w}$ and phase $\unicode[STIX]{x1D719}_{w}$ of the wake need to be determined. For the five optimized gaits with Gaussian envelope and $\unicode[STIX]{x1D706}=1$ presented in § 3.1, we estimated from the vorticity field the phase $\unicode[STIX]{x1D719}_{1}$ at several distances $d$ along the wake. In figure 23(a) we show the phase $\unicode[STIX]{x1D719}_{1}$ as a function of the distance $d$ , as well as the least squares linear fit for each swimming gait. The coefficients for the linear fit are summarized in table 4. For all the gaits, the phase is essentially proportional to the distance $d$ , with a coefficient of proportionality very close to the undulation frequency $f$ . Since $\unicode[STIX]{x1D706}_{w}=fc_{w}$ , where $c_{w}$ is the speed at which the vortices travel in the wake, this result shows that the vortices travel at the same speed as the free stream. Moreover, for the five gaits considered, the phase at $d=0$ is approximately 0.25, which means that the vortices are shed by the fish when the trailing edge has maximum transverse velocity. Finally, we confirm these observations by plotting $\unicode[STIX]{x1D719}_{1}$ as a function of $fd$ in figure 23(b). Assuming $c_{w}=1$ , the least squares estimate ( $\pm$ standard deviation) of the phase $\unicode[STIX]{x1D719}_{w}$ is:

From now on, $\unicode[STIX]{x1D706}_{w}=1/f$ and $\unicode[STIX]{x1D719}_{w}=0.25$ will be used to estimate the phase $\unicode[STIX]{x1D719}_{1}$ encountered by a downstream fish whose leading edge is located at distance $d$ from the upstream fish.

Figure 23. Vorticity phase in the wake of self-propelled undulating fish as a function of (a) distance, $d$ , and (b) frequency times distance, $fd$ . For each frequency, the optimized gait with Gaussian envelope is used.

4.2 Effect of phase and distance for two fish in line

Here we consider two fish-like foils following each other and undulating at frequency $f=1.5$ with the optimized gait for this frequency, as illustrated in figure 25. The amplitude of undulation, $a_{0}$ , is adjusted independently for each fish to ensure that both fish are in a stable position and produce zero net thrust on average. We vary the distance $d$ between the trailing edge of the upstream fish and the leading edge of the downstream fish, as well as $\unicode[STIX]{x1D719}$ , the phase of the downstream fish motion as defined in (2.1). An important parameter will be $\unicode[STIX]{x0394}\unicode[STIX]{x1D719}$ , the phase difference between the motion of the downstream fish leading edge and the vortices it encounters: $\unicode[STIX]{x0394}\unicode[STIX]{x1D719}=\unicode[STIX]{x1D719}-\unicode[STIX]{x1D719}_{1}(d)$ . In order to measure the impact of the pair configuration on each fish, we define $R(C_{P})$ (respectively $R(a_{0})$ ), the ratio of the power coefficient (respectively amplitude) in the pair configuration over the power coefficient (respectively amplitude) for the corresponding gait in open water.

Since the goal is to minimize expended power, $R(C_{P})$ is what should be minimized. However, it is interesting to determine whether a reduction in expended power results from a reduction in drag or from a pure reduction in swimming power. Since both fish are self-propelled, the drag is always zero, so the net drag on the fish cannot be used as an indicator of drag reduction. Instead we use the swimming amplitude as an indicator of drag: fish swim in order to overcome drag; if drag is reduced, then a smaller swimming motion will be sufficient to overcome the drag.

Both fish can benefit from swimming in pair, but the trends are very different. The swimming power and amplitude of the upstream fish is virtually independent of the phase of the downstream fish, as shown in figure 24(a). For large separations $d$ , the downstream fish does not impact the upstream fish whose efficiency is then almost the same as in open water. However, as the downstream fish gets close ( $d<0.5$ ), the high pressure around the leading edge of the downstream fish ‘pushes’ the upstream one, regardless of their phase difference. As a result, the upstream fish can reduce its swimming amplitude, expending less power than it would in open water. At $d=0.25$ , the undulation amplitude is reduced by 10 %, resulting in 28 % energy saving, corresponding to a quasi-propulsive efficiency (based on the towed drag in open water) of $\unicode[STIX]{x1D702}_{QP}=69\,\%$ .

Figure 24. Time-averaged power coefficient $\overline{C_{P}}$ and amplitude $a_{0}$ for (a) the upstream fish as a function of distance $d$ and (b) the downstream fish as a function of phase $\unicode[STIX]{x0394}\unicode[STIX]{x1D719}$ . The solid (respectively dashed) line marks the average value with respect to the phase (respectively distance) and the shaded area indicates the range of values reached across the various distances (respectively phases).

Figure 25. Snapshot of the vorticity field for two fish undulating at $f=1.5$ with separation distance $d=1$ and optimal phase $\unicode[STIX]{x0394}\unicode[STIX]{x1D719}=0.83$ at time $t/T-\unicode[STIX]{x1D719}=0.25~(\text{mod }1)$ . The colour axis is the same as in figure 13.

Figure 26. Snapshot of the velocity and pressure field for two fish undulating at $f=1.5$ with separation $d=1$ and optimal phase $\unicode[STIX]{x0394}\unicode[STIX]{x1D719}=0.83$ . (a,b) $x$ -velocity; (c,d) $y$ -velocity; (e,f) pressure and arrows showing the velocity of the fish. (a,c,e) Upstream fish at $t/T=0.25~(\text{mod }1)$ ; (b,d,f) downstream fish at $t/T-\unicode[STIX]{x1D719}=0.25~(\text{mod }1)$ . The same colour axis as in figure 14 is used for the pressure, and the same as in figure 22(b) is used for the velocity (centred in $0$ for the $y$ -velocity).

For the downstream fish, the situation is very different. Even when the upstream fish is several chord lengths ahead, the downstream fish encounters its wake. The performance of the downstream fish is determined by its interaction with the vortices of the wake. It appears from figure 24(b) that the swimming power of the downstream fish depends primarily on the phase difference $\unicode[STIX]{x0394}\unicode[STIX]{x1D719}$ between its undulation and the encountered vortices. Regardless of the distance $d$ , the swimming power of the downstream fish is low if $\unicode[STIX]{x0394}\unicode[STIX]{x1D719}$ is between 0.7 and 1, and it is high if $\unicode[STIX]{x0394}\unicode[STIX]{x1D719}$ is between 0.1 and 0.5. As for the upstream fish, the reduced swimming power is the result of a reduced undulation amplitude $a_{0}$ but, for the downstream fish, the reduction in $C_{P}$ is weaker than the reduction in amplitude. A maximum energy saving of 24 % is reached at $\unicode[STIX]{x0394}\unicode[STIX]{x1D719}=0.85$ , corresponding to an efficiency of $\unicode[STIX]{x1D702}_{QP}=65\,\%$ .

Figure 25 shows the vorticity field around the two fish undulating with frequency $f=1.5$ at distance $d=1$ for the phase $\unicode[STIX]{x0394}\unicode[STIX]{x1D719}=0.83$ , close to the minimum amplitude and power coefficient. At $t/T-\unicode[STIX]{x1D719}=0.25$ , the downstream fish approaches the negative vortex (at $x=-0.1$ on its left) as it is initiating the rotation of its ‘head’ (leading edge) to the right. This acceleration of the head causes a low pressure on the left side of the head, as shown in figure 26(e). Due to its position, the approaching negative vortex causes an increase in the longitudinal velocity, as shown in figure 26(b), which results in an increased stagnation pressure (figure 26 f). However, this vortex also generates a large transverse velocity with negative sign, as shown in figure 26(d). As a result, the effects of the head motion are amplified by the incoming vortex, displacing the stagnation point downstream on the right side and deepening the low pressure on the left side (figure 26 f). While the energy required by the fish to rotate its head is increased, the drag is decreased, despite the faster flow encountered by the fish.

At the same time, the positive vortex located at $x=0.2$ on the right side of the fish thickens the boundary layer and significantly accelerates the flow in a region where the fish undulation already accelerates it (figure 26 a,b). This interaction between the vortex and the fish results in a very large pressure drop around $x=0.3$ that also contributes to the reduction in drag while increasing the swimming power. The vortices are convected downstream at a speed which is substantially lower than the phase speed of the fish deformation. Further downstream, the distance between the vortices and the fish increases, and their interaction becomes weaker. At the trailing edge, the phase between the vortices and the fish motion is close to $\unicode[STIX]{x03C0}$ , such that the positive vortex reaches $x=1$ as the trailing edge of the fish is at its leftmost position. This vortex will be shed just downstream of the same sign vortex shed by the upstream fish. The resulting wake configuration is unstable and it takes several body lengths for the wake to reorganize into two pairs of opposite sign vortices per cycle.

The results presented in this section are consistent with the experimental results from thrust-producing rigid pitching foils in an in-line configuration (Boschitsch et al. Reference Boschitsch, Dewey and Smits2014). We found that for small separation distance the propulsive efficiency of the upstream fish is greatly improved thanks to drag reduction. We also showed that the efficiency of the downstream fish only weakly depends on the separation distance, the primary parameter being the phase difference between the wake from the upstream fish and the undulating motion of the downstream fish. If the undulation amplitude was fixed, the downstream fish would experience an increased drag and decreased power for $0\leqslant \unicode[STIX]{x0394}\unicode[STIX]{x1D719}\leqslant 0.5$ , whereas it would experience a decreased drag and increase power for $0.5\leqslant \unicode[STIX]{x0394}\unicode[STIX]{x1D719}\leqslant 1$ . For a self-propelled fish, the energetic benefits of a reduced amplitude resulting from a reduced drag overcome the power increase caused by the vortices.

4.3 Effect of frequency for two undulating fish in line

Next, we fix the distance between the two fish to $d=1$ and vary their undulation frequency. For frequencies $f=[1.5,1.8,2.1]$ , their optimized Gaussian envelope is used, for which the parameters are summarized in table 4. Figure 27 shows that most of the conclusions drawn in § 4.2 still hold as the undulation frequency is increased. While the upstream fish is mostly unaffected by the presence and phase of the downstream fish, the undulation amplitude of the downstream fish is largest for $0\leqslant \unicode[STIX]{x0394}\unicode[STIX]{x1D719}\leqslant 0.5$ and smallest for $0.5\leqslant \unicode[STIX]{x0394}\unicode[STIX]{x1D719}\leqslant 1$ . However, the exact value of the optimal phase depends on the frequency, and the correlation between amplitude and power coefficient becomes weaker as $f$ increases. For instance, at $f=2.1$ , the amplitude is minimum for $\unicode[STIX]{x0394}\unicode[STIX]{x1D719}=0.85$ , but the power coefficient is minimum for $\unicode[STIX]{x0394}\unicode[STIX]{x1D719}=0$ . Whereas with $f=1.5$ most phases result in an increased amplitude and power coefficient, with $f=1.8$ and $f=2.1$ , the amplitude and power coefficient of the downstream fish never exceed that of the upstream fish. Therefore, at these frequencies, it is always beneficial to swim in the wake of an undulating fish, despite the jet-like flow encountered in this area.

Figure 27. Ratio of (a) amplitude and (b) power coefficient as a function of frequency $f$ and phase $\unicode[STIX]{x0394}\unicode[STIX]{x1D719}$ for two fish swimming in line at distance $d=1$ . The ratios are defined with respect to the corresponding gait in open water. The black dots show the location of the points that have been used to build the thin-plate smoothing spline (tpaps function in Matlab with smoothing parameter $p=0.999$ ) represented in colour. Note that the thin-plate smoothing is only accurate in regions with a sufficient density of points.

At frequency $f=1.8$ , the results for the optimal phase, shown in figure 28, are very similar to those described in the previous section for $f=1.5$ , with the vortices from the upstream fish reinforcing the effects of the body undulation. However, the wake at this frequency is narrower; therefore the vortices are closer to the fish and they lose more strength through their interaction with the boundary layer. Moreover, since the distance between two consecutive vortices is proportional to the undulation frequency, the vortices are spaced closer to each other. The resulting wake is dominated by two single vortices shed by the downstream fish, each accompanied by weaker vortices of opposite sign from the upstream fish. This is identical to the destructive vortex interaction mode of a flapping foil within a Kármán street in Gopalkrishnan et al. (Reference Gopalkrishnan, Triantafyllou, Triantafyllou and Barrett1994), which has been associated with increased efficiency.

Figure 28. Snapshot of the vorticity, velocity and pressure field for two fish undulating at $f=1.8$ with separation distance $d=1$ and optimal phase $\unicode[STIX]{x0394}\unicode[STIX]{x1D719}=0.87$ . (a,b) Vorticity; (c,d) $x$ -velocity; (e,f) $y$ -velocity; (g,h) pressure and arrows showing the velocity of the fish. (a,c,e,g) Upstream fish at $t/T=0.25~(\text{mod }1)$ ; (b,d,f,h) downstream fish at $t/T-\unicode[STIX]{x1D719}=0.25~(\text{mod }1)$ . The same colour axis as in figure 14 is used for the pressure, and the same as in figures 22(a) and 22(b) for vorticity and velocity (centred in 0 for the $y$ -velocity).

Figure 29 illustrates the case of the downstream fish undulating with the phase providing the worse swimming efficiency for $f=1.8$ . In this configuration, the vortices from the upstream fish counteract the effects of the undulating motion of the downstream fish. As the fish turns its head to the right, displacing the stagnation point to the right and causing a low pressure on the left side of the head, the positive $y$ -velocity caused by the approaching positive vortex has the opposite effect. The high velocity caused by the vortices along the fish cancels the low velocity from the undulating motion. Finally, when the vortices from the upstream fish reach the trailing edge, they merge with the same sign vortices from the downstream fish. The resulting wake is very stable and is a classical reverse Kármán vortex street much wider than the one behind a single fish. This corresponds precisely to the constructive interaction mode of a flapping foil within a Kármán street in Gopalkrishnan et al. (Reference Gopalkrishnan, Triantafyllou, Triantafyllou and Barrett1994), which has been shown to have reduced efficiency. Note that here the organization of the vortices cannot be used to quantify the momentum in the wake, since the total momentum is zero for a self-propelled body, but can be used to quantify the energy shed in the wake.

Figure 29. Snapshots of (a) vorticity, (b) $x$ -velocity and (c) pressure field for two fish undulating at $f=1.8$ with separation distance $d=1$ and phase $\unicode[STIX]{x0394}\unicode[STIX]{x1D719}=0.38$ at $t/T-\unicode[STIX]{x1D719}=0.25~(\text{mod }1)$ . The same colour axis as in figure 22(a) is used for the vorticity and the same as figure 14 for the pressure.

As the frequency increases further, the vortices from the upstream fish lose even more energy through interaction with the boundary layer of the downstream fish; therefore for each period of oscillation the wake behind the two fish contains a pair of single vortices, as shown in figure 30. Moreover, since the efficiency of the downstream fish mostly depends on the phase of the leading edge with respect to the arrival of the upstream fish reverse Kármán street vortices, the phase of the trailing edge with respect to the incoming vortices in the optimal configuration changes with frequency.

Figure 30. Snapshot of the vorticity field for two fish undulating at $f=2.1$ with separation distance $d=1$ and phase $\unicode[STIX]{x0394}\unicode[STIX]{x1D719}=0$ . The same colour axis is used as in figure 22(a).

For all the frequencies considered here, a self-propelled fish can save energy by undulating behind another self-propelled fish undulating at the same frequency, reaching efficiencies close to $\unicode[STIX]{x1D702}_{QP}=70\,\%$ . When the motion of the downstream fish is phased appropriately with respect to the incoming vortex street, the vortices can reinforce the effect of the undulation. Whereas for a fixed amplitude this phase would result in an increased swimming power, the reduction in drag results in an overall decreased swimming power.

4.4 Fish undulating in the reduced velocity region of the wake

We have so far considered the case of a pair of fish swimming in an in-line configuration. Since our fish model has a feedback controller ensuring its stability in $y$ , it is also possible to impose an asymmetric configuration with an offset in the $y$ direction. Indeed, according to Weihs’ theory (Weihs Reference Weihs1973), the only way for a fish to save energy in a school is to swim in the region of reduced velocity located between two wakes. We have already shown that a fish can reduce its drag and save energy by swimming directly in the wake of another self-propelled fish and will now investigate if additional savings are possible by swimming in areas of reduced flow velocity. Figure 22(d) shows that, even with a single fish upstream, the flow on either side of the wake is slower than the free stream: for $f=1.5$ the average flow is slowest at $y=\pm 0.17$ . With the downstream fish offset from the upstream fish by $\unicode[STIX]{x0394}y=0.17$ , we vary the phase difference $\unicode[STIX]{x0394}\unicode[STIX]{x1D719}$ in order to see if the downstream fish can also save energy when swimming at this location.

Figure 31 shows that, even when the downstream fish is offset from the vortex street, its swimming performance greatly varies with the phase. However, it is easier for the fish to save energy in this region of reduced flow velocity than directly behind the upstream fish. Directly behind the upstream fish, only 30 % of the phases result in energy savings, and by using the unsteadiness of the wake, the quasi-propulsive efficiency at $f=1.5$ can be brought up from 50 % to 60 %. When undulating in the region of reduced flow velocity, it is easier to save energy since over 70 % of the phases result in energy savings. The energy savings can even be very large since $\unicode[STIX]{x1D702}_{QP}=80\,\%$ is possible for $\unicode[STIX]{x0394}\unicode[STIX]{x1D719}=0.65$ .

Figure 31. Ratio of undulation amplitude $a_{0}$ and time-averaged power coefficient $\overline{C_{P}}$ as a function of phase for two fish undulating at $f=1.5$ with longitudinal separation distance $d=1$ . In-line fish and fish at an offset $\unicode[STIX]{x0394}y=0.17$ are compared.

Figure 32 shows that, at the optimal phase, the leading edge of the downstream fish reaches its leftmost position at the same time as it reaches a positive vortex. Figure 33(b) shows that the leading edge of the downstream fish exactly passes through the region where the longitudinal velocity is smallest. As a result, the stagnation pressure is greatly reduced (figure 33 d). Moreover, the region of accelerated flow between the negative vortex and the fish ( $x=0.4$ ) reinforces the accelerated region caused by the undulation, which we showed earlier is beneficial.

Figure 32. Snapshot of the vorticity field for two fish undulating at $f=1.5$ with longitudinal separation distance $d=1$ , transverse separation $\unicode[STIX]{x0394}y=0.17$ and optimal phase $\unicode[STIX]{x0394}\unicode[STIX]{x1D719}=0.65$ at time $t/T-\unicode[STIX]{x1D719}=0.1~(\text{mod }1)$ . The colour axis is the same as in figure 13.

Figure 33. Snapshot of the (a,b) $x$ -velocity and (c,d) pressure field for two fish undulating at $f=1.5$ with longitudinal separation distance $d=1$ , transverse separation $dy=0.17$ and optimal phase $\unicode[STIX]{x0394}\unicode[STIX]{x1D719}=0.65$ . (a,c) Upstream fish, $t/T=0.1~(\text{mod }1)$ ; (b,d) downstream fish, $t/T-\unicode[STIX]{x1D719}=0.1~(\text{mod }1)$ . The same colour axis as in figure 14 is used for the pressure, and the same as in figure 22(b) for the velocity.