# Intermittent versus continuous swimming: An optimization tale

Gen Li, Dmitry Kolomenskiy, Hao Liu, Ramiro Godoy-Diana, Benjamin Thiria

► To cite this version:

Gen Li, Dmitry Kolomenskiy, Hao Liu, Ramiro Godoy-Diana, Benjamin Thiria. Intermittent ver- sus continuous swimming: An optimization tale. Physical Review Fluids, 2023, 8 (1), pp.013101. 10.1103/PhysRevFluids.8.013101 . hal-04290424

HAL Id: hal-04290424

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

Submitted on 16 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
# Intermittent vs. continuous swimming:
an optimization tale

Gen Li*

*Japan Agency for Marine-Earth Science and Technology (JAMSTEC), 3173-25, Showa-machi, Kanazawa-ku, Yokohama-city, Kanagawa, 236-0001, Japan*

Dmitry Kolomenskiy

*Center for Design, Manufacturing and Materials (CDMM), Skolkovo Institute of Science and Technology, Bolshoy Boulevard 30, bld. 1 Moscow, 121205, Russia*

Hao Liu

*Graduate School of Engineering, Chiba University, 1-33, Yayoicho, Inage-ku, Chiba-shi, Chiba, 263-8522, Japan*

Ramiro Godoy-Diana

Benjamin Thiria

*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*

## Abstract

Intermittent swimming, also termed as “burst-and-coast swimming”, has been reported as a strategy for fish to enhance their energetical efficiency. Intermittent swimming involves additional control parameters, which complexifies its understanding by means of quantitative and parametrical analysis, in comparison with continuous swimming. In this study, we used a hybrid computational fluid dynamic (CFD) model to assess the swimming performance in intermittent swimming parametrically and quantitatively. A Navier-Stokes solver is applied to construct a database in the multi-dimensional space of the control parameters to connect the undulation kinematics to swimming performance. Based on the database, an indirect numerical approach named “gait assembly” is used to generate arbitrary burst-and-coast gaits to explore the parameter space. Our simulations directly measured the hydrodynamics and energetics under the unsteady added-mass effect during burst-and-coast swimming. The results suggest that the instantaneous power of burst is basically determined by undulatory kinematics. The results show that the energetical performance of burst-and-coast swimming can be better than that of continuous swimming, but also that an unoptimized burst-and-coast gait may become very energetically expensive. These results shed light on the mechanisms at play in intermittent swimming, enabling us to better understand fish behavior and to propose design guidelines for fish-like robots.

<footnote>* ligen@jamstec.go.jp</footnote>
# I. INTRODUCTION

Undulatory propulsion is commonly adopted by fish in their locomotion [1], [2]. A great number of biomechanical and physiological studies on fish locomotion base on a cyclic swimming state, where continuous undulation enables fish to maintain steady cruising velocity and energy expenditure [3]. On the other hand, like many other animals, fish may also perform intermittent locomotion [4], [5]. The intermittent locomotion may be applied in predator-prey interaction, sensing, habitat assessment, and cruising motion [5]–[8]. Intermittent locomotion in fish swimming is also termed burst-and-coast swimming, which is a two-phase periodic behavior consisting of an active propulsive phase followed by a passive gliding phase with the body straightened [4], [9]. As shown in Fig. 1(a), during the burst phase, the fish undulates its body and caudal fin to gain forward momentum, while during the coast phase, the fish keeps its body straight and consumes its forward momentum for traveling distance. By repeating the burst-and-coast process, the average swimming velocity can be sustained at a desired level.

Since the pioneering studies by Weihs et al. [9], [10], intermittent swimming has been investigated from the bio-mechanical perspective. Intermittent swimming is considered as a means to improve the energetic performance in linear swimming (e.g. [9]–[11]). A theoretical basis is provided by the ‘Bone–Lighthill boundary-layer thinning hypothesis’ [3]. According to it, the lateral motion of a fish body element may reduce the thickness of the boundary layer, and thus increase the friction drag. Propulsive movement of the body and appendages are expected to increase the drag by a factor of 2-5 [3], [12]–[14]. Thus, by ceasing its body undulation a fish may exploit the lower drag of a rigid body and eventually reduce the energy consumed to overcome drag [15]. Such advantages are confirmed by experiment and simulation [12], [16]–[18]. It has been hypothesized that the advantage of an intermittent swimming strategy mainly benefits BCF (body and caudal fin) swimmers using sub-carangiform and carangiform swimming, as MPF (median-paired fin) swimmers and, to a lesser extent, thick-bodied thunniform swimmers already maintain a rigid body, thus little advantage would be gained by adopting burst and coast swimming [19]. An analytical study by Blake [20] reports that the frequent burst-and-coast swimmers are characterized by a fineness ratio (the ratio of the length of a body to its maximum width) around 5, and that fish with lower fineness ratios use less burst-and-coast swimming. Nevertheless, in a computational fluid dynamic study by Xia et al. [21], they report that an energetic improvement in the burst-and-coast swimming of virtual thunniform swimmers is possible (our remarks on this study are discussed in § IV. A). Burst-and-coast swimming gait and performance are also correlated with the cruising velocity. Our most recent experimental-numerical study shows that fish are able to optimize their gait for minimal cost of transport, by modulating a unique intrinsic cycle to sustain the demanded velocity [22]. Interestingly, a few recent studies report that in some circumstances intermittent swimming may cost more energy than cyclic swimming (e.g. [23], [24]), which suggests the energy-saving function of burst-and-coast swimming is not inevitable and may require optimization in control parameters. Overall, since the hydrodynamic details of burst-and-coast swimming are usually hard to access by experimental measurement, and the intermittent behavior is associated to a complex control parameter space, there is a strong lack of quantitative and parametrical analyses concerning
intermittent swimming.

The purpose of this paper is to parametrically assess the swimming performance in intermittent swimming and to compare the performance of intermittent and continuous swimming patterns. To accurately examine the highly dynamic instantaneous performance during burst-and-coast swimming, we use a direct computational fluid dynamic (CFD) model of a self-propelled fish based on the Navier-Stokes equations [14], [25], [26]. To parametrically explore the gait and its corresponding swimming performance in intermittent swimming, we use an indirect numerical approach named as *gait assembly* [22], which generates arbitrary intermittent gaits based on the database constructed by multiple direct CFD simulations. The rest of this paper is organized as follows: In §II we investigate the hydrodynamic characteristics of burst and coast swimming using a Navier-Stokes solver; In §III, we optimize the burst-and-coast gait based on the database constructed in §II and explore the parameter space. We investigate the optimal gait in burst-and-coast swimming and make a comparison with a continuous swimming gait. Note that the numerical methods used in §II and §III are respectively explained in the first sub-section of each section. Comparison between our results and previous studies, as well as the inspirations for future research are discussed in §IV.
## II. HYDRODYNAMIC MODELLING OF BURST-AND-COAST SWIMMING

### A. Numerical approach in this section

We use a previously validated three-dimensional computational fluid dynamic (CFD) approach based on an overset-grid finite-volume method to simulate a self-propelled model fish [25]–[27] (for more information including the numerical validation, see Electronic Supplementary Materials (ESM)). The model fish swims freely in the horizontal plane (3 degrees of freedom (DoF), including translational motion in the horizontal plane and rotation about the vertical axis). The approach comprises surface models of a typical carangiform fish shape (Fig. 2(a), model fish body length: 2cm, dimension in CFD mesh: 121 × 97), and local fine-scale body-fitted grids plus a large stationary global grid (Fig. 2(c)) to calculate the flow patterns around the fish with sufficient resolution. As shown in Fig. 2(b), the instantaneous body shape is driven by sinusoidal variation of the midline, cf. [26],

$$H(l, t) = \alpha \cdot l^2 \cdot \sin\left(\frac{2\pi l}{\lambda} - 2\pi f t\right) \qquad (\text{Eq.1})$$

where $l$ is the dimensionless distance from the snout along the longitudinal axis of the fish based on the length of the fish model $L$; $H(l, t)$ is the dimensionless lateral excursion in a frame attached to fish head at time $t$; $\alpha$ is the dimensionless amplitude control factor; $\lambda$ is the length of the body wave. Since we simulate a carangiform swimmer [1] and the wavelength of typical carangiform swimmer is approximately equivalent or greater than $1L$ [28], we define $\lambda = 1.1L$; $f$ is the tail beat frequency. This equation may cause total body length along the midline to vary during the tail beat, which is corrected by a procedure that preserves the lateral excursion $H(l, t)$ and ensures constant body length (see § B-3, ESM).

Three-dimensional Navier-Stokes (NS) equations in an inertial frame of reference are solved. The solving process is implemented using the finite volume method (FVM), based on a multi-block, an overset mesh system and an inter-block communication algorithm. The governing equations for the fluid solution are the three-dimensional, incompressible and unsteady NS equations written in strongly conservative form for mass and momentum [27]. To accelerate the computation and improve the robustness during iteration, the artificial compressibility method is adopted by adding a pseudo time derivative of pressure to the continuity equation [27]. For further details of the NS solver, see § B, ESM.

While the deformation of the central body axis of the fish is prescribed, the center-of-mass (CoM) movements and body orientation are determined by the hydrodynamic and inertial forces generated by the swimming model fish. The forces acting on the body and its motion are obtained through coupling the hydrodynamic and body dynamic solutions, which ensures that the motions of the fish correspond to the hydrodynamic and inertial forces exerted on the fish. The range of Reynolds numbers in this study is below 6000, turbulence models are not used, and the grid resolution at Re=6000 has been justified in our previous study [29].

A swimming fish generates unsteady pressure and shear stress at all locations on its body. The shear component of the stress predominantly converts into skin friction drag, while the
pressure component of the stress, particularly in the posterior fish body, changes its direction rapidly during fish undulation. Hence, in this study, the instantaneous thrust (drag) at each time-step is defined as the sum of the forward (backward) components of pressure and shear stress over all fish surface elements (see Fig. S4, ESM). Each surface element could contribute to thrust or drag at different time steps even within one tail beat cycle, as well as to thrust and drag simultaneously (e.g. generating pressure-based thrust and shear-based drag). Such definition may effectively separate thrust and drag forces during unsteady, undulatory swimming.

In this paper, power refers to “mechanical power”, defined as the sum of the hydrodynamic and body inertial powers:

$$P = P_{\text{hydro}}+P_{\text{body}} \qquad (\text{Eq.2})$$

Hydrodynamic power is calculated as the sum of the hydrodynamic work on the body surface, such that:

$$P_{\text{hydro}} = \oint \oint^{\text{surface}} (\mathbf{f} \cdot \mathbf{U}) dS \qquad (\text{Eq.3})$$

where $P_{hydro}$ is the hydrodynamic power; $ds$ denotes surface element, $\mathbf{f}$ is the hydrodynamic stress vector acting on the surface element; $\mathbf{U}$ is the velocity vector on this surface element.

Body inertial power is computed as the sum of the kinetic energy change rate of all body elements (inside the body), such that:

$$P_{\text{body}} = \oint \oint^{\text{body}} (\rho \cdot \mathbf{a} \cdot \mathbf{U}) dV \qquad (\text{Eq.4})$$

where $P_{body}$ is the body inertial power; $dV$ denotes body volume element, $\mathbf{a}$ and $\mathbf{U}$ are the acceleration and velocity vectors of each body volume element, respectively, and $\rho$ is the local density.

Note that based on Eq. 2, during the coast phase, the mechanical power is zero.

The dimensionless cost of transport (CoT) is defined as

$$\Omega^* = \frac{E}{m g S} \qquad (\text{Eq.5})$$

where $E$ denotes energy consumption, $m$ denotes body mass, $g$ denotes gravitational acceleration and $S$ denotes travelled distance.

## B. Hydrodynamics and energetics of a non-repeated burst-and-coast bout

Using the Navier-Stokes solver introduced in §II.A, we simulated a non-repeated swimming bout consisting of developed burst and coast phases (Fig. 3; Supplementary Video 1). As shown in Fig. 3 and Fig. 4(a), the fish started swimming from rest (body-axis #0 in Fig. 4(a)) with a body deformation driven by Eq. 1 at $f = 10Hz$, $\alpha = 0.11$, based on experimental observations [28], [30]. The cyclic swimming continued for 15 tail-beat cycles (until body-axis #15 in Fig. 4(a)) at the end of which the fish has accelerated to a nominally terminal velocity, then stopped its undulation (body-axis #15.5 in Fig. 4(a)) and the fish body was kept straight. The fish coasted for another 1.5s, slowing down until almost reaching a full stop (body-axis #30 in Fig. 4(a)). The flow field produced by this burst-and-coast bout is shown by Fig. 3 (Video: SV1): the fish
generated jet flows and vortex rings behind itself. Corresponding to the sharp acceleration at the beginning, those initial jets merged into a strong backward jet (① in Fig. 3). As the acceleration rate decreased, the lateral velocity components gradually dominated the jet flows (② in Fig. 3). On the fish body, the anterior part was rigid and surrounded by relatively stable boundary vorticity—a shear layer surrounding the fish body. At the posterior body, the surrounding boundary vorticity was disturbed due to the undulation, and vortices shed in a staggered pattern (③ in Fig. 3) similar to the observation by Wu et al. [12]. As the fish swims further, those starting jets merged into a large backward jet zone (④ in Fig. 3). Unlike the early vortices that merged into two main vortices, those vortices generated later were in a staggered pattern (⑤ and ⑥ in Fig. 3). The angle between the vortex ring axis and the backward direction (opposite to the swimming direction) was observed to change—early vortex rings generated during sharp acceleration have smaller angles, while vortex rings generated during stable swimming contain a larger portion of lateral velocity component, which agreed with the observation by Akanyeti et al. [31]. When the fish stopped undulating and began to coast, its tail stopped shedding staggered-pattern vortices, and the entire body produced an elongated boundary vorticity layer (⑦ in Fig. 3) that covered the gliding trajectory of the fish (⑧ in Fig. 3).

The undulation formed periodic fluctuations in velocity, power, thrust and drag time-sequences (gray curves in Fig. 4(b, d, f, h)). These fluctuations were low-pass-filtered (function *lowpass*, MATLAB R2020b, cut-off frequency = *f*) to clearly demonstrate the average trend of those performance parameters (black curves in Fig. 4(b, d, f, h)). As shown in Fig. 4(b), in the burst phase, the acceleration was rapid at the beginning, with a sharp increase of fish swimming velocity, but it attenuated to zero as the drag matches the thrust. In the coast phase, as the fish stopped undulating, its velocity dropped sharply at the beginning of the coasting phase and then the deceleration rate weakened as the velocity decreases. The shape of the velocity curve (Fig. 4(b)), consisting of a gentle-sloped top of the burst phase and a gentle-sloped bottom of the coast phase, agreed with the theoretical investigation by Videler and Weihs [10]. The trend of this velocity curve was determined by the characteristics of the thrust and drag of the swimming fish. The thrust of fish, as shown by Fig. 4(f), was most strong at the beginning of the burst phase, then gradually decreased as swimming velocity raised. During the coast phase, the thrust immediately dropped to zero as the undulation ceased. The drag on the fish showed positive correlation with the swimming velocity (Fig. 4(h)), however, at one specific velocity, the magnitude of drag differed dramatically between burst and coast phases, since the undulatory propulsion amplified the magnitude of drag [9].

In Fig. 4(i), a drag coefficient curve is drawn to demonstrate the dynamic change in drag during burst-and-coast swim (black solid curve, which is further decomposed into two parts respectively formed by the burst and coast phases). As a reference, a velocity-specified drag coefficient curve of a fish steadily gliding with a straight body is added to the plot (gray dashed line). The ratio between the instantaneous drag coefficients of burst-and-coast swim and steady glide, defined as drag amplifying factor $β$, was not a constant value as applied in previous analytical research (e.g. [10]), but varied in a wide range. At the beginning of the burst phase, $β ≈ 10$, representing a strong amplification in drag due to added mass effect; at the end of the burst phase, $β$ gradually diminished to approximately 2, representing that the added mass effect
is attenuated as the magnitude of acceleration diminishes, while a certain extent of drag amplification was still maintained due to the body undulation; as the fish transited to coasting, $\beta$ quickly became less than 1, which means that during the coast phase the added mass effect was reversed, resulting in a lower drag than that observed in a steady glide with the same velocity.

As shown in Fig. 4(g), during the burst phase, a stabilized velocity condition was achieved as the decreasing thrust and the increasing drag balanced at a specific velocity. In the coast phase, thrust dropped immediately to zero due to the absence of undulation, meanwhile, the drag also dropped sharply but not to zero because the fish moving body had to overcome the hydrodynamic drag and decelerated.

As to the power (Fig. 4(e)), it is worth noting that during the burst the relation between power and forward-swimming velocity differed much from that in steady swimming: the power is not proportional to velocity, instead, it raises to a maximum level immediately when the fish starts the burst, and is slightly reduced as the velocity increases, reaching a relatively stable level at a specific velocity. The power during the burst seems not strongly influenced by the instantaneous velocity. Apparently, the lateral undulating velocity mattered most for the power expenditure during the early stage of burst, rather than the forward component of the CoM velocity. In the coast phase, the power dropped immediately to zero.

The simulation results show that, during the burst, the fish had to spend extra effort to overcome the added-mass effect and lateral-side power expenditure. Here, we define the accumulated cost of transport to represent the cost of transport from the start to present time $t_p$:

$$ \Omega_{\text{accumulated}}^{*}(t_p) = \frac{E(t_p)}{mgS(t_p)} = \frac{\int_{0}^{t_p} P(t)dt}{mg \int_{0}^{t_p} U(t)dt} $$

where $E(t_p)$ denotes the accumulated energy consumption from the beginning to present time $t_p$, $S(t_p)$ denotes the accumulated travelled distance from the beginning to present time $t_p$. As shown in Fig.4(c), $\Omega_{\text{accumulated}}^*$ started extremely high at the beginning of the burst phase, then decreased as the velocity approaches the equilibrium state. This shows that, in this burst-and-coast process, fish invest extra energy during the burst-phase, and earn the benefit during the coast-phase.

### C. Instantaneous velocity and power during the burst phase

In this section, our analysis focuses on the swimming performance during the burst phase. By using the Navier-Stokes solver introduced in §II.A, we directly simulated 40 burst swimming processes from rest to full velocity under various undulatory kinematics: the combinations of tail beat frequency (2, 6, 10, 14 and 18 Hz) and amplitude ($\alpha$ =0.02, 0.06, 0.08, 0.10, 0.12, 0.14, 0.18 and 0.22 in Eq. 1). In all burst swimming processes, as the fish sufficiently accelerates itself and approaches cyclic regime, its thrust, velocity and power are termed as terminal thrust, velocity and power, respectively.

We obtained the time profiles of the forward swimming velocity and hydrodynamic power, which are shown, respectively, in Fig. 5(a) & (b) in the normalized form $U/U_{\text{terminal}}$ and $P/P_{\text{terminal}}$
as functions of the dimensionless time $t^* = t/(A^{*-1/2}f^{-1})$, where $A^* = A/L$. Normalization with the values $U_{terminal}$ and $P_{terminal}$ that correspond the end-point values in each transient reveals the shapes of typical time profiles of $U$ and $P$. In Fig. 5(a), considering the velocity curves passing from the origin and asymptotically approaching $U^* = U/U_{terminal} = 1$, we attempted to use an exponential function $U^*(t^*) = 1 - e^{a t^*} (a < 0)$ to fit the average curve (green curve) of all those burst curves. We obtained an exponential fit line that accurately represents the acceleration process in burst swim:

$$U^*(t^*) = 1 - e^{-6.49 t^*} \qquad (\text{Eq.6})$$

In Fig. 5(b), power curves also form a common pattern: the power raises sharply at the very initial stage of burst, and after that, though the fish is still accelerating, the power curves are predominantly constant. The results suggest that when the kinematics are specified, the burst-swimming velocity and power are basically determined by the burst undulation kinematics (in this study, tail-beat frequency and amplitude). We emphasize that the instantaneous power during the burst is not determined by the instantaneous swimming velocity. Our results suggest that the velocity during a burst can be approximated by an exponential fit line $U^*(t^*) = 1 - e^{a t^*} (a < 0)$, while the power during the entire burst can be approximated by a constant value equal to the terminal power.

The magnitude of $U_{terminal}$ and $P_{terminal}$ vary across the cases in a wide range, but Fig. 5(c) & (d) suggest that both can be approximated as simple functions of tail-beat frequency and amplitude. In Fig. 5(c), $Re_{terminal} \propto Sw^k$, where $Re_{terminal} = U_{terminal}L/v$ and $Sw = fAL/v$ (v: kinematic viscosity) is swimming number, defined as $k$ ranges dynamically between 1.1 and 1.6, for which Gazzola et al. [32] derive $k = 4/3$ at low $Re$ and $k = 1$ at high $Re$. Fig.5(d) demonstrates that $P_{terminal}$ was, in the majority of cases, reasonably well characterized by a constant value of the power coefficient $C_p = 0.65$, where $C_p = P_{terminal}/\rho L^3 f^3 A^2$. In Fig.5(d), wave-speed based Reynolds number was defined as $Re_{wave} = fL^2/v$. Large values of $C_p$ in Fig. 5(d) correspond to those cases where $Re_{wave}$ is small or $Re_{terminal}$ (the inset in Fig. 5(d)). The maps of $Re_{terminal}$, $T_{terminal}$ and $P_{terminal}$ in the $Re_{wave}-A^*$ plane are shown in panels (e), (f) and (g), respectively.
### III. OPTIMIZATION IN BURST-AND-COAST SWIMMING

#### A. Algorithm of burst-and-coast gait assembly and optimization

Fish usually perform an intermittent swimming gait formed by repeated burst-and-coast bouts, which determines a multi-dimensional control parameter space. The control parameters can be limited to four: the tail-beat frequency $f$ during burst swim, the tail-beat amplitude $A$ during burst swim, an upper velocity boundary $U_U$ (transitional point from coast to burst, Fig. 1(b)) and a lower velocity boundary $U_L$ (transitional point from coast to burst, Fig. 1(b)). Fully covering this four-dimensional parameter space by direct simulations is impractical. To reduce the complexity in the parametric analysis, we assume that:

**Assumptions for burst-and-coast gait assembly algorithm**

1) A burst-and-coast bout can be considered as assembled by a burst-phase and a coast-phase. The two phases share the same upper and lower velocity boundaries, while the hydrodynamics of each phase is independent.

2) The swimming performance in the burst phase is determined by the burst kinematics, i.e. the tail-beat frequency and amplitude. A burst phase with velocity range from $U_L$ to $U_U$ can be regarded as a part trimmed from a full burst process with a velocity range from 0 to $U_{terminal}$ ($U_{terminal} \ge U_U$).

3) The coast phase is completely passive. The instantaneous drag only depends on the instantaneous velocity. Thus, an arbitrary coast phase with velocity ranges from $U_U$ to $U_L$ can be regarded as a part trimmed from a coast process covering a sufficiently large velocity range.

4) The transitions between burst- and coast- phases are immediate and can be neglected.

As shown in Fig. 2(e), based on these assumptions, we designed an indirect numerical approach, which uses a database of 41 direct simulations (40 bursts and 1 coast) to assemble arbitrary burst-and-coast gaits. The parameter space was scanned with coarse resolution to find burst-and-coast gaits satisfying the specific velocity requirements, and then further determined the optimal burst-and-coast gait corresponding to minimal CoT.

According to the analysis in sub-section 3.3, the instantaneous velocity during a burst using the tail beat frequency $f$ and amplitude $A$ can be approximated as:

$$U_{\text{burst}}(t) = U_{\text{terminal}}(f, A) \cdot (1 - e^{-6.49t/(A^{-1/2}L^{1/2}f^{-1})}) \quad (\text{Eq.7})$$

where $U_{terminal}(f, A)$ is the terminal velocity of a burst. As to the instantaneous power during the burst, based on the instantaneous power data shown in Fig. 5(c), it is reasonable to approximate the instantaneous power by the terminal power of the burst as:

$$P_{\text{burst}}(t) = P_{\text{terminal}}(f, A) \quad (\text{Eq.8})$$

We mapped $U_{terminal}(f, A)$ and $P_{terminal}(f, A)$ in Fig. 5, so that $U_{terminal}$ and $P_{terminal}$ of any arbitrary burst process can be quantified by interpolation. According to Eqs. 7 and 8, the instantaneous velocity and power can be further calculated.

On the other hand, in coast swimming, the to-be-trimmed coast process was obtained by a single direct simulation, by letting the model fish stop undulating after reaching a velocity higher
than the highest velocity reached by all 40 burst processes. During this coast phase, the body was held straight, and the fish decelerated until the velocity dropped to almost zero. The mechanical power consumption during the coast phase was assumed to be zero.

The full swimming cycle is obtained by concatenating the trimmed burst and coast time sequences (Fig. 2(e)). The procedure is then duplicated to produce a sawtooth-wave time profile of the velocity. For a given set of the four parameters ($f, \alpha, U_U$ and $U_L$), as long as $U_U$ and $U_L$ are within the possible velocity range of a “full burst process”, a unique burst-and-coast swimming gait can be obtained. The average velocity of the generated burst-and-coast swimming gait is defined as $U_{average}$ and calculated numerically by the following equation:

$$U_{average} = \frac{\int_{0}^{t_{B}} U_{burst}(t)dt + \int_{0}^{t_{C}} U_{coast}(t)dt}{t_{B} + t_{C}} \qquad (Eq.9)$$

The cost of transport of the assembled burst-and-coast gait is defined as:

$$\Omega^* = \frac{\int_{0}^{t_{B}} P_{burst}(t)dt}{mg(\int_{0}^{t_{B}} U_{burst}(t)dt + \int_{0}^{t_{C}} U_{coast}(t)dt)} \qquad (Eq.10)$$

In order to find an optimal burst-and-coast swimming gait that would meet the required average velocity $U_{average}$ with the lowest cost of transport, we programmed a MATLAB code to scan the four parameter dimensions ($f$: scan resolution 0.1Hz, scan range 2 ~ 18Hz; $A$: scan resolution 0.0015L, range approximately 0.02~0.26L; the scan resolution in $U_U$ and $U_L$ is less than $10^{-7}$ L/s). To rule out unrealistically short burst-and-coast bouts, in the simulation we require that the burst time is long enough to complete at least one tail-beat stroke, while the coast time is long enough to skip at least one tail-beat stroke, i.e. $T_B \ge 0.5/f$ and $T_C \ge 0.5/f$.

For further details of the scan algorithm, see § C, ESM.

## B. Optimal gait and performance in burst-and-coast swimming

By using the abovementioned gait assembly approach, we specified optimal burst-and-coast gaits with minimal CoT. As shown by the circles (○) in Fig. 6 (a) and (b), the burst-phase frequency of the optimal burst-and-coast gait, represented by $Re_{wave}$, increased with the mean speed, while the optimal peak-to-peak tail-beat amplitude of the optimal burst-and-coast gait appears to be relatively stable at around 0.16L. Between frequency and amplitude, fish in burst-and-coast gait seem to use frequency as the primary control means. Fig. 6(c) demonstrates that the optimal upper- and lower-speed bounds, $U_U$ and $U_L$, respectively, monotonically increase with the average speed. We calculated the ratio of velocity fluctuation as:

$$\delta = \frac{U_U - U_L}{U_{average}} \times 100\% \qquad (Eq.11)$$

At lower speed (1 L/s) $\delta$ was about 27%. In the medium and high-speed regime, $\delta$ became narrower as speed increases, decreasing to about 7% at 5 L/s. Overall, the ratio of velocity fluctuation was less than 30% and tended to become narrower as average speed increases.
### C. Burst-and-coast vs continuous swimming

This sub-section quantitatively compares the optimal kinematics and energetic performance between burst-and-coast and continuous swimming gaits. Here, optimization was also required for the continuous swimming gait for each aimed speed level. The cost of transport of cyclic swimming was calculated as:

$$ \Omega^* = \frac{P_{\infty}(f, A)}{mgU_{\infty}(f, A)} \qquad (\text{Eq.12}) $$

where $P_{\infty}$ and $U_{\infty}$ are respectively the velocity and power of cyclic swimming using the kinematics defined by $(f, A)$, which are respectively equivalent to $P_{\text{terminal}}(f, A)$ and $U_{\text{terminal}}(f, A)$ in Fig. 5(f) and (g). The optimization consisted in finding a minimal $\frac{P_{\infty}(f, A)}{U_{\infty}(f, A)}$ when $U_{\infty}(f, A)$ equals the target speed level, which has been calculated in Li et al. 2021 (Fig. 2 in Li et al. 2021 [14]) and we adopted the results.

The results of the optimal continuous undulation gait are shown in Fig. 6 by rhombus (◆). The predicted optimal frequency of the continuous undulation gait basically raised as the target speed increased (Fig. 6(a)). Compared with the optimal frequency of the burst-and-coast gait, the optimal frequency of the continuous gait was relatively lower. Especially, at the very low velocity of 1 L/s, the continuous gait merely used half the frequency of that used in burst-and-coast gait. The optimal tail-beat amplitude of the continuous gait appeared to be also constant, and very similar to that of the burst-and-coast gait (Fig. 6(b)).

Fig. 6(d) shows the comparison of the predicted optimal CoT values between the burst-and-coast gait and the continuous swimming gait. It is noteworthy that, at all speed levels, the optimized burst-and-coast gait required a lower CoT than that of the continuous undulation gait. In the meantime, our results show that the relative reduction of the CoT by the burst-and-coast gait was most significant at lower speeds, and became less significant as the speed increased (Fig. 6(d)). As the speed increased, the optimal burst-and-coast gait approached the optimal continuous gait in many aspects.

### D. Non-optimized burst-and-coast gaits may be energetically inefficient.

The simulation also shows that non-optimized burst-and-coast gaits result in a rather wide range of CoT values. When scanning the parameter space by the gait-assembly approach, we also obtained extremely large CoT at each velocity levels. Those ‘bad’ kinematics corresponding to extremely high CoT generally occur near the boundary of the parameter space, with both extremely large tail-beat frequency or amplitude (at all velocity levels, maximal-CoT burst-and-coast gaits use extreme kinematic values: maximal $f$ and $A$ in the scanned range).

In Fig. 7 (a), the CoT of non-optimized burst-and-coast gaits can be several times higher than the minimal CoT at each velocity level, especially at lower velocity regimes, where the ratio
between maximum and minimum CoT can be up to almost 8. The maximum CoT of the burst-and-coast gait is also much higher than that of the optimal continuous gait, indicating that when the gait is not sufficiently optimized, the energetical efficiency of burst-and-coast could be more expensive than that of the continuous swimming gait. To provide an example to demonstrate typical “good” and “bad” burst and coast gaits, we present the velocity time sequences of the optimal and the “worst” burst-and-coast gaits at target speed $U_{target} = 3$ L/s in Fig. 7(b). Inefficient burst-and-coast gaits were generally characterized by a wide velocity interval (between lower and upper velocity boundaries) and large bout time, used extremely high tail-beat frequencies and amplitudes in the burst phase. The fish accelerated to terminal speed and then coasted for a long period to fully exhaust the momentum. On the contrary, in the optimal burst-and-coast gait, the velocity interval was narrow and the bout time was much shorter. Optimal burst-and-coast gaits used moderate tail-beat frequency and amplitude, and the velocity time-sequences fluctuated around the average speed.
A. Different understanding from previous investigations

In contrast with most previous studies in the comparison of energetical performance between burst-and-coast and continuous swimming gaits, we have in the present work ensured that: 1) the comparison is based on the same velocity [17], and 2) both burst-and-coast and continuous gaits are respectively optimized parametrically. We consider that these two principles deliver a fair comparison between burst-and-coast and continuous gaits. These principles, to some extent, narrow the difference between the energetical performance between burst-and-coast and continuous gaits, in contrast with previous studies reporting a 50% ~ 60% energy saving by adopting burst-and-coast swimming [9], [18]. It is worth noting that, for both burst-and-coast and continuous gaits, the mechanical CoT is positively correlated to swimming speed. Therefore, when we compare the energetical performance between burst-and-coast and continuous swimming gaits we ensure that both are measured at the same swimming speed. This requirement, together with the carangiform kinematics and moderate Reynolds number led us to a different conclusion from what Xia et al. [21] found in their study that focused on thunniform swimming. Our direct simulation of a self-propelled fish reveals details of the hydrodynamics of burst-and-coast swimming. It shows that the burst-and-coast gait may save energy mainly because the ratio of undulatory-body drag to straight-body drag at a specific velocity is greater than 1 (i.e. $\beta > 1$), which agrees with the understanding of previous analytical models [9]. However, our simulation also demonstrates the complexity of the $\beta$ value, which changes rapidly during acceleration and transition between burst and coast. Especially, $\beta$ is relatively large at the early stage of the burst (Fig. 4(i)), hence estimating $\beta$ as a constant seems inaccurate. Furthermore, the CFD prediction demonstrates that the instantaneous power of the burst phase does not follow a proportional relationship with instantaneous velocity. Instead, it sustains at a constant magnitude determined by the kinematics, which may be approximated by the power of the continuous cyclic swimming under the same undulation kinematics. For both $\beta$ and instantaneous power issues, our CFD simulations show that the burst phase is more energetically expensive than anticipated from previous analytical models.

This study provides insight in the roles of control parameters during the highly unsteady burst process. As shown by the typical time profiles of $U$ and $P$ in Fig.5(a) based on the definition of dimensionless $t^*$, increment in both tail-beat amplitude and frequency can enhance the acceleration toward the terminal velocity (i.e. the rate of $U/U_{terminal}$ approaches 1). Between tail-beat frequency and amplitude, tail-beat frequency plays a relatively stronger role than tail-beat amplitude in acceleration. Li et al. [22] quantitatively showed that the choice of tail-beat amplitude can cause changes in swimming drag, as predicted by the Bone-Lighthill hypothesis, which forces fish to find a proper tail-beat amplitude to maintain the effectiveness in propulsion while avoiding causing excessive drag. Such strategy can extend to the unsteady burst-and-coast. Thus, a fish in burst phase still tends to keep a steady tail-beat amplitude level and adjust tail-beat frequency to control swimming velocity.
## B. Limitations of the current approach

The indirect numerical approach (i.e. the “gait assembly” algorithm), has been used in our previous study to analyze the burst-and-coast gait observed experimentally in the station-keeping swimming of rummy-nose tetrafish [22]. For moderate to high velocities ($U > 1$ BL/s), the parameters that minimize the energy cost of swimming match closely the experimental data, however, the prediction at low-velocity regimes diverges from the observations. The gait assembly algorithm is a compromise approach to explore the four-dimensional parameter space at a feasible cost. The primary source of error in the gait assembly algorithm is the absence of the transition process between the burst and coast phases, while neglecting the fluctuation of the velocity and power is also an important factor of error. The gait assembly algorithm only assembled smoothed velocity trajectories. To examine how the forces and powers differ from direct CFD under the “gait assembly” assumptions, two testing CFD cases are conducted respectively at speeds of 1L/s and 3L/s in §E, ESM. According to this comparison between CFD and the “gait assembly” algorithm, the latter appears, despite its imperfections, to be a reasonable method to explore the highly complex parameter space of the burst-and-coast swimming. To further improve the accuracy and reliability in future works, massive computational resources are needed. Also, beyond the scope of this paper, an optimization algorithm with machine learning that would optimize the strategy parametrically based on both real-time simulation and a historically accumulated database may provide a new perspective. Recent works on the optimization in fish collective swimming based on machine learning [33], [34] provide examples.

Our numerical approach focused on hydrodynamic consumption and neglected physiological consumption. A difference exists between metabolic consumption and hydrodynamic consumption – the former is higher because of the basal metabolic consumption (the rate of energy expenditure at rest is not zero) and the lossy conversion of chemical into mechanical energy [35], [36]. As explained by Li et al. [14], when physiological contributions are included in the consideration, the relation between metabolic consumption and speed is likely U-shaped [36], [37], thus a global optimal swimming velocity for cruising can be found, whereas when physiological contributions are neglected in a computational fluid dynamic investigation, the relation between metabolic consumption and speed is basically monotonic, thus a global optimal swimming velocity cannot be found. Therefore, the current CFD approach can only be used to optimize kinematics at a specified speed, but cannot be used to find a global optimal speed with minimal metabolic cost, unless CFD approaches are combined with models representing the basal metabolic consumption and the conversion of chemical energy into mechanical work by the swimming musculature in future.

## C. Advices for robot fish design

The Bone-Lighthill hypothesis is a basis of the energy saving mechanism of the burst-and-coast swimming. Hence the optimal burst-and-coast gait obtained in this study may only apply to
undulatory (or flapping) swimmers. It remains to be verified if a the burst-and-coast propulsive
mode can be applied to rigid-body swimmers, such as torpedo-shaped autonomous underwater
vehicles (AUVs), where changes in drag between the bursting and coasting phases might be
negligible.

It should be mentioned that, in a realistic situation, the burst-and-coast parameters may
change over the course of a swimming trajectory with the changing environment and navigation
objectives. The gait assembly approach is justified if the time scale of those variations is longer
than the burst-and-coast bout time. Indeed, swimming trajectories consist of approximately
constant gait sequences when fish forage [38] or swim in a stream [22]. The gait assembly
approximation may not be suitable for analyzing rapid maneuvering or strong external
perturbations, but these situations are of interest from the viewpoint of stability and control rather
than energy efficiency.

Our results suggest that, for fish-like robots, the mechanical power during acceleration may
be approximated by the power during continuous cyclic swimming with the same undulation
kinematics. So far, the burst-and-coast style of locomotion can be realized by fish-like robots and
artificial swimmer designs [39], [40]. This research demonstrates that the energy-saving function
in a burst-and-coast gait is feasible. Nonetheless, it also warns that unoptimized burst-and-coast
gaits may be extremely inefficient, with energetical efficiency much worse than the continuous
gait with the same velocity, which partly agrees with a numerical study by Ashraf *et al.* [24].
Therefore, when introducing the burst-and-coast swimming to fish-like robot systems, developers
ought to be aware that burst-and-coast swimming does not necessarily lead to energy saving, and
be cautious when choosing the burst-and-coast gait. As suggested by Fig. 7, a burst-and-coast gait
with short bout time and small velocity interval may be relatively more economical than a burst-
and-coast gait with long bout time and large velocity interval.
<figure><img src="image_2.png" /><figcaption>FIG. 1. Features of burst-and-coast swimming. (a) Example observation on burst-and-coast swimming; (b) Schematic description of burst-and-coast swimming with regards to time and velocity, redrawn from Videler and Weihs [10].</figcaption></figure>
<figure><img src="image_3.png" /><figcaption>FIG. 2. Numerical methods. (Upper) Direct numerical approach: (a) Fish surface model; (b) Kinematic model; (c) Computational mesh; (d) An example of CFD flow field; (e) Indirect numerical approach: Burst-and-coast *gait assembly*.</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>**FIG. 3.** Flow field features of a non-repeated burst-and-coast bout. The burst of cyclic swimming lasts for 15 tail-beat cycles during which the fish has accelerated to a relatively stable velocity. It then stops its undulation and the fish body is kept straight during the coast phase that lasts for another 1.5s (equivalent to a 15 tail-beat cycle time), reaching an almost static condition. Flow features are marked by numbers and explained in the text. Panels (a), (b), (c) and (d) correspond to 5, 10.5, 17 and 30 tail-beat cycles, respectively.</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>FIG. 4. Simulation results of a single burst-and-coast bout. (a) Body axis time sequences; (b) *U*\* vs tail beat cycles; (c) Ω\*<sub>accumulated</sub>; (d) *P*\* and η vs tail-beat cycles. I *P*\* and η vs *U*\*; (f) *T*\* vs tail-beat cycles; (g) *T*\* and *D*\* vs *U*\*; (h) vs tail-beat cycles; (i) Drag coefficient *C*<sub>D</sub>. In (i), as a reference, a speed-specific drag coefficient curve of the fish steadily gliding with a straight body is presented by gray dash line. Note that this</figcaption></figure>
reference curve does not represent a dynamic process. Each point on the reference curve denotes a drag coefficient obtained in a simulation when the fish constantly gliding at a given speed. (U*: dimensionless velocity, defined as $U* = U/U_{ref}$; P*: dimensionless power, defined as $P* = P/(\rho U_{ref}^3 L^2)$; T*: dimensionless thrust, defined as $T* = T/(\rho U_{ref}^2 L^2)$; D*: dimensionless drag, defined as $D* = D/(\rho U_{ref}^2 L^2)$; C_D: drag coefficient, defined as $C_D = 2D/(\rho U^2 S)$; P: instantaneous power; ρ: the water density; U: instantaneous velocity; $U_{ref}$: the reference velocity, defined as the velocity at the end of burst, 3.9L/s; L: the body length; S: the wetted area; T: instantaneous thrust; D: instantaneous drag; P: instantaneous power; η: Froude efficiency, where $η = TU/P$).
<figure><img src="image_6.png" /><figcaption>FIG. 5. Simulation results in 40 full burst simulations with various undulatory kinematics: the combinations between tail beat frequency and amplitude. (a) Relationship between instantaneous velocity and dimensionless time for all burst simulations, where the dimensionless instantaneous velocity in each simulation is normalized by the terminal velocity $U_{\text{terminal}}$ (b) Relationship between instantaneous power and dimensionless time in all burst simulations, where the instantaneous dimensionless power in each simulation is normalized by the terminal power $P_{\text{terminal}}$. The dimensionless time in (a) and (b) is given by $t^* = t/(A^{-1/2}L^{1/2}f^{-1})$. (c) Relationship between $Sw$, $Re_{\text{wave}}$ and $Re_{\text{terminal}}$ in burst simulations, where $Sw = fAL/v$, $Re_{\text{wave}} = fL^2/v$, $A^* = A/L$ and $Re_{\text{terminal}} = UL/v$; (d) Relationship between dimensionless burst-undulation frequency $Re_{\text{wave}}$ and power coefficient ($C_p = P_{\text{terminal}}/(\rho L^3 f^3 A^2)$), the relationship between $Re_{\text{terminal}}$ and $C_p$ is also shown in the inset; (e) $Re_{\text{terminal}}$ map in the</figcaption></figure>
$Re_{wave-A^*}$ plane; (f) $T_{terminal}^*$ map in the $Re_{wave-A^*}$ plane; (g) $P_{terminal}^*$ map in the $Re_{wave-A^*}$ plane. (e), (f) and (g) are recalculated from Li et al. [14]. The dimensional data of all burst simulations are shown in §D, ESM.
<figure><img src="image_7.png" /><figcaption>FIG. 6. Comparison between optimal burst-and-coast and continuous gaits at various speed levels. (a) optimal tail-beat frequency; (b) Optimal tail-beat amplitude; (c) Upper and lower speed boundaries of optimal burst-and-coast swimming, and the fluctuation range around the average speed. Note that upper and lower velocity boundaries are not applicable for the continuous gait; (d) Cost of transport of optimal burst-and-coast gait and continuous gait. 1L/s = 0.02m/s.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>FIG. 7. Optimal gaits at various aim velocity levels with burst-and-coast and continuous gaits, as well as information about the worst burst-and-coast gaits scanned. (a) Optimal cost of transport in burst-and-coast gait (○) and continuous gait (♦), as well as the CoT range between optimal and worst burst-and-coast gaits scanned. (b) Examples of scanned optimal and worst burst-and-coast gaits at $U_{\text{target}} = 3$ L/s. The worst gait adopted an extremely powerful burst to reach a high velocity and then coasted for a long period to exhaust momentum.</figcaption></figure>
## AUTHOR CONTRIBUTIONS

Conceptualization: Gen Li, Dmitry Kolomenskiy, Benjamin Thiria, Ramiro Godoy-Diana.

Data analysis: Gen Li, Dmitry Kolomenskiy, Benjamin Thiria, Ramiro Godoy-Diana.

Methodology: Gen Li, Hao Liu, Dmitry Kolomenskiy

Project administration: Gen Li

Software: Gen Li, Hao Liu.

Validation: Gen Li, Dmitry Kolomenskiy, Benjamin Thiria, Ramiro Godoy-Diana

Visualization: Gen Li, Dmitry Kolomenskiy

Writing – original draft: Gen Li; Revision: all authors

## ACKNOWLEDGEMENTS

Gen Li is funded by the Japan Society for the Promotion of Science (20K14978).
REFERENCES

[1] M. Sfakiotakis, D. M. Lane, and J. B. C. Davies, "Review of fish swimming modes for aquatic locomotion," *IEEE J. Ocean. Eng.*, vol. 24, no. 2, pp. 237-252, 1999, doi: 10.1109/48.757275.

[2] J. R. Hunter and J. R. Zweifel, "Swimming speed, tail beat frequency, tail beat amplitude, and size in jack mackerel, *Trachurus symmetricus*, and other fish," *Fish. Bull. Fish Wildl. Serv. U.S.*, vol. 69, pp. 253-267, 1971.

[3] M. J. Lighthill, "Large-amplitude elongated-body theory of fish locomotion," *Proc. R. Soc. London. Ser. B. Biol. Sci.*, vol. 179, no. 1055, pp. 125-138, 1971, doi: 10.1098/rspb.1971.0085.

[4] J. J. Videler, *Fish swimming*, 10th ed. pringer Science & Business Media., 1993.

[5] D. L. Kramer and R. L. McLaughlin, "The behavioral ecology of intermittent locomotion," *Am. Zool.*, vol. 41, no. 2, pp. 137-153, 2001, doi: 10.1668/0003-1569(2001)041[0137:TBEOIL]2.0.CO;2.

[6] S. P. Windsor, D. Tan, and J. C. Montgomery, "Swimming kinematics and hydrodynamic imaging in the blind Mexican cave fish (Astyanax fasciatus)," *J. Exp. Biol.*, vol. 211, no. 18, pp. 2950-2959, 2008, doi: 10.1242/jeb.020453.

[7] A. P. Soto and M. J. McHenry, "Pursuit predation with intermittent locomotion in zebrafish," *J. Exp. Biol.*, vol. 223, 2020, doi: 10.1242/jeb.230623.

[8] A. McKee, A. P. Soto, P. Chen, and M. J. McHenry, "The sensory basis of schooling by intermittent swimming in the rummy-nose tetra (Hemigrammus rhodostomus): Schooling by intermittent swimming," *Proc. R. Soc. B Biol. Sci.*, vol. 287, no. 1937, 2020, doi: 10.1098/rspb.2020.0568rspb20200568.

[9] D. Weihs, "Energetic Advantages of Burst Swimming of Fish," *J. Theor. Biol.*, vol. 48, no. 1, pp. 215-229, 1974.

[10] J. J. Videler and D. Weihs, "Energetic advantages of burst-and-coast swimming of fish at high speeds.," *J. Exp. Biol.*, vol. 97, pp. 169-178, 1982.

[11] J. J. Videler, "Swimming Movements, Body Structure and Propulsion in Cod Gadus morhua," *Symp. zool. Soc. Lond.*, vol. 48, no. January 1981, pp. 1-27, 1981.

[12] G. Wu, Y. Yang, and L. Zeng, "Kinematics, hydrodynamics and energetic advantages of burst-and-coast swimming of koi carps (Cyprinus carpio koi)," *J. Exp. Biol.*, vol. 210, no. 12, pp. 2181-2191, 2007, doi: 10.1242/jeb.001842.

[13] E. J. Anderson, W. R. McGillis, and M. A. Grosenbaugh, "The boundary layer of swimming fish," *J. Exp. Biol.*, vol. 204, no. 1, pp. 81-102, 2001.

[14] G. Li, H. Liu, U. K. Müller, C. J. Voesenek, and J. L. van Leeuwen, "Fishes regulate tail-beat kinematics to minimize speed-specific cost of transport," *Proc. R. Soc. B Biol. Sci.*, vol. 288, no. 1964, 2021, doi: 10.1098/rspb.2021.1601.

[15] D. Webb, P. W., Weihs, *Fish biomechanics*. Praeger, 1983.

[16] D. Floryan, T. Van Buren, and A. J. Smits, "Forces and energetics of intermittent swimming," *Acta Mech. Sin. Xuebao*, vol. 33, no. 4, pp. 725-732, 2017, doi: 10.1007/s10409-017-0694-3.

[17] L. Dai, G. He, X. Zhang, and X. Zhang, "Intermittent locomotion of a fish-like swimmer driven by passive elastic mechanism," *Bioinspiration and Biomimetics*, vol. 13, no. 5, 2018, doi: 10.1088/1748-3190/aad419.
[18] E. Akoz and K. W. Moored, "Unsteady propulsion by an intermittent swimming gait," *J. Fluid Mech.*, vol. 834, pp. 149–172, 2018, doi: 10.1017/jfm.2017.731.

[19] F. E. Fish, "Swimming strategies for energy economy," in *Fish locomotion. An eco-ethological perspective*, 2010, pp. 90–122.

[20] R. W. Blake, "Functional design and burst-and-coast swimming in fishes.," *Can. J. Zool.*, vol. 61, no. 11, 1983, doi: 10.1139/z83-330.

[21] D. Xia, W. shan Chen, J. kao Liu, and X. Luo, "The energy-saving advantages of burst-and-glide mode for thunniform swimming," *J. Hydrodyn.*, vol. 30, no. 6, pp. 1072–1082, 2018, doi: 10.1007/s42241-018-0120-8.

[22] G. Li et al., "Burst-and-coast swimmers optimize gait by adapting unique intrinsic cycle," *Commun. Biol.*, vol. 4, no. 1, pp. 1–7, 2021, doi: 10.1038/s42003-020-01521-z.

[23] P. Han, J. Wang, and H. Dong, "Effects of intermittent swimming gait in fish-like locomotion," in *AIAA Scitech 2020 Forum*, 2020, vol. 1 PartF, doi: 10.2514/6.2020-1779.

[24] I. Ashraf, S. van Wassenbergh, and S. Verma, "Burst-and-coast swimming is not always energetically beneficial in fish (Hemigrammus bleheri)," *Bioinspiration and Biomimetics*, vol. 16, no. 1, p. 16002, 2021, doi: 10.1088/1748-3190/abb521.

[25] G. Li, U. K. Müller, J. L. Van Leeuwen, and H. Liu, "Body dynamics and hydrodynamics of swimming fish larvae: A computational study," *J. Exp. Biol.*, vol. 215, no. 22, pp. 4015–4033, 2012, doi: 10.1242/jeb.071837.

[26] G. Li, U. K. Müller, J. L. Van Leeuwen, and H. Liu, "Fish larvae exploit edge vortices along their dorsal and ventral fin folds to propel themselves," *J. R. Soc. Interface*, vol. 13, no. 121, 2016, doi: 10.1098/rsif.2016.0585.

[27] H. Liu, "Integrated modeling of insect flight: From morphology, kinematics to aerodynamics," *J. Comput. Phys.*, vol. 228, no. 2, pp. 439–459, 2009, doi: 10.1016/j.jcp.2008.09.020.

[28] C. S. J. J. Wardle, J. Videler, and J. Altringham, "Tuning in to fish swimming waves: body form, swimming mode and muscle function," *J. Exp. Biol.*, vol. 198, pp. 1629–1636, 1995, doi: 10.1242/jeb.198.8.1629.

[29] G. Li, U. K. Müller, J. L. Van Leeuwen, and H. Liu, "Escape trajectories are deflected when fish larvae intercept their own C-start wake," *J. R. Soc. Interface*, vol. 11, no. 101, 2014, doi: 10.1098/rsif.2014.0848.

[30] G. Li, D. Kolomenskiy, H. Liu, B. Thiria, and R. Godoy-Diana, "On the energetics and stability of a minimal fish school," *PLoS One*, vol. 14, no. 8, pp. 1–20, 2019, doi: 10.1371/journal.pone.0215265.

[31] O. Akanyeti, J. Putney, Y. R. Yanagitsuru, G. V. Lauder, W. J. Stewart, and J. C. Liao, "Accelerating fishes increase propulsive efficiency by modulating vortex ring geometry," *Proc. Natl. Acad. Sci. U. S. A.*, vol. 114, no. 52, pp. 13828–13833, 2017, doi: 10.1073/pnas.1705968115.

[32] M. Gazzola, M. Argentina, and L. Mahadevan, "Scaling macroscopic aquatic locomotion," *Nat. Phys.*, vol. 10, no. 10, pp. 758–761, 2014, doi: 10.1038/nphys3078.

[33] S. Verma, G. Novati, and P. Koumoutsakos, "Efficient collective swimming by harnessing vortices through deep reinforcement learning," *Proc. Natl. Acad. Sci. U. S. A.*, vol. 115, no. 23, pp. 5849–5854, 2018, doi: 10.1073/pnas.1800923115.
[34] S. L. Brunton, B. R. Noack, and P. Koumoutsakos, "Machine Learning for Fluid Mechanics," *Annu. Rev. Fluid Mech.*, vol. 52, pp. 477-508, 2020.

[35] P. W. Webb, "The swimming energetics of trout," *J. Exp. Biol.*, vol. 55, pp. 521-540, 1971.

[36] S. P. Gerry and D. J. Ellerby, "Resolving shifting patterns of muscle energy use in swimming fish," *PLoS One*, vol. 9, no. 8, 2014, doi: 10.1371/journal.pone.0106030.

[37] V. Di Santo, C. P. Kenaley, and G. V. Lauder, "High postural costs and anaerobic metabolism during swimming support the hypothesis of a U-shaped metabolism-speed curve in fishes," *Proc. Natl. Acad. Sci. U. S. A.*, vol. 114, no. 49, pp. 13048-13053, 2017, doi: 10.1073/pnas.1715141114.

[38] B. François, "Physical aspects of fish locomotion: an experimental study of intermittent swimming and pair interaction.," Thèse de doctorat de l'Université de Paris, soutenue le 30 mars 2021., 2021.

[39] X. Ye, Y. Su, S. Guo, and L. Wang, "Design and realization of a remote control centimeter-scale robotic fish," 2008, doi: 10.1109/AIM.2008.4601629.

[40] S. Verma, P. Hadjidoukas, P. Wirth, and P. Koumoutsakos, "Multi-objective optimization of artificial swimmers," *2017 IEEE Congr. Evol. Comput. CEC 2017 - Proc.*, pp. 1037-1046, 2017, doi: 10.1109/CEC.2017.7969422.
# Electronic Supplementary Materials

## Intermittent vs. continuous swimming: an optimization tale

Gen Li<fnref n="1" />*, Dmitry Kolomenskiy<fnref n="2" />, Hao Liu<fnref n="3" />, Ramiro Godoy-Diana<fnref n="4" />, Benjamin Thiria<fnref n="4" />

## Contents

Part A Definitions of symbols and abbreviations..................................................................................2

Part B Computational model..................................................................................................................4

B-1 Overview ..................................................................................................................................................4

B-2 Computational grids and fish model deformation........................................................................5

B-3 Body length correction algorithm ....................................................................................................7

B-4 Fluid solution............................................................................................................................................8

B-5 Motion solution........................................................................................................................................9

B-6 Definition of thrust and drag........................................................................................................11

B-7 Difference between amplitude control parameter and resultant tail-beat amplitude............12

B-8 Information of CFD solver validation..................................................................................14

Part C “Gait Assembly” Flow Chat........................................................................................................15

Part D Supplementary Computational Results ..................................................................................16

Part E “Gait assembly” vs CFD ........................................................................................................17

References .........................................................................................................................................................19

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

<footnote n="2">² Center for Design, Manufacturing and Materials (CDMM), Skolkovo Institute of Science and Technology, Moscow, Russia</footnote>

<footnote n="3">³ Graduate School of Engineering, Chiba University, Chiba, Japan</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é de Paris, 75005, Paris, France</footnote>

<footnote>* ligen@jamstec.go.jp</footnote>
# Part A
## Definitions of symbols and abbreviations

<table><caption>**Table S1.** List of symbols and abbreviations.</caption><thead><tr><th>Symbols and abbreviations</th><th>Physical meaning</th></tr></thead><tbody><tr><td>A</td><td>tail-beat amplitude</td></tr><tr><td>α</td><td>tail-beat amplitude control factor</td></tr><tr><td>β</td><td>drag amplifying factor</td></tr><tr><td>C<sub>D</sub></td><td>drag coefficient</td></tr><tr><td>CFD</td><td>computational fluid dynamics</td></tr><tr><td>CoM</td><td>centre of mass</td></tr><tr><td>CoT or Ω</td><td>cost of transport</td></tr><tr><td>D</td><td>cycle-averaged drag; D(t) denotes instantaneous drag</td></tr><tr><td>DoF</td><td>degrees of freedom</td></tr><tr><td>f</td><td>frequency</td></tr><tr><td>L</td><td>length of the fish model, used as reference length</td></tr><tr><td>l</td><td>length</td></tr><tr><td>m</td><td>mass of the fish model</td></tr><tr><td>P</td><td>total mechanical power</td></tr><tr><td>Re</td><td>Reynolds number</td></tr><tr><td>S</td><td>wetted area of the fish model</td></tr><tr><td>t</td><td>time</td></tr><tr><td>T</td><td>cycle-averaged thrust; T(t) denotes instantaneous thrust</td></tr><tr><td>t<sub>p</sub></td><td>period</td></tr><tr><td>U</td><td>swimming speed; U(t) denotes instantaneous swimming speed</td></tr><tr><td>λ</td><td>length of the body wave</td></tr><tr><td>ν</td><td>kinematic viscosity of water</td></tr><tr><td>η</td><td>Froude propulsive efficiency</td></tr><tr><td>ρ<sub>w</sub></td><td>density of water</td></tr></tbody></table>
<table><caption>**Table S2.** Reference values.</caption><thead><tr><td>Parameter</td><td>Value</td></tr></thead><tbody><tr><td>density of water (ρ<sub>w</sub>)</td><td>1.00×10<sup>-3</sup> kg·m<sup>-3</sup></td></tr><tr><td>kinematic viscosity of water (ν)</td><td>8.71×10<sup>-7</sup> m<sup>2</sup>·s<sup>-1</sup></td></tr><tr><td>length of fish (L)</td><td>2.0×10<sup>-2</sup> m</td></tr><tr><td>wetted area of fish model (S)</td><td>1.55×10<sup>-4</sup> m<sup>2</sup></td></tr><tr><td>mass of fish model (m)</td><td>7.9×10<sup>-5</sup> kg</td></tr><tr><td>volume of fish model</td><td>7.9×10<sup>-8</sup> m<sup>3</sup></td></tr></tbody></table>
# Part B
## Computational model

### B-1 Overview

In the simulation, the model fish swims freely in the horizontal plane (3 degrees of freedom (DoF)). The centre-of-mass (CoM) movements and body orientation are not prescribed, but are determined by the hydrodynamic and inertial forces generated by the swimming model fish. The forces acting on the body and motion of the body are obtained through coupling the hydrodynamic and body motion solutions (Fig. S1), which ensures that the motions of the fish correspond to the hydrodynamic and inertial forces exerted on the fish.

<figure><img src="image_9.png" /><figcaption>**Figure S1.** Flow chart of the numerical approach. Black solid lines represent the flow of the computational process, and blue dashed lines show the flow of information.</figcaption></figure>
## B-2 Computational grids and fish model deformation

As shown in Fig. S2, our CFD approach comprises surface models of the changing fish shape (dimension: 121 × 97), and local fine-scale body-fitted grids (dimension: 97×121×20 at $Re<1000$, 97×121×60 at $Re>1000$) plus a large stationary global grid (dimension: various) to calculate the flow patterns around the fish with sufficient resolution.

The global grid surrounded the body-fitted grids and covered a sufficiently large domain to enclose the swimming fish and their wake. The ensemble of body-fitted grids and global grid was set up as a multi-blocked, overset-grid system based on a chimera grid scheme [1]. During the simulation, the body-fitted and global grids share values on their interfaces through inter-grid communication algorithm.

The body was modelled on the silhouette of a Red nose tetra fish (*Hemigrammus bleheri*), with a body length of 2 cm. At this length, the range of the Reynolds number in this study is below 6000, turbulence models are not used, and the grid resolution at $Re=6000$ has been justified in our previous study [2]. All cross-sections of the fish were modelled as ellipses. To reduce the complexity in modelling and computation, we assume that the hydrodynamic influence of all fins other than the tail fin is relatively minor, and neglect them in the model. Also, for the same reason, the gap of the fork-shaped tail fin is neglected, and the fish model has a triangle-shaped fin instead.

We simulated 40 full burst-process simulations, with 5 different tail beat frequency $f = 2, 6, 10, 14$ and $18$ Hz and 8 different amplitude control factor $\alpha = 0.02, 0.06, 0.08, 0.10, 0.12, 0.14, 0.18$ and $0.22$. Note that the amplitude control factor $\alpha$ controls the tail beat amplitude of fish in a reference frame attached to fish head, while the resultant tail beat amplitude in world frame of reference A is solved by simulation. The relation between resultant tail beat amplitude $A$ and tail beat amplitude control factor $\alpha$ is shown in Table S3 and Fig. S5.
<figure><img src="image_10.png" /><figcaption>**Figure S2.** Computational fluid dynamics model. (a) Surface model of Red nose tetra fish (dimension: 97×121); (b) A function drives the instantaneous body shape. Variation of body length caused by this driving function was corrected to keep lateral excursion and body length constant at 1L. (c) Multi-blocked computational grid system consists of local fine-scale body-fitted grid (dimension: 97×121×20 at Re<1000, 97×121×60 at Re>1000) plus a large stationary global grid (dimension: variant).</figcaption></figure>
## B-3 Body length correction algorithm

The sinusoidal functions driving the fish deformation causes the total body length along the midline to vary during the tail beat. This variation is corrected at every time step by a procedure that preserves the lateral excursion while ensuring that the body length remains constant. Fig. S3 explains the procedure of this correction:

An axis (green axis in Fig. S3) is first generated according to the sinusoidal function. The length of this axis exceeds the actual body length. Then, in the longitudinal direction the green axis is linearly contracted until it reaches the actual length, while the lateral excursion is preserved. The corrected axis is used as the body shape of the fish in the simulation at this time step.

<figure><img src="image_11.png" /><figcaption>**Figure S3.** Body length variation caused by deformation is corrected by a procedure that preserves the lateral excursion while ensuring a constant body length.</figcaption></figure>
## B-4 Fluid solution

The governing equations for the fluid solution are the three-dimensional, incompressible and unsteady NS equations written in strongly conservative form for mass and momentum [3]. To accelerate the computation and improve the robustness during iteration, the artificial compressibility method is adopted by adding a pseudo time derivative of pressure to the continuity equation [3]. The governing equations for an arbitrary deformable control volume V are

$$ \int_V \left( \frac{\partial \mathbf{Q}}{\partial t} + \frac{\partial \mathbf{q}}{\partial \tau} \right) dV + \int_V \left( \frac{\partial \mathbf{F}}{\partial x} + \frac{\partial \mathbf{G}}{\partial y} + \frac{\partial \mathbf{H}}{\partial z} + \frac{\partial \mathbf{F}_v}{\partial x} + \frac{\partial \mathbf{G}_v}{\partial y} + \frac{\partial \mathbf{H}_v}{\partial z} \right) dV = 0 $$

where

$$ \mathbf{Q} = \begin{bmatrix} u \\ v \\ w \\ 0 \end{bmatrix}, \quad \mathbf{q} = \begin{bmatrix} u \\ v \\ w \\ p \end{bmatrix}, \quad \mathbf{F} = \begin{bmatrix} u^2 + p \\ uv \\ uw \\ \lambda u \end{bmatrix}, \quad \mathbf{G} = \begin{bmatrix} vu \\ v^2 + p \\ vw \\ \lambda v \end{bmatrix}, \quad \mathbf{H} = \begin{bmatrix} wu \\ wv \\ w^2 + p \\ \lambda w \end{bmatrix}, $$

$$ \mathbf{F}_v = -\frac{1}{Re} \begin{bmatrix} 2u_x \\ u_y + v_x \\ u_z + w_x \\ 0 \end{bmatrix}, \quad \mathbf{G}_v = -\frac{1}{Re} \begin{bmatrix} v_x + u_y \\ 2v_y \\ v_z + w_y \\ 0 \end{bmatrix}, \quad \mathbf{H}_v = -\frac{1}{Re} \begin{bmatrix} w_x + u_z \\ w_y + v_z \\ 2w_z \\ 0 \end{bmatrix}, $$

and u, v, and w are velocity components respectively in X-, Y- and Z-directions of the Cartesian coordinate system; velocity components with subscripts x, y, and z are velocity gradients respectively in X-, Y-, and Z-directions of the Cartesian coordinate system; p is the pressure; λ is the pseudo-compressibility coefficient; the set of equations modified from the incompressible Navier-Stokes equations can be solved implicitly by marching in pseudo time: t denotes the physical time, τ denotes the pseudo time, and the term q associated with the pseudo time is designed for an inner-iteration inside each physical time step, and will vanish when the divergence of velocity is driven to zero so as to satisfy the equation of continuity. By introducing the Reynolds transport theorem and employing the Gauss integration theorem, an integrated form of the Navier-Stokes equations in general curvilinear coordinate system corresponding to the FVM structural mesh is gained as

$$ \frac{\partial}{\partial t} \int_V \mathbf{Q} dV + \int_V \frac{\partial \mathbf{q}}{\partial \tau} dV + \oint_S (\mathbf{f} - \mathbf{Q} \mathbf{U}_V) \cdot \mathbf{n} dS = 0, $$

where $\mathbf{f} = (\mathbf{F} + \mathbf{F}_v, \mathbf{G} + \mathbf{G}_v, \mathbf{H} + \mathbf{H}_v)$; the control volume V is a hexahedral cell, while S denotes the surfaces of the hexahedral cell; $n = (n_x, n_y, n_z)$ is the unit outward normal vector corresponding to all the surfaces of the hexahedral cell; $\mathbf{U}_V$ is the local velocity of the moving cell surface caused by displacement and deformation of the cell. For a structured, three-dimensional mesh system (ξ-, η-, and ζ-dimensions respectively represent each of the three dimensions of mesh) and cell-centered storage architecture, A semi-discrete form can be further derived:

$$ \frac{\partial}{\partial t} (\mathbf{Q} \mathbf{V})_{ijk} + \left( \frac{\partial \mathbf{q}}{\partial \tau} \right)_{ijk} V_{ijk} + \mathbf{R}_{ijk} = 0, $$

where
$$ \mathbf{R}_{ijk} = (\hat{\mathbf{F}} + \hat{\mathbf{F}}_v)_{i+\frac{1}{2},j,k} - (\hat{\mathbf{F}} + \hat{\mathbf{F}}_v)_{i-\frac{1}{2},j,k} + (\hat{\mathbf{G}} + \hat{\mathbf{G}}_v)_{i,j+\frac{1}{2},k} - (\hat{\mathbf{G}} + \hat{\mathbf{G}}_v)_{i,j-\frac{1}{2},k} \\ + (\hat{\mathbf{H}} + \hat{\mathbf{H}}_v)_{i,j,k+\frac{1}{2}} - (\hat{\mathbf{H}} + \hat{\mathbf{H}}_v)_{i,j,k-\frac{1}{2}} $$

$$ \hat{\mathbf{F}} + \hat{\mathbf{F}}_v = (\mathbf{f} - \mathbf{Q}\mathbf{U}_v) \cdot \mathbf{S}_n, \quad \hat{\mathbf{G}} + \hat{\mathbf{G}}_v = (\mathbf{f} - \mathbf{Q}\mathbf{U}_v) \cdot \mathbf{S}_n, \quad \hat{\mathbf{H}} + \hat{\mathbf{H}}_v = (\mathbf{f} - \mathbf{Q}\mathbf{U}_v) \cdot \mathbf{S}_n, \\ \mathbf{S}_n = [S_{nx}, S_{ny}, S_{nz}] = \mathbf{n} \cdot \mathbf{S} $$

i, j, and k are cell indexes respectively in ξ-, η-, and ζ-dimensions in the boundary-fitted curvilinear coordinate system, while “+ 1/2” and “- 1/2” denote the surface location (e.g. i + 1/2 denotes a surface locates on the positive ξ-direction of the hexahedral cell V); $\hat{\mathbf{F}}$, $\hat{\mathbf{G}}$ and $\hat{\mathbf{H}}$ denote inviscid flux in ξ-, η-, and ζ-directions, respectively; $\hat{\mathbf{F}}_v$, $\hat{\mathbf{G}}_v$ and $\hat{\mathbf{H}}_v$ denote viscous flux in ξ-, η-, and ζ-directions, respectively; For each surface of the hexahedral cell, $\mathbf{S}_n$ denotes a vector consisting of projected areas of present surface in X-, Y-, and Z-directions, respectively; $\mathbf{n}$ is the unit outward normal vector of the present surface in the Cartesian coordinate system, and $\mathbf{S}$ is the area of the present surface.

The computational mesh consists of two grids as shown in Fig. S2: one is surrounding the fish body, deforming and moving with the fish and constructed in a body-fitted curvilinear coordinate system (ξ-, η-, and ζ-dimensions are circumferential, axial and radial directions relative to the fish body, respectively); the other is a static background grid, of which the ξ-, η-, and ζ-dimensions, respectively, overlap the X-, Y-, and Z-dimensions of the Cartesian coordinate system.

The NS equations are solved in each grid (block). The boundary conditions for the NS equations are: 1) in the fish-body-fitted grid, the non-slip condition is applied to the innermost cells on the surface of the fish body; 2) at the outside boundaries of the background grid, a zero-gradient condition is used; 3) at the interfaces between fish-body-fitted and background grids, the two grids provide boundary conditions to each other through interpolations.

## B-5 Motion solution

Initially, the model fish was assumed to be of uniform density with a density equal to water. Based on the body surface grid, the body was divided into inner-body cells. The mass of each inner-body cell was calculated when the fish was straight, and it was assumed to be constant during swimming regardless of cell deformation. This assumption ensured that fish mass remained constant throughout the simulated swimming episode.

The centre of mass (CoM) was refreshed at each time step during the computation:

$$ \mathbf{r}_{\text{CoM}} = \frac{\sum_{i}^{\text{body}} (m_i \cdot \mathbf{r}_i)}{\sum_{i}^{\text{body}} m_i} $$
where $m_i$ is the initial mass of the i-th inner-body cell, $r_i$ is the position vector of the centre of the i-th inner-body cell. The moment of inertia of the j-th segment about the segmental vertical centre line as a function of time was obtained by

$$J_{\text{seg},j} = \sum_{i}^{\text{seg}} (m_i \cdot (\mathbf{r}_i - \mathbf{r}_{\text{seg},j})^2)$$

where $\mathbf{r}_i$ is the position vector of the centre of the i-th inner-body cell in this section, and $\mathbf{r}_{\text{seg},j}$ is the central vertical axis through segment j.

The moment of inertia of segment j with respect to the CoM was computed as:

$$J_{\text{CoM},j} = J_{\text{seg},j} + m_{\text{seg},j} \cdot (\mathbf{r}_{\text{CoM}} - \mathbf{r}_{\text{seg},j})^2$$

where $m_{\text{seg},j}$ is the mass of segment j.

The time-dependent instantaneous moment of inertia of the body with respect to the instantaneous location of the CoM was computed as:

$$J_{\text{body}} = \sum_{j}^{\text{body}} J_{\text{CoM},j} = \sum_{i}^{\text{body}} (m_i \cdot (\mathbf{r}_i - \mathbf{r}_{\text{CoM}})^2)$$

where $m_i$ is the initial mass of the i-th inner-body cell, $r_i$ is the position vector of the centre of the i-th inner-body cell.

We define angle $\alpha_j$ for body segment j as

$$\alpha_j = \begin{cases} \arctan\left(\frac{Y_{\text{seg},j}}{X_{\text{seg},j}}\right) & \text{for } X_{\text{seg},j} > 0 \\ \arctan\left(\frac{Y_{\text{seg},j}}{X_{\text{seg},j}}\right) + \pi & \text{for } X_{\text{seg},j} < 0 \end{cases}$$

where $\mathbf{r}_{\text{seg},j} = (X_{\text{seg},j}, Y_{\text{seg},j})$ are the central coordinates of segment j with respect to the location of the CoM. This definition provides a four-quadrant arctan function, with $-\frac{\pi}{2} \le \alpha_i \le \frac{3\pi}{2}$, without a jump in $\alpha_i$ value. Following the derivation in [4], the body angle is defined as:

$$\alpha_{\text{body}} = \sum_{j}^{\text{body}} \frac{\alpha_j \cdot J_{\text{CoM},j}}{J_{\text{body}}} - \alpha_{0,\text{body}}$$

where $\alpha_{0,\text{body}}$ is the body angle in the earth frame of reference at the beginning of simulation, defined as zero in this study.
## B-6 Definition of thrust and drag

As shown in Fig. S4, the instantaneous thrust at each time-step was defined as the sum of the forward components of both pressure and shear force over all fish surface elements, while the instantaneous drag was defined as the sum-up of backward components of both pressure and shear force over all fish surface elements. Hence, each surface element could contribute to thrust or drag at different time steps even within one tail beat cycle, as well as to thrust and drag simultaneously (e.g. generating pressure-based thrust and shear-based drag).

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

**Figure S4.** Flow chart of the computational procedure calculating instantaneous thrust and drag. Pressure and shear stress force components in each surface element were calculated and evaluated if they contribute to thrust or drag, respectively, and added up accordingly. Here, $n_{\text{pres}}$, $n_{\text{shear}}$ and $X$, respectively, denote the unit direction vectors of pressure force, shear stress force and unit direction vector in the forward direction. $S_{\text{pres}}$ and $S_{\text{shear}}$ respectively denote the magnitude of pressure force and shear force in one surface element.
B-7 Difference between amplitude control parameter and resultant tail-beat amplitude

The amplitude control parameter or curvature control parameter is used to manipulate the fish deformation. However, this deformation is defined in a reference frame attached to fish head. When the fish deforms according to the prescribed input curvature, it will experience lateral recoil owing to the free-swimming CFD model, and the resultant tail-beat amplitude in the world frame of reference will be smaller than what we defined in the fish-head frame of reference, as shown in Fig. S5.

Therefore, while the resultant tail-beat frequency is solely determined by the input *f*, the resultant tail-beat amplitude *A* is not solely determined by amplitude/curvature control parameter, but is also affected by frequency. In other words, we define the model fish tail motion in a reference frame attached to the fish, while the resultant fish tail motion ('emerging' from the simulation) in a world reference frame is narrower, and influenced by tail-beat frequency. The relation between tail beat amplitude and tail beat amplitude control factor is shown in Table S3.
<figure><img src="image_13.png" /></figure>

**Figure S5.** We defined the model fish tail motion in a reference frame attached to the fish, while the resultant fish tail motion ('emerging' from the simulation) in a world reference frame is influenced by recoil (i.e., overall body rotation) and becomes narrower.

**Table S3.** The relation between tail beat amplitude *A* and tail beat amplitude control factor *α*. Note that *A* basically depends on *α* while slightly influenced by tail beat frequency *f*.

<table><thead><tr><th rowspan="2"></th><th rowspan="2"></th><th colspan="5">frequency (Hz)</th></tr><tr><th>2</th><th>6</th><th>10</th><th>14</th><th>18</th></tr></thead><tbody><tr><td rowspan="6">α</td><td>0.02</td><td>0.023</td><td>0.023</td><td>0.023</td><td>0.023</td><td>0.023</td></tr><tr><td>0.06</td><td>0.069</td><td>0.070</td><td>0.070</td><td>0.070</td><td>0.071</td></tr><tr><td>0.08</td><td>0.089</td><td>0.093</td><td>0.096</td><td>0.099</td><td>0.100</td></tr><tr><td>0.10</td><td>0.118</td><td>0.123</td><td>0.127</td><td>0.131</td><td>0.132</td></tr><tr><td>0.12</td><td>0.148</td><td>0.154</td><td>0.157</td><td>0.160</td><td>0.162</td></tr><tr><td>0.14</td><td>0.171</td><td>0.184</td><td>0.190</td><td>0.194</td><td>0.197</td></tr><tr><td></td><td>0.18</td><td>0.232</td><td>0.244</td><td>0.252</td><td>0.258</td><td>0.263</td></tr><tr><td></td><td>0.22</td><td>0.287</td><td>0.309</td><td>0.322</td><td>0.331</td><td>0.337</td></tr></tbody></table>

↑Amplitude, unit: *L*, 1*L*=0.02m
B-8 Information of CFD solver validation

<table><caption>**Table S4.** Locations of additional information regarding the computational approach.</caption><thead><tr><th></th><th>specific information</th><th>relevant publication</th></tr></thead><tbody><tr><td rowspan="5">validations</td><td>mesh density independence test</td><td>[2] (in supplementary materials)<br>Note: in this study the resolution based on<br>grid density is higher than that in [2]</td></tr><tr><td>mesh size independence test</td><td>[5] (in supplementary materials)</td></tr><tr><td>validation on hydrodynamic solution<br>on an oscillating cylinder, compared<br>with experiment</td><td>[2] (in supplementary materials)</td></tr><tr><td>validation on flow field<br>on swimming fish, compared with<br>PIV</td><td>[5,6]</td></tr><tr><td>validation on motion solution<br>on swimming fish, compared with<br>experiment</td><td>[2,6]</td></tr></tbody></table>
Part C “Gait Assembly” Flow Chat

<figure><img src="image_14.png" /><figcaption>**Figure S6.** "Gait Assembly" Flow Chat.</figcaption></figure>
Part D Supplementary Computational Results

Full burst process results

<figure><img src="image_15.png" /><figcaption>**Fig. S7** Simulated time-sequences of speed and power, respectively in 40 full burst processes and one coast process. *f* is tail beat frequency and *α* is the tail beat amplitude control factor, the non-linear relation between amplitude *A* and *α* is shown in Table S3</figcaption></figure>
# Part E “Gait assembly” vs CFD

The “gait assembly” approach bases on several simplifying assumptions. To examine how the forces and powers differ from direct CFD under the assumptions, two testing CFD cases are conducted respectively at speeds of 1L/s and 3L/s.

In each CFD case, the gait parameters predicted by “gait assembly” approach were applied to the CFD fish. In “gait assembly” approach, there are four control parameters: the tail-beat frequency *f* during burst swim, the tail-beat amplitude *A* during burst swim, an upper velocity boundary *U*<sub>U</sub> and a lower velocity boundary *U*<sub>L</sub>. However, since *U*<sub>U</sub> and *U*<sub>L</sub> refer to “low-pass filtered speeds”, which ought to be calculated afterwards and not applicable in real-time simulation, in CFD we use bout time (duration of the burst-and-coast swimming cycle) and duty cycle (the ratio of burst time to bout time) instead of *U*<sub>U</sub> and *U*<sub>L</sub> [7]. The values of the four control parameters in the two cases are shown in Table S4.

<table><caption>**Table S4.** The four gait control parameters in the two CFD cases.</caption><thead><tr><th rowspan="2">CFD cases &<br>target speed</th><th colspan="4">CFD setting corresponding to "gait assembly" results</th></tr><tr><th>frequency (Hz)</th><th>amplitude control factor</th><th>bout time (s)</th><th>duty cycle</th></tr></thead><tbody><tr><td>Case 1: 1L/s</td><td>7.1</td><td>0.118</td><td>0.527</td><td>14.6%</td></tr><tr><td>Case 2: 3L/s</td><td>9.5</td><td>0.125</td><td>0.174</td><td>51.7%</td></tr></tbody></table>

In CFD Case 1, the fish started burst-and-coast swim from static condition and reached a balanced average speed of 0.91L/s (see Fig. S8(a) and Table S5). As shown in Table S5, the CFD mean power and dimensionless power are greater than that of the "gait assembly" results.

In CFD Case 2, the fish started burst-and-coast swim from static condition and reached a balanced average speed of 3.21L/s (see Fig. S8(c) and Table S5). As shown in Table S5, the CFD mean power and dimensionless power are again greater than those of the "gait assembly" results.

The comparison between CFD speed curve and "gait assembly" predicted speed curve are respectively shown by Fig. S8(b) and (d). The CFD predictions of the mean speed in both cases were accurate, however, the speed fluctuation intervals predicted by the CFD were larger than that predicted by the "gait assembly". The CFD predicted power was greater than that predicted by "Gait Assembly", we notice that much of the CFD power was wasted in sudden deformation at the transition between burst and coast phase. The CFD power consumption may be potentially reduced when the transitional deformation is optimized or when the fish body is elastic.

Overall, these tests show that the "gait assembly" approach predictions basically agree with direct CFD simulation results.
<table><caption>**Table S5.** Comparison between CFD and "gait assembly" results.</caption><thead><tr><th colspan="3">"gait assembly" results</th><th colspan="3">CFD results</th></tr><tr><th>target speed (L/s)</th><th>mean power (W)</th><th>dimensionless power</th><th>Case</th><th>mean speed (L/s)</th><th>dimensionless power</th></tr></thead><tbody><tr><td>1</td><td>2.32E-07</td><td>7.25E-02</td><td>Case 1</td><td>0.91</td><td>2.63E-07</td><td>1.08E-01</td></tr><tr><td>3</td><td>2.17E-06</td><td>2.51E-02</td><td>Case 2</td><td>3.21</td><td>2.76E-06</td><td>2.61E-02</td></tr></tbody></table>

<figure><img src="image_16.png" /><figcaption>**Table S8.** Comparison between CFD and "gait assembly" speed results. (A) Speed curve of CFD Case 1; (B) Comparison of speed curves between CFD and "gait assembly" at 1L/s; (C) Speed curve of CFD Case 2; (D) Comparison of speed curves between CFD and "gait assembly" at 3L/s.</figcaption></figure>
References

[1] N. C. Prewitt, D. M. Belk, and W. Shyy, *Parallel Computing of Overset Grids for Aerodynamic Problems with Moving Objects*, Progress in Aerospace Sciences.

[2] G. Li, U. K. Müller, J. L. Van Leeuwen, and H. Liu, *Escape Trajectories Are Deflected When Fish Larvae Intercept Their Own C-Start Wake*, J. R. Soc. Interface **11**, (2014).

[3] H. Liu, *Integrated Modeling of Insect Flight: From Morphology, Kinematics to Aerodynamics*, J. Comput. Phys. **228**, 439 (2009).

[4] J. L. van Leeuwen, C. J. Voesenek, and U. K. Müller, *How Body Torque and Strouhal Number Change with Swimming Speed and Developmental Stage in Larval Zebrafish*, J. R. Soc. Interface **12**, 20150479 (2015).

[5] G. Li, U. K. Müller, J. L. Van Leeuwen, and H. Liu, *Fish Larvae Exploit Edge Vortices along Their Dorsal and Ventral Fin Folds to Propel Themselves*, J. R. Soc. Interface **13**, (2016).

[6] G. Li, U. K. Müller, J. L. Van Leeuwen, and H. Liu, *Body Dynamics and Hydrodynamics of Swimming Fish Larvae: A Computational Study*, J. Exp. Biol. **215**, 4015 (2012).

[7] G. Li, I. Ashraf, B. François, D. Kolomenskiy, F. Lechenault, R. Godoy-Diana, and B. Thiria, *Burst-and-Coast Swimmers Optimize Gait by Adapting Unique Intrinsic Cycle*, Commun. Biol. **4**, 1 (2021).