The fluid domain, governed by the Naiver–Stokes equations, was solved using the finite volume method. The incompressible flow governing equations for continuity and momentum are given in Eqs. 1, 2, with u → = u , v the fluid velocity vector, p fluid pressure, μ fluid kinematic viscosity, and ρ fluid density. ∇ ∙ u → = 0 , (1) ∂ u → ∂ t + u → ∙ ∇   u → = − 1 ρ ∇ p + μ ρ ∇ 2   u . → (2)

The transient equations are computed by a pressure–velocity coupling using a non-iteration time-advancement (NITA) scheme, which allows direct solutions with only one outer iteration and sub-iterations to support convergence. The momentum and continuity equations are decoupled in the selected NITA pressure-based solver fractional-step method (FSM). A relaxation factor was applied. For the transient terms, a first-order implicit time marching scheme was selected. Spatial discretisation and discretisation of the diffusive term were achieved by applying the least squares cell-based approach for the gradient and second-order upwind scheme. To improve accuracy, second-order pressure interpolation was selected.

The unstructured CFD mesh dynamically adjusts to the global motion of each surface node defined by the body motion of the swimmer. For a continuous flexible body, the lateral position of the surface points were linearly interpolated between the local coordinate system of the current and ascending body. To maintain a high mesh quality throughout the computational domain, the dynamic meshing function available in Ansys Fluent was selected. For small displacements, the mesh was adjusted by diffusion-based smoothing. For large displacements, remeshing of the unstructured mesh was applied. The unsteady time step and mesh sizing were selected based on the Courant convergence criteria.

The multi-body algorithm describes the internal and external dynamics of a multi-body system by solving the generalised equation of motion in the Newton–Euler form shown in Eq. 3 via three recursive loops through the elements of the multi-body system. The components of Eq. 3 are the generalised force vector β 0 * , including fluid forces and inertial forces, the generalised mass matrix M 0 * , and the acceleration vector η ₀ of the reference body B ₀. Thus β 0 * = η ˙ 0 M 0 * . (3)

A simplified fish geometry was considered a multi-body system in the shape of a 2D NACA0012 foil comprising body segments connected via actuated links along the central line that resemble muscle actuation. The NACA0012 geometry was separated into 10 equal cord length body segments with nine joints.

The local coordinate system B 0 was the starting point of the recursive loop. The undulation of the body was achieved through rotational motions around the joints (including a sinusoidal motion for undulation and an offset to curve the body with respect to the central line) shown in Figure 3. The length of the individual segments l i is evenly distributed and fixed, resulting in a central line length of L = 0.1   m .

The rotational motion of a segment around its front linkage at location s is given by r s , t = A s sin − 2 π   f t + φ s + C s , (4) where f is the undulation frequency and t is the time. The amplitude envelope A s and phase φ s are defined by A s = c 1 s L + c 2 c 3 , (5) φ s = s L , (6) with c 1 , c 2 , a n d   c 3 being kinematic coefficients. The offset term C s was introduced to achieve the turning manoeuvre. The present study includes two curvature envelopes. First, a constant curvature envelope with a uniform offset at all joints and magnitude coefficient c 4 is described by C s = 0.15   c 4 . (7)

Second, a linear curvature envelope from the leading to trailing edge is described by C s = 0.15   c 5 s L + c 6 , (8) where coefficients c 5 and c 6 define the slope and offset of the linear envelope respectively.

The total curvature of the swimmer is defined as the sum of the time-independent joint displacement over the swimmer length, which is as follows κ = ∑ i = 0 9 r i ¯   / l i . (9)

Tables 1, 2 summarise the total curvature and coefficient values used in an investigation into the curvature magnitude in Section 5.1 and comparison of envelopes in Section 5.2 respectively. The coefficients c 5 and c 6 were chosen for linearly increasing and decreasing envelopes to observe performance differences of predominantly head or tail curvings.

To investigate the turning performance at different forward velocities while maintaining a horizontal heading of the swimmer before turning, two linear feedback controllers were employed. Linear and angular swimming velocities are determined by the swimmer’s undulation amplitude and curvature. Therefore, to implement speed and heading control, control variables c a and c S are added to Eq. 4 so that the joint motion is described by r s , t = c a A s sin − 2 π   f t + φ s + c S C s . (10)

Abrupt changes to c a and c S will lead to unintended chaotic motion. Instead, a cosine-based transition function which enables smooth transitions between different states over one oscillation period ∆ t l = T = 1 f was employed. Let c c o n t r o l represent either c a or c s , we have c c o n t r o l = c t l − 1 − c t l − 1 − c t l 2 1 − cos π t 0,1 , (11) where c t l − 1 and c t l are the current control variables sampled at the beginning of the previous transition period and the new control variable sampled at the beginning of the current transition period respectively. t 0,1 runs correspondingly to the general time step but resets to 0 when reaching 1.

Equation 10 does not require a global reference to produce a travelling wave function; hence, it can be considered an open loop with control inputs c a and c s . To enable setpoint tracking for speed control and steering towards waypoints, a closed-loop PID controller was applied to calculate c a and c s as functions of the error defined as e t n = f ∑ k = 0 n e t k ∆ t k , (12) with e t k being the control error at sampling time. To remove periodicity, the signal was averaged over the oscillation cycle T. Based on the time discrete solution of a CFD simulation, the PID controller was implemented in a time discrete recursive form given by Eqs. 13, 14. u t k = u t k − 1 + ∆ u t k , (13) ∆ u t k = K p e t k − e t k − 1 + K I T s e t k + K d T s e t k − 2 e t k − 1 + e t k − 2 . (14)

There are three individual time steps. First, the numerical simulation time step denoted as ∆ t s = t s + 1 − t s . Second, the controller sample time ∆ t k = t k + 1 − t k . Third, the controller update interval at ∆ t l = T = 1 f . The periods were defined so that ∆ t s ≤ ∆ t k < ∆ t l .

Two separate controllers for speed and steering control were used. The first PID controller controlled the speed by varying the undulation amplitude. The control coefficient c a is calculated based on the error between a speed setpoint and the swimmer’s velocity. The swimmer’s heading directed velocity is calculated as the root square of the cycle-averaged global velocities in X e and Y e , as shown in Eq. 15. V s w i m = V ¯ x 2 + V ¯ y 2 . (15)

The control action is limited to 0 ≤ c a ≤ 1 .

To achieve tracking of a set velocity V s e t , we defined the control error at sampling time as follows e v e l t k = V s e t − V s w i m . (16)

In this study, V_set is selected to achieve the desired Reynolds numbers Re = 2,000, 1,500 and 1,000; thus, V_set = 0.02, 0.015 and 0.01 m/s respectively.

The second PID controller adjusted the curvature magnitude by calculating the control coefficient c S based on the error between the line of sight (LOS) angle and the heading angle. LOS navigation was used to provide a reference angle θ between the swimmer’s current position and a waypoint. Using the origins of the local coordinate systems B 0 x B 0 , y B 0 and B 1 x B 1 , y B 1 together with the coordinates of the waypoint ( x w p , y w p ) (all measured in the global frame), θ can be calculated as follows θ = atan y B 1 − y w p x B 1 − x w p , (17) and the heading of the fish α was calculated by α = atan y B 1 − y B 0 x B 1 − x B 0   . (18)

To achieve LOS waypoint tracking the control error is defined as follows e L O S t k = θ − α   . (19)

The control action was limited to − 0.2 ≤ c s ≤ 0.2 .

Both speed and steering controllers were tuned to achieve a critically damped response, i.e., control gains K P , K I , a n d   K D according to Eq. 14 were chosen manually to achieve a fast-converging system response without significant overshoot. Initial control gains were found in [5]. The final used control gains for speed and steering were as follows: K P = 5 , K I = 5 , K D = 55 and K P = 0.03 , K I = 0 , K D = 0.1 . These tuning parameters further avoid unstable fluctuations; consequently, the controller can be considered stable according to the limit cycle behaviour [24]. It also demonstrates that reduced overshoot and faster convergence can be achieved by choosing proper control gains.

As shown in Figure 4, the computational domain is 25 L by 8 L, in which the swimmer’s leading edge is originally located at 5 L by 4 L from the bottom right corner. Herein, L is the fish body length. The unstructured CFD mesh sizing is ∆ x y = 1 3333 at the swimmer boundary and increases to ∆ x y = 1 33 in the far field. The time step was set to ∆ t s = T / 250 , where T is the undulating period. Coefficients of the amplitude function of Eq. 5 were chosen to follow an anguilliform pattern described in the literature [25] with coefficients c 1 = 0.125 , c 2 = 0.0315 a n d   c 3 = 1.03125 . Initial control coefficients were chosen as c a t s = 0 = 0 , c a t s = 1 = 0.1 and c s = 0.0 for t s < 2 T . The velocity inlet condition is set to zero so that the swimmer starts in still water.

The Froude efficiency describes the useful output power over the input energy and is often used to assess engineering systems. With respect to the simplified BCF fish swimming scenario, it describes the relation between energy required to realise the body undulation and the swimming work. Herein, the swimming work is defined as the product of a longitudinal force on the body and swimming velocity. For a fish following a straight trajectory at constant speed, i.e., moving at a quasi-steady state, the cycle-averaged longitudinal force is zero due to the balance between the thrust and drag. The resulting Froude efficiency is zero and therefore not applicable. First mentioned as Self-Propelled Fitness [26] and later proposed as an efficiency measure [27], the quasi-propulsive efficiency is defined as η Q P = R U s / P i n ¯ , where the output power of a swimming fish is approximated by the product of the swimming speed U s and the resistance force R on a rigid body towed at the same speed. Alternatively, the cost of travel (CoT), the ratio of the energy spent per unit distance travelled, is often used in life science as a quantitative measure of energetics in locomotion. The performance of turning fish was evaluated based on the minimum turning radius derived from the trajectory of the centre of mass [6]. Importantly, for an identical centre of mass trajectories the required turning space may differ depending on the swimmers’ body flexibility [28].

The linear manoeuvrability number (LMN), first mentioned to assess the manoeuvrability of an on-land hexapod [29], has also been used to quantify fish manoeuvrability [6]. The LMN is defined as the ratio of the time integral of the force impulse perpendicular to the forward momentum.

Ideally, a quantitative measure of manoeuvrability must provide a combined assessment of power consumption as well as linear and angular displacement. In that sense, CoT and LMN are not suitable. In this study, a modified cost of travel function for turning is applied. The cost of manoeuvring (CoM) is defined as the ratio of the cycle averaged input power to global angular velocity C o M = P i n ¯ ω g l o b a l = P i n ¯ U s , ⊥ r t . (20)

The average input power was calculated as the sum of all joints’ cycle-averaged power (defined as the product of cycle-averaged torque and angular velocity), which is given as follows P i n ¯ = ∑ j = 0 n τ j ¯ r ˙ j ¯ . (21)

The global angular velocity was calculated using the radius r t of the turning trajectory and perpendicular forward velocity U s , ⊥ . The turning radius was calculated using the MATLAB function [30]. To filter out the curvature of the instantaneous undulation trajectory, the curvature was calculated using points of the cycle-averaged trajectory sampled at an interval of one undulation cycle period T.

Analogous to CoT, CoM relates the energy spent per unit distance travelled; therefore, a performance increase is indicated by a smaller CoM. With its current definition, CoM is not suitable to assess straight swimming.

Undulatory swimming shows a characteristic periodicity stemming from vortex shedding during the peak and trough of the body motion. Straight-line swimming at a constant cycle-averaged speed is described as the quasi-steady state at which body forces are in balance over one undulation cycle. To initiate turning, force symmetry is broken by the curving of the body’s central line. Changes in the position of the centre of mass, moment of inertia, and resulting biased cycle-averaged loads lead to a net moment and subsequent angular acceleration. Figure 5 plots the yaw moment of an accelerating and turning swimmer, according to the setup of Figure 1 for Re = 2,000 and constant curvature. Three distinct periods are visible, i.e., a transition period from a static state to a quasi-steady state (t = 0–25 s), a transition period during which the central line curves into an equally distributed curvature of κ = 8.1 rad/m (t = 25–27 s) and, finally, a period of a quasi-steady turning state (t = 27–35 s).

The vorticity contour of each period is shown in Figure 6. Beginning at the quasi-steady state, the reverse Karman vortex street remains horizontal and periodically symmetric. During the transition stage, the body curves and changes the heading direction. Meanwhile, vortices are shed at a non-zero angle. In the third stage, the re-oriented swimmer continues undulation around the curved central line and reaches the quasi-steady turning state. It is found that the time-averaged angular velocity appears to be zero during straight swimming ( t < 25  s). After that, the swimmer maintains a time-averaged angular velocity of 0.2 rad/s. Surprisingly, turning has only a negligible effect on the heading directed velocity.

The input power to fulfil the undulation motion reflects the effort of the swimmer to accelerate against the surrounding fluid. As shown in Figure 7, during the initial acceleration, the longitudinal force peak coincides with an input power peak at approximately t = 7 s. Herein, the longitudinal force and power peak also correspond to the acceleration peak. From there onwards, the swimmer continues to accelerate to the targeted velocity but at a decreasing rate. At the quasi-steady turning state from approximately t = 27 s, the cycle-averaged power converges to a stable value close to zero. The power consumptions during the quasi-steady state straight-line swimming and turning are close to each other, which suggest that the extra effort required for turning, once it is initiated, is insignificant.

Observations of the aforementioned three states are consistent across different curvature magnitudes and investigated Reynolds numbers (Re = 1,000, 1,500 and 2,000). A linear increase in the magnitude of the constant curvature envelope results in a close to linear increase in power consumption. The resulting turning angle decreases with increased curvature yet flattens out at the end, as shown in Figure 8. This results in a decreasing CoM for an increased curvature. For all Reynolds numbers, a similar pattern was observed, where the radius is the dominant variable in the CoM calculation due to its reduction with an increased curvature.

Fish turn both their head and bend their tail in a turning manoeuvre [11]. In the following, the contribution of the head and tail movement towards turning was investigated by comparing three curvature envelopes with the same total curvature. The three selected envelopes are a constant envelope (equal offset across all joints), a linearly increasing envelope denoted as l i n u (predominantly tail curved) and a linearly decreasing envelope denoted as l i n d (predominantly head curved). The kinematic parameters of these cases are shown in Table 2. The Reynolds number in this part is fixed at 2,000.

As shown in Figure 8A, the cycle-averaged drag shows two distinguishable phases during the turning transition period. These are highlighted as phases I and II. Following general observations, negative drag forces correspond to swimming velocity acceleration and positive drag forces correspond to swimming velocity deceleration in the heading direction. Therefore, Phase I may be associated with the initial body curving against the longitudinal moving flow, leading to an increase in drag forces and subsequent deceleration. Likewise, Phase II may be associated with the first full sweep with a curved centre line, increasing the thrust and thus acceleration.

When comparing the constant envelope with the average of both linear envelopes, the curves show close agreement (see Figure 9). It leads to the conclusion that the constant envelope can be segregated into head and tail contributions. According to Figure 9B, the predominantly tail turning envelope achieves a sharper turn with a smaller turning radius. This, together with the close match between a constant envelope and averaged linear envelopes, suggests that head turning has a negative effect on the turning performance.

A comparison of the heading angle and passing fluid velocity angle may provide an explanation. As shown in Figure 10, in the predominantly tail curving envelope, the heading angle leads the fluid angle (i.e., the relative angle of attack of the flow). In contrast, for the predominantly head curving envelope, after an initial peak, the heading angle is behind the fluid angle. For the constant envelope, the angles are closely aligned, with the heading angle slightly leading the velocity angle.

A leading heading angle may positively influence the turning performance, where it provides additional moment for the rotation and reduces the counter rotation moment. Figure 11 shows a schematic representation to highlight the mechanism. A strong tail sweep may provide sufficient moment and energy to turn the swimmer in front of the passing fluid stream by providing additional pressure force to create a moment in the turning direction. Meanwhile, as the tail now sweeps in the opposite direction, the fluid force acting on the head may provide a dampening effect that reduces the counter turning moment. Evidence of this can be seen in the reduced counter turning moment amplitude shown in Figure 12. On the other hand, a curved head may negatively influence the turning performance where it does not create sufficient moment to rotate the swimmer in front of the fluid stream. Additionally, the reactive force during curving of the head leads to a counter turning moment and, as a result, reduces the positive turning effect of the tail curvature. Furthermore, a curved head may provide less resistance during counter turning undulation.

The negative contribution of head turning on the turning performance is also visible in the cycle-averaged power curve in Figure 13.