Far-field hydrodynamic interaction in a group of swimmers

Abstract 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.


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;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.

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 Γ .
The Stokes stream function of the ring can be expressed in terms of complete elliptic integrals K and E (Batchelor 2000), (2.1) is the elliptic modulus, and The axial velocity u = r −1 ∂ψ/∂r can also be reduced to the elliptic integrals using the identities When r is small, numerical evaluation of (2.4) becomes problematic.The formula for u ring at the axis is derived by keeping only the leading-order terms E(k 2 ) ≈ π/2 and K(k 2 ) ≈ π/2 in the respective series expansions.The axial velocity component at r = 0 becomes equal to Γ Υ ring (x, 0), with (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 y/a -6 -4 -2 0 x/a x/a momentumless wake velocity decays exponentially with r; see § 2.2).Let Γ 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): where x j = x nose + ( j − 1/2)l/n r is the position of the jth vortex ring, and x nose corresponds to the front point of the body.The circulation distribution coefficients γ j are calculated as follows.First, dimensional circulations Γ 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 Γ is found by summing all vortex ring circulations Γ j .Finally, the coefficients are evaluated as γ j = Γ j /Γ .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 Γ .

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 where it is postulated that the radial profile is a function of the similarity variable η = r a ξ p Re p , (2.9) with ξ = x/a, and the Reynolds number Re = U ∞ a/ν being based on U ∞ , the far-field velocity that is equal in magnitude to the swimming speed, and ν, 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 ξ l/a.The radial profile of the velocity, g(η), depends on the spatial distribution of the momentum excess and deficit in the wake.We are interested particularly in the case of g(η) changing sign twice (see Appendix B for more detail): where η 1 , η 2 and η 3 are constant parameters.Since the wake is momentumless and axisymmetric, it is required that It follows that only two of the three parameters can be adjusted freely, since (2.12) The far-field approximation to the velocity due to both the bound vorticity and the momentumless wake is constructed as (2.13) The postulate that the momentumless wake intensity scales linearly with the bound circulation Γ 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 Γ .

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: 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, where Υ (x i ) is calculated as defined in (2.13).
The circulation Γ 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 : where Γ , 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 (2.18) When positions x 1 , . . ., x N are fixed a priori, (2.18) constitutes a system of N linear equations with N unknowns Γ 1 , . . ., Γ N .By solving this system and substituting the solution into (2.15),we evaluate the induced velocities u ij for all i = 1, . . ., N and j = 1, . . ., 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 (2.19)

Optimization objective
We aim to minimize the normalized school-average mechanical power expenditure, 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 , . . ., 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 where σ p is a parameter that controls the minimum allowed distance between swimmers as a multiple of a, and A p 1 is a constant penalty parameter (we use A p = 10 4 ).The pairwise lateral separation distance between the swimmer centroids (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 (2.23)

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ρU 3 0 πa 2 , where ρ is the water density.
The total circulation Γ 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 = Γ γ j were treated as unknown coefficients, thus their values were found using the linear regression.The total circulation Γ | U 0 = −16.5U0 a was found by the summation of G j over j from 1 to n r .Finally, the distribution coefficients γ j = G j /Γ 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 /ν = 2.8 × 10 5 .The radial profile g(η) was adjusted to the data at distance 3.26a downstream from the tail.This resulted in the parameter values η 1 = 2.15, η 2 = 5.5, η 3 = 3.34, and a profile shown in figure 1(c).The parameter values C = 0.516 and p = −0.176ensure 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.

Parameters of optimal schools
The optimization problem consists in finding the swimmer positions x 2 , . . ., 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 σ p , which enters in the penalty function (2.21).All optimal formations that we found activated the constraint d min = σ p a. Figure 2(a) shows how the school average power Pmin 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 σ 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 Pmin suggests that reducing the lateral distance in a school can enhance energy-saving effects.Due to this, d min /a would tend to 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 Pmin 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 − Pmin 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 any swimmer has at most two neighbours in its vicinity.Therefore, the induced velocities remain small.

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 ρ and the polar angle θ 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, ρ/a = 10.8, θ = 34 • 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 . 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 2 x y z is obtained from S 1 xyz by translation of the origin from the first swimmer to the second swimmer.The planes S 2 x y and S 1 xy coincide, because the z-coordinate of the second swimmer is zero.The angle θ lead is between the negative semi-axis S 1 x and the line that connects S 1 with S 2 .Point O marks the perpendicular projection of the third swimmer on the S 2 x axis.Axis O y is parallel to S 2 y , and it resides on the S 2 x y plane.Angle θ hind is the angle between the negative semi-axis S 2 x and the line that connects S 2 with S 3 .Angle φ 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.
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): ρ lead is slightly less and θ 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 ρ hind and θ hind coincide with ρ lead and θ lead , respectively.The φ 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.

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 ρ lead , θ lead , ρ hind and θ hind are on the rear updraft line, as in figures 5(a,b) at the required 974 A34-11 d min , and φ hind = 120 • .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 φ angle is −120 • .Thus we ensure that the fourth swimmer is in the updraft of all those in front.The fifth swimmer is added with φ = 120 • , 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 Pmin is smaller for larger schools.To estimate the lower bound of Pmin , 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 Pmin are indistinguishable from the (3N − 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 Pmin 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 ≈ 15 swimmers are practically as energy-saving as the largest schools.

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 ≤ 5).
The optimal schools save energy by exploiting the following properties: (i) ∂P/∂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 ≈ 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.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.

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 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 matching a single vortex ring with the circulation Γ .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 γ j .In the intermediate range r s ∈ [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.

Figure 1 .
Figure 1.(a) Schematic representation of the vortex ring and momentumless wake model.Dimensional circulations Γ j are related to the dimensionless coefficients γ j by Γ j = γ j Γ .(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.

Figure 2 .
Figure 2. (a) School average power Pmin 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 − Pmin 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 − Pmin = 0.220(d min /a) −1.97 for N = 2, 1 − Pmin = 0.286(d min /a) −1.90 for N = 3, 1 − Pmin = 0.354(d min /a) −1.82 for N = 7, and 1 − Pmin = 0.392(d min /a) −1.83 for N = 11.(c) Plots of Pmin obtained using the sequential increment of N, as functions of N, for d min /a = 3, 4 and 5.

Figure 3 .
Figure 3. Optimized schools of N = 11 individuals with d min /a = 5: (a) front 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.

Figure 5 .
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).

Figure 8 .Figure 9 .
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.

Table 1 .
Parameters of baseline individual swimmer simulations.

Table 2 .
Computational results of individual swimmer performance.