# Far-field hydrodynamic interaction in a group of swimmers

Gen Li, Lei Duan, Jörn Sesterhenn, Ramiro Godoy-Diana, Benjamin Thiria, Dmitry Kolomenskiy

► To cite this version:

Gen Li, Lei Duan, Jörn Sesterhenn, Ramiro Godoy-Diana, Benjamin Thiria, et al. Far-field hydrodynamic interaction in a group of swimmers. Journal of Fluid Mechanics, 2023, 974, pp.A34. 10.1017/jfm.2023.802. hal-04285008

HAL Id: hal-04285008

https://hal.science/hal-04285008v1

Submitted on 14 Nov 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
<figure><img src="image_2.png" /></figure>

# Far-field hydrodynamic interaction in a group of swimmers

Gen Li<fnref n="1" />, Lei Duan<fnref n="2" />, Jörn Sesterhenn<fnref n="3" />, Ramiro Godoy-Diana<fnref n="4" />, Benjamin Thiria<fnref n="4" /> and Dmitry Kolomenskiy<fnref n="5" />,†

(Received 27 January 2023; revised 25 August 2023; accepted 12 September 2023)

Three-dimensional schools of hydrodynamically axisymmetric swimmers self-propelling at a constant velocity are studied. We introduce a low-order model for the induced velocity based on the far-field approximation. We inquire if, by holding suitable relative positions in the three-dimensional space, the swimmers can reduce the overall energy consumption of the school in comparison with the same number of isolated individuals at the same velocity. We find a considerable (several per cent) energy saving achievable by chain formations. The benefit increases asymptotically with the number of individuals, towards a finite limit that is a function of the minimum allowed spacing between each pair of neighbours.

**Key words:** collective behaviour, propulsion, swimming/flying

## 1. Introduction

Living organisms often move as a group. Besides social interactions that drive such collective locomotion, fluid-dynamic coupling has been recognized as an important factor for formation and gait selection in groups of swimming and flying organisms (Herskin & Steffensen 1998; Weimerskirch et al. 2001; Brumley et al. 2014; Ashraf et al. 2016;

†Email address for correspondence: dkolom@gmail.com

<footnote n="1">¹Center for Mathematical Science and Advanced Technology, Japan Agency for Marine-Earth Science and Technology (JAMSTEC), Yokohama 236-0001, Japan</footnote>

<footnote n="2">²School of Naval Architecture, Ocean and Civil Engineering, Shanghai Jiao Tong University, Shanghai 200240, PR China</footnote>

<footnote n="3">³Lehrstuhl für Technische Mechanik und Strömungsmechanik, Universität Bayreuth, 95447 Bayreuth, Germany</footnote>

<footnote n="4">⁴Laboratoire de Physique et Mécanique des Milieux Hétérogènes (PMMH), CNRS UMR 7636, ESPCI Paris-PSL University, Sorbonne Université, Université Paris Cité, 75005 Paris, France</footnote>

<footnote n="5">⁵Center for Materials Technologies (CMT), Skolkovo Institute of Science and Technology, Moscow 121205, Russia</footnote>
Li *et al.* 2021; Yuan *et al.* 2021). It can reduce the metabolic cost in birds and fish alike (Herskin & Steffensen 1998; Weimerskirch *et al.* 2001). However, the energy harvesting mechanisms differ. Birds may acquire extra lift by surfing updraft induced by the leader's trailing vortices. But neutrally buoyant underwater swimmers (biological and artificial alike) spend no energy to support their own body weight, therefore they can save energy only by reducing body drag and propelling forwards more efficiently.

Important theoretical results on the hydrodynamics of fish schooling were obtained using the two-dimensional (2-D) approximation. Weihs (1975) considered a planar layer of a special kind of three-dimensional (3-D) school. He postulated that the latter consisted of a large number of identical layers that could be homogenized in the direction perpendicular to the selected plane. Then he applied a point vortex model for the wake and derived the diamond formation. More recently, Gazzola *et al.* (2016) and Filella *et al.* (2018) combined the hydrodynamic interaction and the behavioural models. The velocity field induced by an individual fish was approximated by that of a vortex dipole. The dipole strength was related to the swimming speed using a self-propelled dipole model (Kanso & Tsang 2014). Thus hydrodynamic interactions between individuals reduced to superposition of the induced flow velocity.

While the vortex dipole model offered a computationally tractable alternative to the numerical solution of the Navier–Stokes equations, a recent study by Zhang, Peterson & Porfiri (2023) pinpointed its deficiencies, such as the absence of a wake structure and the poor approximation when used for elongated bodies. The latter problem was solved in the same study by redistributing the dipole vortex circulations over two sheets of vorticity oriented along the body. However, the amended model remained a 2-D approximation.

The current exploration of energy-saving mechanisms in three dimensions, on the other hand, focuses on near-field interactions in a pair of fish: the leader and the follower. The follower smartly adjusts its position and tail beat phase to benefit from the unsteady flow perturbation without penalizing the leader (Li *et al.* 2020; Seo & Mittal 2022). Obviously, any large school can be split into pairs of leaders and followers that interact through the near field. But the question remains open whether the many-to-many hydrodynamic interaction through the far field can be of any significance. It follows from the Biot–Savart formula that the induced velocity of a bounded region of vorticity decays more slowly in the 2-D space than in the 3-D space. For that reason, energy-saving estimates obtained from 2-D models may be overly optimistic. In fact, one may even doubt if there can be any appreciable far-field interference in a 3-D school. On the other hand, a self-propelled body in the inertia-dominated flow regime produces a momentumless vortical wake that persists over a much longer distance than the size of the body and may be harvested by fellow individuals. Hereby, we present a far-field superposition model that accounts for the essential 3-D traits of the induced flow of individual self-propelled swimmers. By using the external data for the correlation between the induced flow strength with the individual mechanical power expenditure, the model provides quantitative estimates of the mechanical power efficiency of a school. To keep the model simple, we restrict our scope to axisymmetric individual wakes (of e.g. squid or rigid-body artificial underwater vehicles, rather than fish that produce lateral waves) and the steady-state approximation.

From the standpoint of minimizing the energy cost of one selected swimmer in a group, it is obvious that the swimmer should place itself in the low-speed region of its neighbour, because the hydrodynamic power increases with the perceived inflow velocity. But this behaviour may place an extra energy cost burden on the neighbours. It is not self-evident if there exist school arrangements that minimize the hydrodynamic power of the entire group by exploiting the far-field wake interaction in three dimensions. In this paper, we
show constructively that they exist. We develop a model of such optimal schools, and provide quantitative estimates of energy savings depending on the minimum allowed distance between swimmers. The latter parameter needs to be determined by external considerations. The model that we construct is minimalist: the swimmers are substituted by a linear superposition of their steady-state far-field asymptotic representations. Yet it is sufficient to conclude that energy saving is possible without unsteady wake capture or phase matching, and the benefit accumulates with the number of swimmers in the school. The model is described in § 2. The optimization results are presented and discussed in § 3. Section 4 contains concluding remarks.

## 2. Mathematical formulation of the far-field interaction model

### 2.1. Velocity field of the bound vorticity

Let us first consider a circular vortex ring with an infinitesimally thin core. It is convenient to use a cylindrical polar coordinate system with the origin at the centre of the ring and the x-axis being the symmetry axis of the ring. Let $r$ denote the radial coordinate. The azimuthal coordinate does not appear in the equations. The radius of the vortex ring is $a$, and its strength (circulation) is $\Gamma$.

The Stokes stream function of the ring can be expressed in terms of complete elliptic integrals $K$ and $E$ (Batchelor 2000),

$$ \psi(x, r) = \frac{\Gamma\sqrt{ar}}{2\pi} \left\{ \left(\frac{2}{k} - k\right) K(k^2) - \frac{2}{k} E(k^2) \right\}, \quad \text{where } k = \left( \frac{4ar}{x^2 + (r+a)^2} \right)^{1/2} \qquad (2.1) $$

is the elliptic modulus, and

$$ K(k^2) = \int_0^1 \frac{dt}{\sqrt{(1-t^2)(1-k^2t^2)}}, \quad E(k^2) = \int_0^1 \frac{\sqrt{1-k^2t^2}}{\sqrt{1-t^2}} dt. \qquad (2.2a,b) $$

The axial velocity $u = r^{-1} \partial\psi/\partial r$ can also be reduced to the elliptic integrals using the identities

$$ \frac{dK(k^2)}{d(k^2)} = \frac{E(k^2) - (1-k^2)K(k^2)}{2(1-k^2)k^2}, \quad \frac{dE(k^2)}{d(k^2)} = \frac{E(k^2) - K(k^2)}{2k^2}, \qquad (2.3a,b) $$

yielding

$$ u = \Gamma\Upsilon_{\text{ring}}, \qquad (2.4) $$

with

$$ \Upsilon_{\text{ring}}(x, r) = \frac{a^{1/2}k}{4\pi r^{3/2}} \left\{ \frac{1}{1-k^2} \left( 1 - (2-k^2) \frac{r+a}{2a} \right) E(k^2) + \frac{r}{a} K(k^2) \right\}. \qquad (2.5) $$

When $r$ is small, numerical evaluation of (2.4) becomes problematic. The formula for $u_{\text{ring}}$ at the axis is derived by keeping only the leading-order terms $E(k^2) \approx \pi/2$ and $K(k^2) \approx \pi/2$ in the respective series expansions. The axial velocity component at $r = 0$ becomes equal to $\Gamma\Upsilon_{\text{ring}}(x, 0)$, with

$$ \Upsilon_{\text{ring}}(x, 0) = \frac{1}{2} \frac{a^2}{(x^2 + a^2)^{3/2}}. \qquad (2.6) $$

At lateral distance $r$ much larger than the swimmer's body length $l$, any axisymmetric steady swimmer induces the velocity asymptotically matching a vortex ring (note that the
<figure><img src="image_3.png" /><figcaption>Figure 1. (a) Schematic representation of the vortex ring and momentumless wake model. Dimensional circulations $\Gamma_j$ are related to the dimensionless coefficients $\gamma_j$ by $\Gamma_j = \gamma_j \Gamma$. (b) Postulated bound circulation distribution. (c) Postulated normalized radial profile of the momentumless wake axial velocity component. (d) Axial velocity component $u$ over the azimuthal symmetry plane, obtained using the far-field approximation. The black dashed line connects the points of minimum $u$ at every $d_{min}$ in the rear part of the field. The black dash-dotted line connects the points of minimum $u$ at every $d_{min}$ in the front part of the field.</figcaption></figure>

momentumless wake velocity decays exponentially with $r$; see § 2.2). Let $\Gamma$ now stand for the circulation of that ring. For slender bodies, however, it will be practical to consider an intermediate range of distances greater than the body half-width $a$, but possibly smaller than $l$. This velocity field can be approximated as induced by a finite number of elementary vortex rings, each substituting for the bound vorticity over a small portion of the body (see figure 1a):

$$u = \Gamma \Upsilon_{\text{bound}}, \quad \text{where} \quad \Upsilon_{\text{bound}} = \sum_{j=1}^{n_r} \gamma_j \Upsilon_{\text{ring}}(x - x_j, r), \qquad (2.7)$$

where $x_j = x_{\text{nose}} + (j - 1/2)l/n_r$ is the position of the $j$th vortex ring, and $x_{\text{nose}}$ corresponds to the front point of the body. The circulation distribution coefficients $\gamma_j$ are calculated as follows. First, dimensional circulations $\Gamma_j$ are determined from a least mean square fit of the velocity profile along a line parallel to the axis and at a constant distance from it. Then $\Gamma$ is found by summing all vortex ring circulations $\Gamma_j$. Finally, the coefficients are evaluated as $\gamma_j = \Gamma_j/\Gamma$.

The above approach can be regarded as a variant of the vortex lattice method, which substitutes the solid body in a potential flow by a perpendicular lattice of vortices. This approximation is also suitable for real flows if the boundary layers are thin and there is no flow separation. In the axisymmetric case, the longitudinal vortices have zero circulation. This leads us to the approximation of the body by a system of coaxial vortex rings. If we calculate the velocity at a test point so far from the body that the distances between the vortex rings are much less than the distance from the body to the test point, then it becomes possible to neglect the distance between the coaxial vortex rings, and associate all vorticity with only one vortex ring of radius $a$ and circulation $\Gamma$.
## 2.2. Velocity field of a momentumless wake

An analytical solution of the incompressible Navier-Stokes equations that describes a steady laminar momentumless wake can be found in Korobko, Shashmin & Shul'man (1986). However, that solution shows the velocity decay downstream as $x^{-2}$, whereas turbulent momentumless wakes rather decay as $x^{-1.5}$ (Sirviente & Patel 2000). Since there is no analytical solution for the mean flow velocity of the turbulent momentumless wake, let us modify the laminar functional relation (Korobko et al. 1986) to become consistent with the turbulent rate of decay. For this purpose, we introduce adjustable parameters C and p. Also, we regularize the singularity at $x = 0$. Thus the axial velocity component is given by

$$u = \Gamma \Upsilon_{\text{wake}}, \quad \Upsilon_{\text{wake}} = \frac{C}{a} \xi \frac{1}{1 + \xi^{2-2p}} g(\eta), \qquad (2.8a,b)$$

where it is postulated that the radial profile is a function of the similarity variable

$$\eta = \frac{r}{a} \frac{\xi^p}{Re^p}, \qquad (2.9)$$

with $\xi = x/a$, and the Reynolds number $Re = U_\infty a/\nu$ being based on $U_\infty$, the far-field velocity that is equal in magnitude to the swimming speed, and $\nu$, the kinematic viscosity of the fluid. Since this is a far-field approximation, the virtual wake origin is set at $x = 0$. This simplification has little effect on the accuracy when $\xi \gg l/a$.

The radial profile of the velocity, $g(\eta)$, depends on the spatial distribution of the momentum excess and deficit in the wake. We are interested particularly in the case of $g(\eta)$ changing sign twice (see Appendix B for more detail):

$$g(\eta) = - \left(1 - \frac{\eta^2}{\eta_1^2}\right) \left(1 - \frac{\eta^2}{\eta_2^2}\right) \exp\left(-\frac{\eta^2}{\eta_3^2}\right), \qquad (2.10)$$

where $\eta_1$, $\eta_2$ and $\eta_3$ are constant parameters.

Since the wake is momentumless and axisymmetric, it is required that

$$\int_{0}^{\infty} 2\pi\eta g(\eta) d\eta = 0. \qquad (2.11)$$

It follows that only two of the three parameters can be adjusted freely, since

$$2 \frac{\eta_3^2}{\eta_1^2} \frac{\eta_3^2}{\eta_2^2} - \frac{\eta_3^2}{\eta_1^2} - \frac{\eta_3^2}{\eta_2^2} + 1 = 0. \qquad (2.12)$$

The far-field approximation to the velocity due to both the bound vorticity and the momentumless wake is constructed as

$$u = \Gamma \Upsilon, \quad \text{where} \quad \Upsilon = \begin{cases} \Upsilon_{\text{bound}} & \text{if } x < 0, \\ \Upsilon_{\text{bound}} + \Upsilon_{\text{wake}} & \text{otherwise.} \end{cases} \qquad (2.13)$$

The postulate that the momentumless wake intensity scales linearly with the bound circulation $\Gamma$ is justified as follows. The drag region of the wake is formed by the separated boundary layer, which carries the bound circulation. The jet region has momentum equal in magnitude but of opposite sign, to ensure that the total momentum vanishes. Therefore, the wake velocity must be proportional to $\Gamma$.
## 2.3. Condition of multiple swimmers propelling at the same velocity

Let us focus on a situation when N swimmers self-propel in a fluid otherwise at rest. All swim at the same speed $-U_0$ in the same direction, such that their relative positions do not vary. Hence, in the reference frame moving with the swimmers, the velocity field of the fluid is steady. The inflow velocity at the location of the ith swimmer is the superposition of $U_0$ and the velocity induced by all companion swimmers:

$$U_i = U_0 + \sum_{\substack{j=1 \\ j \neq i}}^{N} u_{ij}, \qquad (2.14)$$

where $u_{ij}$ is the velocity induced by the jth swimmer at the location of the ith swimmer $x_i$. These inflow locations are taken at the swimmer body centres, assuming that the induced velocity has little variation over the body length. Let us use the far-field approximation derived in the previous subsections,

$$u_{ij} = \Gamma_j \Upsilon(x_i), \qquad (2.15)$$

where $\Upsilon(x_i)$ is calculated as defined in (2.13).

The circulation $\Gamma$ and the hydrodynamic power $P$ of a swimmer are functions of the local inflow velocity of the surrounding fluid $U$. These functions can be linearized for a small deviation of $U$ from the reference velocity $U_0$:

$$\Gamma(U) = \Gamma|_{U_0} + \left. \frac{\partial \Gamma}{\partial U} \right|_{U_0} (U - U_0), \qquad (2.16)$$

$$P(U) = P|_{U_0} + \left. \frac{\partial P}{\partial U} \right|_{U_0} (U - U_0), \qquad (2.17)$$

where $\Gamma$, $P$ and their derivatives at $U_0$ are evaluated using computational fluid dynamics (CFD) simulations, as explained below. For a swimmer surrounded by neighbours, it is assumed that $U$ differs from $U_0$ by a small amount of induced velocity, warranted for the Taylor series approximation (2.16) and (2.17). Figure 1(d) confirms that $|u_{ij}| < 0.02U_0$ for $y > 2a$. Substituting (2.14) and (2.15) into (2.16), we obtain

$$\Gamma_i = \Gamma|_{U_0} + \left. \frac{\partial \Gamma}{\partial U} \right|_{U_0} \sum_{\substack{j=1 \\ j \neq i}}^{N} \Gamma_j \Upsilon(x_i), \quad i = 1, \dots, N. \qquad (2.18)$$

When positions $x_1, \dots, x_N$ are fixed a priori, (2.18) constitutes a system of $N$ linear equations with $N$ unknowns $\Gamma_1, \dots, \Gamma_N$. By solving this system and substituting the solution into (2.15), we evaluate the induced velocities $u_{ij}$ for all $i = 1, \dots, N$ and $j = 1, \dots, N$. Then the effective inflow velocities $U_i$ of each swimmer are evaluated using (2.14). Finally, using (2.17), we calculate the mechanical power expenditure of each swimmer as

$$P_i = P|_{U_0} + \left. \frac{\partial P}{\partial U} \right|_{U_0} (U_i - U_0), \quad i = 1, \dots, N. \qquad (2.19)$$
## 2.4. Optimization objective

We aim to minimize the normalized school-average mechanical power expenditure,

$$ \bar{P} = \frac{1}{P|_{U_0} N} \sum_{i=1}^{N} P_i, \qquad (2.20) $$

under the condition that all swimmers move in the same direction at the same constant speed $U_0$. Minimization is achieved by finding the best relative position of the swimmers. Thus, without loss of generality, we assume $x_1 = 0$. The Cartesian scalar components of $x_2, \dots, x_N$ provide $3N - 3$ optimization parameters.

The far-field approximation fails in the vicinity of the swimmer. Therefore, only solutions with large enough spacing between the swimmers should be admitted as optimal. This is achieved by introducing a penalty function

$$ V_p = \max \left\{ 0, A_p \left( 1 - \frac{d_{min}}{\sigma_p a} \right) \right\}, \qquad (2.21) $$

where $\sigma_p$ is a parameter that controls the minimum allowed distance between swimmers as a multiple of $a$, and $A_p \gg 1$ is a constant penalty parameter (we use $A_p = 10^4$).

The pairwise lateral separation distance between the swimmer centroids $x_i = (x_i, y_i, z_i)^T$, $i = 1, \dots, N$, is calculated as

$$ d_{min} = \min_{\substack{i=1,\dots,N \\ j=1,\dots,i-1}} \sqrt{(y_j - y_i)^2 + (z_j - z_i)^2}. \qquad (2.22) $$

The penalty function (2.21) ensures that the lateral separation between swimmers is no less than $d_{min}$. The axial coordinate is not included in the formula, therefore the optimizer cannot choose to place swimmers either in the immediate proximity or in the wake. Thus we disallow placing swimmers at locations where the flow is vorticity-dominated (as known from the CFD data). The minimization objective function is

$$ J = \bar{P} + V_p. \qquad (2.23) $$

# 3. Results of the numerical optimization

## 3.1. Individual swimmer parameter fitting

The far-field interaction model requires knowledge of the velocity field of an isolated swimmer, as well as its hydrodynamic power expenditure. We derived representative values of the necessary parameters from CFD simulations of a real-scale streamlined-body artificial swimmer using the SIEMENS PLM software Simcenter STAR-CCM+ 2020.1. The CFD methodology is described in a precursor study by Li et al. (2023), and the simulations used in this study are detailed in the Supplementary Information. The swimmer has body half-width $a = 0.175$ m and body length (BL) $l = 20a$. In all calculations, we place the origin of the attached coordinate system at a distance $9.4a$ downstream from the nose. Thus the position of the tail in this coordinate system is $x_{tail} = 10.6a$. The swimmer cruise speed is $U_0 = 0.466$ BL s$^{-1}$. The numerical dimensional output data from the CFD simulations are converted in this section into dimensionless coefficients multiplied by proper scaling factors, for greater generality. Thus the hydrodynamic power at cruise is evaluated in STAR-CCM+ as $P|_{U_0} = 0.107 \rho U_0^3 \pi a^2$, where $\rho$ is the water density.
The total circulation $\Gamma$ of the far-field vortex ring model was determined as follows. The CFD velocity was sampled at the grid points distributed on a line that was offset from the symmetry axis by 2.27a. Then the far-field formula (2.7) was evaluated at the same points. The bound vorticity was approximated by $n_r = 21$ rings. The products $G_j = \Gamma \gamma_j$ were treated as unknown coefficients, thus their values were found using the linear regression. The total circulation $\Gamma|_{U_0} = -16.5U_0a$ was found by the summation of $G_j$ over $j$ from 1 to $n_r$. Finally, the distribution coefficients $\gamma_j = G_j/\Gamma$ were determined; see figure 1(b).

Further, the momentumless wake model coefficients were fitted to the CFD data at $U_0$. The corresponding Reynolds number was $Re = aU_0/v = 2.8 \times 10^5$. The radial profile $g(\eta)$ was adjusted to the data at distance 3.26a downstream from the tail. This resulted in the parameter values $\eta_1 = 2.15$, $\eta_2 = 5.5$, $\eta_3 = 3.34$, and a profile shown in figure 1(c). The parameter values $C = 0.516$ and $p = -0.176$ ensure that the decay along the wake axis agrees with the CFD data. Figure 1(d) visualizes the velocity reconstructed from the far-field approximation.

The velocity field of the momentumless wake decays in the lateral direction $r$ much faster (exponentially) than the velocity field of the vortex rings (that decays algebraically). Therefore, when the velocity is sampled on a line parallel to the axis, if the offset distance is sufficiently large for the far-field vortex ring approximation to be valid, then the momentumless wake contribution at that station is negligible. Likewise, on the axis behind the swimmer, the momentumless wake velocity decays more slowly with the downstream distance $x$ than the vortex-ring-induced velocity. Therefore, coefficients of the bound vorticity model and coefficients of the momentumless wake model can be determined independently.

Calculation of the derivatives $(\partial\Gamma/\partial U)|_{U_0}$ and $(\partial P/\partial U)|_{U_0}$ required additional simulations at two slightly different speeds. These simulations used a coarser mesh. At $U_-=0.9U_0$, we obtained the circulation $\Gamma_- = -15.9U_0a$ and power $P_- = 0.079\rho U_0^3 \pi a^2$. At $U_+ = 1.1U_0$, $\Gamma_+ = -19.4U_0a$ and $P_+ = 0.140\rho U_0^3 \pi a^2$. These data enable finite-difference approximations to the derivatives, $(\partial\Gamma/\partial U)|_{U_0} = -17.4a$ and $(\partial P/\partial U)|_{U_0} = 0.304\rho U_0^2 \pi a^2$.

## 3.2. Parameters of optimal schools

The optimization problem consists in finding the swimmer positions $x_2, \dots, x_N$ relative to the first swimmer such as to minimize the objective function (2.23). We use the CMA-ES genetic optimization algorithm (Hansen, Müller & Koumoutsakos 2003). To accelerate the search, the swimmer coordinates vary with discrete increments that are multiples of $0.01a$, and the whole domain is a cube with side length $2Nl$. In all cases, the population size parameter was equal to 500, and a sufficient number of iterations were performed to ensure that the objective function does not change by more than $10^{-4}$ during the last 100 iterations.

Since $u$ decays with the lateral spacing between the swimmers, an important control parameter is $\sigma_p$, which enters in the penalty function (2.21). All optimal formations that we found activated the constraint $d_{min} = \sigma_p a$. Figure 2(a) shows how the school average power $\bar{P}_{min}$ that corresponds to the minimized cost function $J_{min}$ varies with the minimum dimensionless separation distance between swimmers $d_{min}/a$, which is controlled by $\sigma_p$. The latter is an external control parameter. Its value may be selected, for example, from the accuracy considerations of the far-field approximation or from collision avoidance considerations. The relationship between $d_{min}/a$ and $\bar{P}_{min}$ suggests that reducing the lateral distance in a school can enhance energy-saving effects. Due to this, $d_{min}/a$ would tend to
Far-field hydrodynamic interaction in a group of swimmers

<figure><img src="image_4.png" /><figcaption>Figure 2. (a) School average power $\bar{P}_{min}$ as a function of the minimum dimensionless lateral separation distance between swimmers, $d_{min}/a$, for four different school sizes $N$. Continuous lines correspond to the results of simultaneous optimization of the positions of $N-1$ swimmers. Markers correspond to the results obtained by sequentially adding one swimmer to the school and optimizing its position, shown for $d_{min}/a = 2, 3$ and $4$ as ‘○’, ‘×’ and ‘+’ signs, respectively. (b) The quantity $1 - \bar{P}_{min}$ shown as a function of $d_{min}/a$ on the logarithmic scale, for different school sizes $N$. The continuous lines show the optimization results, and the dashed lines are the linear regression lines on those data: $1 - \bar{P}_{min} = 0.220(d_{min}/a)^{-1.97}$ for $N = 2$, $1 - \bar{P}_{min} = 0.286(d_{min}/a)^{-1.90}$ for $N = 3$, $1 - \bar{P}_{min} = 0.354(d_{min}/a)^{-1.82}$ for $N = 7$, and $1 - \bar{P}_{min} = 0.392(d_{min}/a)^{-1.83}$ for $N = 11$. (c) Plots of $\bar{P}_{min}$ obtained using the sequential increment of $N$, as functions of $N$, for $d_{min}/a = 3, 4$ and $5$.</figcaption></figure>

0 under unrestricted $d_{min}/a$. In the context of optimization, the role of the penalty function is to impose constraints on the dimensionless lateral separation. It is necessary to impose such constraints on $d_{min}/a$ during the optimization process to avoid swimmer overlap and to match real-world situations. Energy savings of several per cent (relative to $P|_{U_0}$) may be realized at $d_{min}$ as small as $2a$ to $5a$. Greater lateral separation leads to vanishing energy saving. Trends for schools of different sizes are similar, but lower $\bar{P}_{min}$ of larger schools shows that they are slightly more efficient in saving energy. To elucidate the trend with $d_{min}/a$, figure 2(b) shows the power economy $1 - \bar{P}_{min}$ on the logarithmic scale. It decays at a rate close to $-2$ for all $N$, and the power-law constant is larger for greater $N$.

Figure 3 shows an example of a typical optimal school. The frontal projections occupy a small region of the yz plane, such that the lateral distances between neighbours are almost as small as possible, i.e. about $5a$ (0.25 BL). From the side view, each pair of neighbours forms a leader–follower pair with $11.4a$ (0.57 BL) longitudinal offset. Incidentally, the best energy-saving formations of two tetra fish were previously found to have respective separation distances of about 0.35 BL and 0.75 BL (Li et al. 2019), and the optimum was explained from the mean flow field considerations. It appears that similar effects apply to an elementary leader–follower pair in a school of axisymmetric swimmers. Note that
<figure><img src="image_5.png" /><figcaption>Figure 3. Optimized schools of $N = 11$ individuals with $d_{min}/a = 5$: (a) front view and (b) side view of a school obtained by numerical optimization; (c) front view and (d) side view of a school obtained by a direct-rule-based method. Individuals are labelled by numbers from 1 to 11.</figcaption></figure>

any swimmer has at most two neighbours in its vicinity. Therefore, the induced velocities remain small.

### 3.3. Two- and three-swimmer optimal schools

The optimality of chain formation motivates us to revisit the groups of two and three swimmers. They can be considered as elementary blocks of larger schools, because the nearest neighbours interact most strongly, even if in our model they rely only on the far field such that unsteady effects are not included.

The minimal school consists of two swimmers ($N = 2$); see figure 4(a). The position of the second swimmer relative to the first swimmer is optimized. Without loss of generality, the first swimmer is located at the origin of the coordinate system, and the second swimmer is on the *xy* plane. Thus only the distance $\rho$ and the polar angle $\theta$ are non-zero. The optimal values are shown in figure 5, as functions of the minimum dimensionless lateral separation distance $d_{min}/a$. Both quantities increase monotonically. This means that the optimally selected lateral separation increases at a greater rate than the longitudinal separation.

As conjectured in § 1, the follower is optimally placed in the rear updraft region of the leader's induced flow: the blue circles in figures 5(a,b), which show the position of the follower relative to the swimmer, almost coincide with the black dashed lines that visualize the leader's induced field minimum velocity point coordinates as functions of relative separation distance $d_{min}/a$. Importantly, this formation also minimizes the leader's energy cost, because the leader's position relative to the follower is in the front updraft region of the follower's induced flow: the blue circles are also very close to the black dash-dotted lines.

Similar staggered formations with lateral separation 0.5l and longitudinal offset 0.75l were found to minimize the cost of transport of a pair of tetra fish (Li *et al.* 2019). This arrangement corresponds to $d_{min}/a = 6$, $\rho/a = 10.8$, $\theta = 34^\circ$ in polar coordinates. It was explained in Li *et al.* (2019) by the same mechanism of mean velocity updraft as considered in this study. However, in Li *et al.* (2019), the fish undulated laterally, and the amount of power economy varied with the phase difference between the leader and the follower. The greatest power economy would be achieved at the phase matching of the
Far-field hydrodynamic interaction in a group of swimmers

<figure><img src="image_6.png" /><figcaption>Figure 4. Geometrical parameters of (a) a two-swimmer school, and (b) a three-swimmer school. Points $S_1$, $S_2$ and $S_3$ correspond to the first, second and third swimmers, respectively. The coordinate system $S_2x'y'z'$ is obtained from $S_1xyz$ by translation of the origin from the first swimmer to the second swimmer. The planes $S_2x'y'$ and $S_1xy$ coincide, because the z-coordinate of the second swimmer is zero. The angle $\theta_{lead}$ is between the negative semi-axis $S_1x$ and the line that connects $S_1$ with $S_2$. Point $O''$ marks the perpendicular projection of the third swimmer on the $S_2x'$ axis. Axis $O''y''$ is parallel to $S_2y'$, and it resides on the $S_2x'y'$ plane. Angle $\theta_{hind}$ is the angle between the negative semi-axis $S_2x'$ and the line that connects $S_2$ with $S_3$. Angle $\phi_{hind}$ is the angle between the line $O''S_3$ and the $O''y''$ axis. Schematic side views: (c) a two-swimmer school; (d) a three-swimmer school. Respective top views: (e) a two-swimmer school; (f) a three-swimmer school. The minimum dimensionless lateral distance between neighbours is $d_{min}/a = 3$.</figcaption></figure>

follower undulation with the incoming wave from the leader, such as that described by Li *et al.* (2020). Thus benefit from the far-field interaction and unsteady mechanisms can accumulate.

In a school of three swimmers ($N = 3$, figure 4b), the optimal configuration of the leading pair is close to that found in the two-swimmer case ($N = 2$): $\rho_{lead}$ is slightly less and $\theta_{lead}$ is slightly larger than the respective values for $N = 2$. The points fall on the rear updraft line. The position of the third swimmer relative to the second swimmer is given by three non-zero polar coordinates. The values $\rho_{hind}$ and $\theta_{hind}$ coincide with $\rho_{lead}$ and $\theta_{lead}$, respectively. The $\phi_{hind}$ angle is 120°. The two pairs are, essentially, two two-swimmer optimal formations combined and rotated such that the lateral separation distance between the first and the third swimmer equals $d_{min}$. Thus the follower in each pair maximizes the beneficial interaction with its leader, and the two pairs adjust their relative orientation to gain extra benefit.

3.4. Rules for large schools

The previous three-swimmer analysis motivates a direct-rule-based method for constructing larger energy-saving schools. The first three swimmers are placed such that $\rho_{lead}$, $\theta_{lead}$, $\rho_{hind}$ and $\theta_{hind}$ are on the rear updraft line, as in figures 5(a,b) at the required
<figure><img src="image_7.png" /><figcaption>Figure 5. Optimal values of the geometrical parameters of two- and three-swimmer schools: (a) dimensionless distance between neighbours; (b) elevation; (c) bank. The black dashed lines show the polar coordinates ρ and θ of the minimum (most negative) induced velocity in the rear part of the swimmer's field, and they correspond to the black dashed line in figure 1(d). The black dash-dotted lines show the polar coordinates ρ and −θ of the minimum (most negative) induced velocity in the front part of the swimmer's field, and they correspond to the black dash-dotted line in figure 1(d).</figcaption></figure>

$d_{min}$, and $\phi_{hind} = 120^{\circ}$. The fourth swimmer is added behind the third one such that the second, third and fourth swimmers form a three-swimmer school obeying the same rule, but the $\phi$ angle is $-120^{\circ}$. Thus we ensure that the fourth swimmer is in the updraft of all those in front. The fifth swimmer is added with $\phi = 120^{\circ}$, and so on. Figures 3(c,d) show a formation obtained for $N = 11$, and the pattern extends intuitively to any large $N$.

The results in figure 2(a) suggest that $\bar{P}_{min}$ is smaller for larger schools. To estimate the lower bound of $\bar{P}_{min}$, for a given $d_{min}/a$, we use the above rules to construct schools with $N$ ranging from 3 to 100. The markers in figure 5 show that the schools obtained using this method are close to optimal. Visually, the values of $\bar{P}_{min}$ are indistinguishable from the (3$N$ − 3)-parameter numerical optimization results. The discrepancy is either positive or negative, which means that both algorithms may deliver slightly suboptimal results. Figure 2(c) displays $\bar{P}_{min}$ obtained using the direct-rule-based method, as a function of $N$. At $d_{min}/a = 3$, increasing $N$ up to 100 may provide energy benefit, but the far-field approximation may fail at this small $d_{min}/a$. For more trustworthy values $d_{min}/a = 4$ and 5, the optimization results show that schools of $N \approx 15$ swimmers are practically as energy-saving as the largest schools.

4. Concluding remarks

Our model suggests that self-propelled 3-D swimmers can exploit each other's induced far field such that the entire school acquires a considerable net energy benefit. This gain is smaller than in two dimensions, but it is non-negligible (several per cent) if the swimmers are not far away from each other ($d_{min}/a \le 5$).
Far-field hydrodynamic interaction in a group of swimmers

The optimal schools save energy by exploiting the following properties: (i) $\partial P/\partial U$ of an isolated swimmer is positive (see Appendix A); (ii) the flow is slowed down everywhere upstream of the head of a swimmer, and everywhere downstream of the tail of a swimmer, except in the wake (see figure 1d); (iii) the total flow is constructed by superposition of individual swimmer flows in the regions exterior to the boundary layers and wakes. The first two properties make it possible to find an arrangement where every swimmer saves energy: arrange them in a chain so that they are all upstream (downstream) of each other's head (tail). The second and third properties make it so that increasing the number of swimmers will increase the energy savings. By adding another swimmer in a chain formation, all the other swimmers will feel an even slower flow, thereby increasing energy savings. However, the third property implies that the distances between the swimmers are large enough to ensure that the inflow velocities evaluated at their central points are representative of the average inflow velocity over the body.

The interactions based on the far-field interference are distinct from the vortex phase matching (Li et al. 2020) and the leading-edge vortex enhancement (Seo & Mittal 2022) mechanisms of two fishes that swim near each other. The latter require kinematic adjustment of the follower's tail beat to harvest the leader's wake with the maximum efficiency. The far-field interaction mechanism, on the other hand, depends on the time-averaged velocity. It is weak if we consider only a pair of swimmers ($N = 2$). But the average per capita energy benefit increases considerably with the number of individuals in the school, up to $N \approx 15$.

Optimal formations found by the optimizer in this study are chain formations. This formation exists in fish and squid schools, although it is not the most typical one. Some further enhancement of the model may be necessary, but it is also possible that energy optimization is not always the main objective.

Our approach is in many respects similar to well-known low-Reynolds-number models. But while the velocity superposition is exact in the limit of small Re, it is only a linearized approximate method in our situation. The predictive capacity of the model is limited further by the simplifications of using only one control point per individual, only the axial velocity component and the steady flow approximation. Therefore, the exact amounts of the energy benefit need to be verified further by e.g. Navier-Stokes simulations. But the model suggests clearly that individual swimmers within a school can occupy positions according to the perceived time-averaged velocity such that the school would expend less energy than the same number of isolated individuals swimming at the same speed over the same distance.

**Acknowledgements.** D.K. and J.S. thank the Bayreuth Humboldt Centre for the support of this project. Three anonymous reviewers are thanked for constructive comments and suggestions.

**Funding.** G.L. is funded by Japan Society for the Promotion of Science (JP20K14978).

**Declaration of interests.** The authors report no conflict of interest.

**Author ORCIDs.**

© Gen Li https://orcid.org/0000-0002-4423-3266;

© Jörn Sesterhenn https://orcid.org/0000-0003-1309-2306;

© Ramiro Godoy-Diana https://orcid.org/0000-0001-9561-2699;

© Benjamin Thiria https://orcid.org/0000-0002-2449-1065;

© Dmitry Kolomenskiy https://orcid.org/0000-0003-0107-6894.
<table><caption>Table 1. Parameters of baseline individual swimmer simulations.</caption><thead><tr><th>Parameters</th><th>Values</th></tr></thead><tbody><tr><td>Time step</td><td>0.0015625 s</td></tr><tr><td>Inner iterations</td><td>10</td></tr><tr><td>Computational domain</td><td>Length 12 m, Diameter 6 m</td></tr><tr><td>Total mesh grid number</td><td>Static 5.2 million, Rotational 0.9 million</td></tr><tr><td>Boundary layer thickness</td><td>Static 0.015 m, Rotational 0.005 m</td></tr><tr><td>Maximum grid size</td><td>Static 0.16 m, Rotational 0.01 m</td></tr></tbody></table>

<table><caption>Table 2. Computational results of individual swimmer performance.</caption><thead><tr><th>Case</th><th>Propeller rotation (rpm)</th><th>Cruising speed (m s<sup>-1</sup>)</th><th>Drag (N)</th><th>Net force (N)</th><th>Thrust (N)</th><th>Power (W)</th></tr></thead><tbody><tr><td>90 % speed</td><td>288</td><td>1.467</td><td>15.770</td><td>0.046</td><td>15.724</td><td>32.985</td></tr><tr><td>98 % speed</td><td>313.6</td><td>1.597</td><td>18.390</td><td>0.039</td><td>18.351</td><td>41.959</td></tr><tr><td>100 % speed</td><td>320</td><td>1.630</td><td>19.075</td><td>0.033</td><td>19.042</td><td>44.436</td></tr><tr><td>102 % speed</td><td>326.4</td><td>1.663</td><td>19.772</td><td>0.034</td><td>19.739</td><td>46.995</td></tr><tr><td>110 % speed</td><td>352</td><td>1.793</td><td>22.680</td><td>0.029</td><td>22.651</td><td>58.211</td></tr></tbody></table>

## Appendix A. The CFD simulation on an individual swimmer

The CFD simulation on an individual swimmer has been introduced earlier in Li *et al.* (2023). In this study, we have performed Reynolds-averaged Navier-Stokes simulations using commercial software Siemens Simcenter STAR-CCM+ 2020.1. The swimmer is an autonomous underwater vehicle with a propeller. Its body is axisymmetric, of 3.4 m length and 0.35 m diameter. The geometrical centre of the body is located 1.64 m downstream from the nose.

The CFD software and turbulence model were validated previously in a study focused on rotational wind blades (Fang *et al.* 2021). The propeller zone, which rotates, features the finest computational mesh, while the other parts are static regions covered by a Cartesian mesh. The static mesh is refined at the swimmer surface and rear zone. The boundary condition at the frontal surface of the computational domain is set as a velocity inlet at various values in each case. The rear surface of the computational domain is set as a zero-pressure outlet. The side surfaces are set as symmetry boundaries. The baseline parameters used in the simulations are listed in table 1.

A cruising propeller rotational speed 320 rpm was established based on realistic data. Through multiple iterative simulations, the equilibrium condition (swimmer's drag equal to thrust) at cruising speed $U_0 = 1.63 \text{ m s}^{-1}$ was determined. At this cruising speed, the propeller power and the flow field around the swimmer were recorded. Additionally, further simulations were conducted to compute the derivatives of power and flow field with respect to the cruising speed $U_0$. Using the same mesh, multiple simulations were implemented to approach the equilibrium condition at speeds of 90%, 98%, 102% and 110% of the cruising speed $U_0$. The results are presented in table 2.

To validate the mesh resolution and obtain detailed information about the wake field, an additional simulation with refined wake-zone mesh was implemented (parameters shown in table 3). Though the wake-zone mesh resolution in the additional simulation is several times higher than that of the baseline simulations, the extra refinement of the
<table><caption>Table 3. Parameters of swimmer simulations using a fine mesh.</caption><thead><tr><th>Parameters</th><th>Values</th></tr></thead><tbody><tr><td>Time step</td><td>0.0015625 s</td></tr><tr><td>Inner iterations</td><td>10</td></tr><tr><td>Computational domain</td><td>Length 24 m, Diameter 12 m</td></tr><tr><td>Total mesh grid number</td><td>Static 23.5 million, Rotational 0.9 million</td></tr><tr><td>Boundary layer mesh thickness</td><td>Static 0.015 m, Rotational 0.005 m</td></tr><tr><td>Maximum grid size</td><td>Static 0.16 m, Rotational 0.01 m, Wake zone 0.01 m</td></tr></tbody></table>

<table><caption>Table 4. Comparison between computational results using two different meshes.</caption><thead><tr><th>Case</th><th>Fixed propeller<br>rotation (rpm)</th><th>Fixed cruising<br>speed (m s<sup>-1</sup>)</th><th>Drag (N)</th><th>Net force<br>(N)</th><th>Thrust<br>(N)</th><th>Power<br>(W)</th></tr></thead><tbody><tr><td>Baseline mesh</td><td>320</td><td>1.630</td><td>19.075</td><td>-0.033</td><td>19.042</td><td>44.436</td></tr><tr><td>Fine mesh</td><td>320</td><td>1.630</td><td>19.021</td><td>0.021</td><td>19.042</td><td>44.440</td></tr></tbody></table>

<figure><img src="image_8.png" /><figcaption>Figure 6. Power expenditure as a function of swimming speed. The speed is normalized by $U_0$, and the power is normalized by $P_0 = P|_{U_0}$. Asterisk markers correspond to the values obtained from the CFD simulations. The blue dashed line is a power-law fit $P/P_0 = (U/U_0)^{2.7947}$. The red solid line is the linear fit used in the text. Over the interval $U/U_0 \in [0.9, 1.1]$, the discrepancy is less than 3 %.</figcaption></figure>

mesh does not bring any significant influence in the swimmer propulsive performance at the cruising speed (table 4), suggesting that the baseline mesh resolution was already sufficient to predict the propulsive performance. The detailed wake information obtained by this additional simulation was input into our theoretical model to describe the wake structure.

The power expenditure curve is obtained from CFD simulations of an isolated swimmer. We performed a series of simulations at different speeds. The results are shown in figure 6. The power $P$ increases nonlinearly with the swimming speed $U$. However, in the neighbourhood of $U_0$, this nonlinear dependence can be approximated by a linear dependence, following the usual procedure of linearization.
<figure><img src="image_9.png" /><figcaption>Figure 7. Velocity excess/deficit in the wake. Dashed lines show the CFD results, and solid lines show the far-field approximation.</figcaption></figure>

## Appendix B. Wake velocity profile fitting

The velocity profile in Korobko et al. (1986) consists of one jet region (excess velocity in the inner part of the profile) and one wake region (velocity deficit in the outer part due to the boundary layer separated from the swimmer's body). However, in our CFD simulation, the jet is produced by a propeller that has a central body, and the central body produces drag (velocity deficit). Then the propeller blades produce an annular jet (excess velocity). Finally, in the outer region, there is still a velocity deficit region due to the separated boundary layer. In order to have these three regions in the velocity profile, we introduce a simple (perhaps the simplest possible) modification of Korobko's formula. Note that while Korobko's formula is an exact result for a laminar momentumless wake, our intent is to have a semi-empirical formula suitable for the turbulent wake in our particular situation. Figure 7 displays a comparison between the fit and raw CFD data. The staircase in the CFD data is due to the low order of interpolation on the non-uniform grid.

Figure 8 displays the CFD output data and the algebraic fit of the velocity sampled along the axis of the wake. The vertical line is drawn at the axial position $x_{flex}$ where the line has an inflection point. This location is used for determining the fitting coefficients. In the central part of the CFD curve, the fit is very close to the CFD data, therefore choosing any $x$ in that region for measuring velocities does not affect the results of fitting.

## Appendix C. Vortex ring fitting

Figure 9 shows the results of a sensitivity analysis, when the velocity profiles are fitted on a line at different distances from the swimmer axis. The results in the main text are obtained using the distance $r_s^* = 4.36a$, while figure 9 shows $\gamma_j$ and $\Gamma$ calculated using several smaller or larger values of the distance. When the distance $r_s$ is small, the far-field approximation fails, and the fitting is inaccurate. This happens at $r_s = 0.6r_s^*$. When the distance $r_s$ is large, the distance between the $j$th and $(j+1)$th vortex rings is much smaller than $r_s$. In this situation, the entire system of vortex rings induces the velocity closely
Far-field hydrodynamic interaction in a group of swimmers

<figure><img src="image_10.png" /><figcaption>Figure 8. CFD output data and algebraic fit of the velocity on the axis of the wake. The vertical line corresponds to the axial position $x_{flex}$ where the line has an inflection point.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Figure 9. Sensitivity of (a) the vortex ring circulation coefficients $\gamma_j$, and (b) total circulation $\Gamma$, on the distance $r_s$ from the axis to the fitting line.</figcaption></figure>

matching a single vortex ring with the circulation $\Gamma$. For that reason, the linear regression matrix becomes ill-conditioned; small numerical errors in the velocity sampling amplify, and lead to spurious dispersive error in the distribution of $\gamma_j$. In the intermediate range $r_s \in [0.8, 1.2]$, the above-mentioned sources of error vanish, the results are not sensitive to $r_s$, and this is the correct fit.

The results in figure 9 were obtained using $n_r = 21$ vortex rings. Smaller values of $n_r$ deteriorate the approximation accuracy at distances of the same order of magnitude as the vortex ring spacing or less. Increasing $n_r$ at a constant $r_s$ provokes the spurious dispersive error as discussed above. We found that the value $n_r = 21$ is well adapted for the parameters of our CFD data set.

REFERENCES

ASHRAF, I., GODOY-DIANA, R., HALLOY, J., COLLIGNON, B. & THIRIA, B. 2016 Synchronization and collective swimming patterns in fish (*Hemigrammus bleheri*). *J. R. Soc. Interface* **13** (123), 20160734.

BATCHELOR, G.K. 2000 *An Introduction to Fluid Dynamics*, Cambridge Mathematical Library edn, pp. 521–522. Cambridge University Press.
BRUMLEY, D.R., WAN, K.Y., POLIN, M. & GOLDSTEIN, R.E. 2014 Flagellar synchronization through direct hydrodynamic interactions. *eLife* **3**, e02750.

FANG, Y., LI, G., DUAN, L., HAN, Z. & ZHAO, Y. 2021 Effect of surge motion on rotor aerodynamics and wake characteristics of a floating horizontal-axis wind turbine. *Energy* **218**, 119519.

FILELLA, A., NADAL, F., SIRE, C., KANSO, E. & ELOY, C. 2018 Model of collective fish behavior with hydrodynamic interactions. *Phys. Rev. Lett.* **120**, 198101.

GAZZOLA, M., TCHIEU, A.A., ALEXEEV, D., DE BRAUER, A. & KOUMOUTSAKOS, P. 2016 Learning to school in the presence of hydrodynamic interactions. *J. Fluid Mech.* **789**, 726–749.

HANSEN, N., MÜLLER, S.D. & KOUMOUTSAKOS, P. 2003 Reducing the time complexity of the derandomized evolution strategy with covariance matrix adaptation (CMA-ES). *Evol. Comput.* **11** (1), 1–18.

HERSKIN, J. & STEFFENSEN, J.F. 1998 Energy savings in sea bass swimming in a school: measurements of tail beat frequency and oxygen consumption at different swimming speeds. *J. Fish Biol.* **53** (2), 366–376.

KANSO, E. & TSANG, A.C.H. 2014 Dipole models of self-propelled bodies. *Fluid Dyn. Res.* **46** (6), 061407.

KOROBKO, V.I., SHASHMIN, V.K. & SHUL’MAN, Z.P. 1986 Contribution to the theory of linear momentumless wake. *Izv. Akad. Nauk SSSR* **2**, 35–39.

LI, G., DUAN, L., GODOY-DIANA, R. & THIRIA, B. 2023 Group formation of autonomous underwater vehicles that optimizes energetic efficiency in cruising. In *IEEE Underwater Technology Tokyo, Japan*, pp. 1–6. IEEE.

LI, G., KOLOMENSKIY, D., LIU, H., THIRIA, B. & GODOY-DIANA, R. 2019 On the energetics and stability of a minimal fish school. *PLoS ONE* **14** (8), e0215265.

LI, L., NAGY, M., GRAVING, J.M., BAK-COLEMAN, J., XIE, G. & COUZIN, I.D. 2020 Vortex phase matching as a strategy for schooling in robots and in fish. *Nat. Commun.* **11**, 5408.

LI, L., RAVI, S., XIE, G. & COUZIN, I.D. 2021 Using a robotic platform to study the influence of relative tailbeat phase on the energetic costs of side-by-side swimming in fish. *Proc. R. Soc. A* **477** (2249), 20200810.

SEO, J.-H. & MITTAL, R. 2022 Improved swimming performance in schooling fish via leading-edge vortex enhancement. *Bioinspir. Biomim.* **17** (6), 066020.

SIRVIENTE, A.I. & PATEL, V.C. 2000 Wake of a self-propelled body. Part 1. Momentumless wake. *AIAA J.* **38**, 613–619.

WEIHS, D. 1975 Some hydrodynamical aspects of fish schooling. In *Swimming and Flying in Nature: Volume 2* (ed. T.Y.-T. Wu, C.J. Brokaw & C. Brennen), pp. 703–718. Springer.

WEIMERSKIRCH, H., MARTIN, J., CLERQUIN, Y., ALEXANDRE, P. & JIRASKOVA, S. 2001 Energy saving in flight formation. *Nature* **413** (6857), 697–698.

YUAN, Z.-M., CHEN, M., JIA, L., JI, C. & INCECIK, A. 2021 Wave-riding and wave-passing by ducklings in formation swimming. *J. Fluid Mech.* **928**, R2.

ZHANG, P., PETERSON, S.D. & PORFIRI, M. 2023 Dipole- and vortex sheet-based models of fish swimming. *J. Theor. Biol.* **556**, 111313.