Thrust generation and propulsive efficiency in dolphin-like swimming propulsion

Given growing interest in emulating dolphin morphology and kinematics to design high-performance underwater vehicles, the current research effort is dedicated to studying the hydrodynamics of dolphin-like oscillatory kinematics in forward propulsion. A computational fluid dynamics method is used. A realistic three-dimentional surface model of a dolphin is made with swimming kinematics reconstructed from video recording. The oscillation of the dolphin is found to enhance the attachment of the boundary layer to the posterior body, which then leads to body drag reduction. The flapping motion of the flukes is found to generate high thrust forces in both the downstroke and the upstroke, during which vortex rings are shed to produce strong thrust jets. The downstroke jets are found to be on average stronger than the upstroke jet, which then leads to net positive lift production. The flexion of the peduncle and flukes is found to be a crucial feature of dolphin-like swimming kinematics. Dolphin-inspired swimming kinematics were created by varying the flexion angle of the peduncle and flukes, which then resulted in significant performance variation. The thrust benefits and propulsive efficiency benefits are associated with a slight decrease and slight increase of the flexion of the peduncle and flukes, respectively.


Introduction
Locomotion in water places high demands on aquatic animals. Natural evolution has produced many species of fish and mammals with high aquatic propulsive performance, and the efforts are growing in engineering to emulate kinematics and morphology of fishes and marine mammals (dolphins, whales, seals, sea lions) towards designing high-performance underwater vehicles. Many such efforts have been successful in producing bio-inspired underwater vehicles that greatly outperform their conventional propellerpowered counterparts [1][2][3][4][5][6][7][8]. Specifically, there has been a growing interest in dolphin-like robots [3][4][5][6][7]. Differing from many fishes that exhibit undulatory and oscillatory kinematics in the lateral plane for forward propulsion, the dolphin exhibits a dorsoventral lift-based oscillatory propulsion kinematics [9][10][11], granting it unique thrust producing abilities and maneuvering habits. Nakashima et al were among the first to design a dolphin inspired robot, with a dolphin-like ellipsoidal body shape, and whose dorsoventral oscillating motion is actuated through two joints, one at mid-body, and the other at the peduncle [7]. Their 2-link design did not only achieve high propulsive efficiency, but also demonstrated impressive maneuverability. Concurrently, Yu et al developed several iterations of dolphin-like multilink propulsive systems and documented their efforts well [4][5][6]. Yu's design divided the dolphin-like body shape into three sections including a rigid anterior body beginning at the dolphin's rostrum and ending roughly at the trailing edge of the dorsal fin, a flexible posterior body from the trailing edge of the dorsal fin to the junction of the peduncle, and the flukes. The flexible oscillatory motion of the posterior body is actuated through a sequence of joints, whose vertical excursion from the center line was governed by a sinusoidal function [4]: h(x, t) = h T (0.21 − 0.66 x LB + 1.1( x LB ) 2 + 0.35( x LB ) 8 ) sin(2π f t), with h being the vertical excursion of the flukes whose maximum is h T , x being the longitudinal coordinate from the nostrum, L B being the body length of the animal, f being the tail beat frequency, and t being time. The resulting midline kinematics and fluke trajectory resemble those of a dolphin. The original prototype had rigid flippers (pectoral fins) and dorsal fin. The rigid flippers were later replaced with flippers with 3-Degree-of-Freedom (3-DOF) control surfaces capable of leadlag, feathering, and up-down motions, with crosssectional shape following the shape of a NACA0018 airfoil, resembling that of a spotted dolphin. The 3-DOF flippers granted the dolphin robot impressive pitch and yaw maneuverability, leading to its ability to leap out of water resembling the 'porpoising' behavior of many species of dolphin [5,6,12]. Additional to posterior body joints, Liu et al employed a system of two push-pull tendons running along most of the dolphin robot's body length to ensure continuum of motion [3].
Design of dolphin-like swimming robots can greatly benefit from biological insight. Though aforementioned dolphin robots have demonstrated impressive controllability and maneuverability, real dolphins' ability to sustain long periods of fast swimming and agility to transition between different speeds have been unmatched. The secret to dolphins' high performance lies in many aspects, including body morphology, skin property, and kinematics [8,11,[13][14][15][16]. Much work has been done in studying the effect of dolphins' body morphology and skin property. Dolphins have a bulky body with a large cross section which is the primary source of drag, and the dolphin's ability to overcome this drag and maintain high swimming speeds has fascinated biologists [13,17]. Studies have attacked the problem from two different perspectives: drag reduction and thrust production. The reduction of pressure drag relies on delaying the separation of boundary layer [18]. This separation can cause the pressure to decrease near the posterior body and thus creating an adverse pressure gradient. In terms of body morphology, drag reduction can be achieved through dolphin's streamlined shape [13], especially in the plane of oscillation [15]. The body reaches maximum thickness near the leading edge of its dorsal fin [13]. The narrowing of the body posterior to the dorsal fin can cause flow transition from laminar to turbulent, which can cause flow separation [19]. But it has been observed both in bio-luminescent environments and in simulation of dolphin-like body that the flow over the dolphin body became turbulent well before reaching the point of maximum thickness, and the turbulent boundary layer stays well attached to the thickest portion of the body and thereafter [13,20]. The early transition of the boundary layer from laminar to turbulent has been conjectured as a result of dolphin's skin texture [13,14,21]. Though measurements have shown the dolphin skin to be incredibly smooth and unlikely to alter the boundary layer [16], the dense packing of dermal papillae with cutaneous ridges that result in skin smoothness can grant dolphins heightened sensitivity to flow disruption and ability to detect boundary layer separation [14], prompting the modification of body kinematics. Drag reduction through body movements was hypothesized to be achieved through the acceleration of the flow near the posterior body region [14]. In Fish and Hui's review study it was found that dolphin's body oscillation is confined to its posterior third [13]. This has largely been realized and followed by roboticists to produce successful results in mimicking the body kinematics of dolphins. However, in such kinematics modeling, the attention has largely been dedicated to the kinematics of the body as a whole, even though the flukes and the peduncle have been seen to significantly vary its flexion in different maneuvers, such as when performing a tail stand [22], which had force production requirements different from those of forward swimming.
The propulsive performance and hydrodynamics associated with dolphin's oscillatory motion have been investigated both experimentally and computationally. It has been commonly held that dolphins propel themselves through the oscillation of the tail flukes [13,23], which has been found to achieve an astonishing propulsive efficiency of 0.75-0.9, surpassing that of most manufactured propellers [24]. In terms of muscle control of the fluke, it was found that the dolphin can control the flexion of the fluke peduncle (junction between the body and fluke) and the cambering of the fluke itself through a continuous bundle of muscle fiber [24]. Through digital particle imaging velocimetry (DPIV), the wake vortices of a dolphin in forward propulsion were visualized by Fish et al finding that the reverse von-Kármán vortex streets in the plane of fluke oscillation produce backward momentum jets that propel the dolphin forward [25].
In Tanaka et al's study of the dolphin kinematics in rapid acceleration, an asymmetry of thrust generation and vortex wake generation was discovered between the up-and downstrokes, wherein downstrokes were found to induce more acceleration than the upstrokes [26]. This result was opposite to the expected asymmetry of the stroke based on muscle mass. The dorsal propulsive musculature of dolphins is larger than the ventral muscles inferring a more powerful upstroke [11].
Three-dimensional vortex wakes produced by dolphins have been visualized computationally, and diverging rows of vortices were found advected from dolphin models performing forward swimming and leap [27][28][29]. In Han et al's computational study of dolphin in forward swimming whose realistic oscillatory kinematics were reconstructed from high-speed video recordings, four vortex rings were identified in the wake of a dolphin in forward propulsion, produced during the up-and downstrokes and the transitions between them [28]. Using similar techniques, Wang et al analyzed the hydrodynamics of a dolphin swimming at three different speeds, finding a similar vortex structure, as well as significant variation in the effective angle of attack of the fluke [29], indicating the important role that the flukes play in dolphin's forward propulsion. In multiple computational studies, a thrust-producing distribution of pressure was found in the posterior region of the body [26,28,29]. This corresponds well with the observation by Fish that the flow is observed to accelerate over the posterior body, leading to flow attachment and drag reduction [14].
Given the broad interest in dolphin-like swimming, and the gap in knowledge in flow modification by body kinematics, specifically the effect of varying fluke flexion kinematics, and in order to better understand the mechanisms of thrust generation of dolphin-like propulsion, this work offers a detailed analysis of the hydrodynamics of a dolphin in forward propulsion, and produces dolphin-inspired swimming kinematics through the scaling of the bending of the flukes. The wake structure and flow around the dolphin are visualized with a direct-numerical simulation (DNS) computational fluid dynamics (CFD) solver. The varying effectiveness of different fluke flexion in forward propulsion is examined. The subsequent contents are organized as such. The process of kinematic data collection through high-speed video recording the dolphin swimming, and the reconstruction of a computational dolphin model and oscillatory kinematics will be described in section 2.1. In the same section, a brief discussion of the kinematic characteristics of the dolphin will be given, to highlight key features that set it apart from fish swimming. The CFD solver and the computational domain setup will be described in section 2.2. A case is set up removing the oscillatory kinematics; the dolphin model is towed at a constant speed that is the same as the oscillating dolphin in a neutral posture. The goal for this study comparing the dolphin with and without oscillation kinematics is two folds. 1. By isolating the kinematics, the role it plays in dolphin propulsion and key benefits can be discovered. 2. By removing the kinematics, the role of dolphin morphology on its performance can be isolated and drag reduction mechanisms related to the shape of the dolphin can be studied. Section 3.1 will examine the simulation results from the dolphin with and without oscillatory kinematics. Section 3.2 will discuss the results of parametric study varying the flexion magnitude of the peduncle and fluke.

Swimming kinematics reconstruction
High-speed, high-resolution video recordings were collected on the swimming of real dolphins to reconstruct the swimming kinematics. Seven bottlenose dolphins (Tursiops truncatus; two males and five females) were video recorded at 60 frames s −1 with a Canon camera as they swam freely, or swam rapidly under trainer command and performed tail stands in a large tank at the National Aquarium in Baltimore, MD. Two Canon EOS 5D Mark III cameras were positioned on tripods in front of a large viewing window. The cameras were set up to collect video from the lateral view normal to the dolphin, as shown in figures 1(a1) and (b1). Additionally, an oblique-view camera was set up to record from slightly behind the dolphin as it swam or was performing a tail stand. A three-dimensional calibration device was placed in the field of view of both cameras before and after animal trials. The calibration device consisted of six aluminum rods (9.7 mm diameter, 0.61 m long) that were screwed into each face of a 50 mm nylon cube, providing an orthogonal arrangement of the rods. Two 38 mm diameter nylon balls were attached to each rod with one ball at the end of the rod and another positioned equidistantly on the rod at 0.3 m between the centers of the terminal ball and the central cube. The research on the dolphins was approved by the West Chester University Institutional Animal Care and Use Committee under 201 701 and F2011. A 3D surface model of the dolphin was then created in Autodesk Maya. A Catmull-Clark subdivision method used in the meshing of the 3D surface ensured the smoothness of the surface of the dolphin model. A simplified surface mesh that subdivides the body is shown in figures 1(a2) and (b2). This basic, rectangular surface mesh was then subdivided and triangulated to achieve the surface smoothness shown in figure 1(c). The high level of surface smoothness ensured compatibility of immersed boundary method (IBM) used in the DNS [30]. The surface model was rigged up with an anatomically realistic virtual skeleton, as shown in figure 1(c), to control the model dolphin's posture and corresponding skin deformation in a realistic manner. A virtual scene was set up mimicking the setup of the two video cameras in the real scene, with the rigged surface model placed in the middle. The use of two different camera angles allowed for calibration of the virtual scene to the real scene. The dolphin swimming kinematics was thus reconstructed by deforming the surface model via the virtual skeleton to match the image of the dolphin frame by frame. This technique has been used to accurately produce biologically realistic motions of jackfish and manta ray swimming [31,32] allowing for an accurate reconstruction of the body kinematics from the video recording. The sequence of frames dis-   [10]. The head and fluke oscillation amplitudes of the dolphin in this study fall within range of experimental measurements [33], and are larger than most fish [10]. The headfluke oscillation amplitude ratio comes to 0.21, falling well within fish's range of head-tail oscillation amplitude ratio [10].
One key characteristic of the dolphin forward propulsion kinematics that sets it apart from fishes is the asymmetry of its tail stroke. Discounting small variations between the two strokes (left-to-right and right-to-left) in a typical cycle, fish's midline kinematics is largely symmetric [9,10], due to a similar left-right symmetry of fish's body morphology and propulsive musculature. However, for the dolphin, the upper envelope of the mid-lines shows more pronounced curvatures than the lower envelope, as shown in figure 2(a). Higher curvature on the upper envelope of the mid-lines indicate the pronounced arching back of the dolphin during its swimming, especially during the transition from upstroke to downstroke. This can be attributed to the asymmetry in the body morphology, as the dolphin's dorsal surface is marked by the existence of a dorsal fin of a significant size compared to its body size, whereas the ventral side contains a pair of pectoral flippers, which are more anterior along the body than the dorsal fin. The body itself also shows a thoracic region that is bulkier than the abdominal region, contributing to the asymmetry of oscillation.
Pronounced bending of the flukes can also be observed in the center-line kinematics schematic, as exemplified by the sharp turn on the upper envelope, at the point corresponding to the peduncle near the flukes in figure 2(a). From the instantaneous midline curves over one propulsive cycle, the cycle average curvature at various points along the body was also estimated, and plotted in figure 2(b), with the cycle standard deviation also shown. The calculation of the curvature was taken as the inverse of bending radius. Most of the body lays relatively flat, except for a small bump in the middle of the body, corresponding to a location roughly just behind the dorsal fin. However, a notable peak can be seen at the peduncle near the flukes. Directly following it is a smaller yet still significant peak indicating pronounced bending in the middle of the fluke. Fish et al found that the bending of the peduncle at the anterior insertion of the flukes and cambering of the flukes proper are connected through a continuous bundle of muscle fiber and associated tendons [14,24]. With the maximum flexion corresponding to the position of the ball vertebra near the fluke base, it was determined that the cambering of the fluke is a passive maneuver, linked tightly to the flexion of the peduncle [14,23,24]. The two peaks in the curvature plot occurring at the anterior position of the insertion of the flukes on the is multiplied by a flexion ratio (FR) to set up cases whose kinematics varied from the original, baseline case. Cases in this parametric study have instantaneous flexion angles that can be described as θ f = FR × θ f0 . FR for the original kinematics of the dolphin was 1.00. We varied FR by one-third increments, increasing the FR to 1.33 times and 1.67 times the original, and decreasing to 0.67 and 0.33 times the original. Higher FR leads to enhanced instantaneous θ f . The instantaneous θ f of the cases in the parametric study are shown in figure 2(d). As a result of changing FR, the fluke tip trajectory over one oscillating cycle is also slightly altered, and the instantaneous vertical excursion of the fluke tip over one oscillating cycle is shown in figure 2(e) for different parametric cases. Decreasing the FR had an overall increasing effect on the maximum vertical excursion. During the transition from upstroke to downstroke (t/T ≈ 0.05 to t/T ≈ 0.15), the fluke tip reaches peak vertical excursion y/L B . As shown in figure 2(e), the peak fluke tip y/L B is the highest for FR = 0.33 and the lowest for FR = 1.67. This is caused by the fluke tips extending upward in y-direction at the top of an upstroke (figure 2(f)) with reduced flexion (low FR), contrasting the case of enhanced flexion (high FR), in which the flukes bent downwards at the top of an upstroke and hence reducing the vertical excursion of the fluke tips. Similarly, at the bottom of the downstroke, reduced flexion caused the fluke tips to extend further downward, and hence arriving at lower minimum y/L B values, in the low-FR case. Increasing FR from a value of 1 also cause the fluke to extend further downward at the bottom of the downstroke; as the peduncle spent more time to transitioning from downstroke to upstroke, the flukes with enhanced flexion (such as FR = 1.67) also had time to transition from bending upward to bending downward, and thus helping the fluke tips to reach a minimum y/L B that was lower than the minimum   (1), in which f is the oscillation frequency, A is the fluketip amplitude (figure 2(a)), W ∞ is the incoming flow velocity, L B is the total length of the dolphin, and ν is the kinematic viscosity. The key kinematic parameters and their values are listed in table 1. The selection of Re ensures reasonable computational load, even though the real dolphin swam at a higher Re of 2 × 10 6 . The appropriate Strouhal number was subsequently selected to ensure that the sum of streamwise forces was approximately zero. Though a small difference in fluke tip amplitude is seen among the parametric cases, as shown in figure 2(e), the Strouhal number is kept constant for all parametric cases (1)

The direct numerical solver
The fluid dynamics of dolphin swimming was solved using a DNS developed in-house. An IBM allowed the solver to resolve the surface boundaries of the dolphin model and compute the fluid behavior around the body. With kinematics described in the previous subsection, the oscillating dolphin model was placed in a computational domain of size (15L B ) 3 , set up in a Euclidean space with an x-y-z coordinate system. The dolphin model was oriented along the z axis, with its rostrum pointing in the negative z direction and the dorsal fin pointing in the positive y direction. The dolphin was thus oscillating in place, while an inflow boundary condition with a constant velocity W ∞ was applied at the negative-z boundary, an outflow boundary condition was applied to the positive-z boundary, and no-gradient conditions were applied to the flow velocity at the positive-and negative-x and y boundaries. The pressure at all the boundaries were governed by the Neumann boundary condition. A Cartesian mesh was used to discretize the domain spatially. A schematic of the spatial discretization is shown in figure 3(a). The dolphin is situated 10 L B 's from the inlet, in the center of the xy-plane. An abundance of space is given to the flow upstream of the dolphin's rostrum to allow a uniform flow to fully develop. A total length of 4L B left behind the dolphin's flukes to capture its full wake. A dense region with uniform grid spacing ∆x = ∆y = ∆z = 0.0118L B was set up around the dolphin body. Stretching grids surrounded the dense region, with the stretching ratio in the wake region behind the dolphin set to 0.65, so that the grid spacing at the outlet (∆z = 0.05) was small enough to resolve vortices of comparable dimensions. Nearfield of the body, two nested tree-topological local mesh refinement (TLMR) blocks [34] were used to further refine the dense-region grid, to ensure that the boundary of the dolphin model and the pressure and viscous forces acting on it could be fully resolved. A zoomed-in view of the spatial grid in the TLMR region is shown in figure 3(a). The dolphin body is fully contained within the inner-most TLMR block (pink), which is then fully bounded by the intermediate TLMR block (purple). The intermediate TLMR block extends 0.5 L B into the wake region downstream of the fluke tip. To ensure visibility, only every other grid point is shown in the TLMR-region schematics. The grid size in the TLMR region is uniformly 0.003 L B . The total grid count of the entire computational domain including the base mesh and the TLMR blocks is 16.0 million.
Non-dimensional coefficients of forces and power were calculated using equation (2), where C T , C L and C PW are the coefficients of thrust, lift and power, respectively, F Z and F Y are the z-and y-directional pressure forces acting on the dolphin body, with thrust and lift pointing in the negative-z and positivey directions respectively, P W is the hydrodynamic power, ρ is the fluid density, S is the surface area of the dolphin body and W ∞ is the inlet free-stream velocity. The hydrodynamic power, P W is defined as the negative dot product of the body element's velocity and force vectors, integrated over the surface of the body, the flukes or the entire model. A positive C PW then indicates the consumption of power, as the body mesh elements move in opposition to hydrodynamic forces acting on them, and negative C PW indicates forces assisting the motion of the body or the flukes. The feasibility of the current grid density in resolving forces acting on the dolphin body was studied through a grid-independent study. A coarse mesh with minimum grid spacing of ∆z = 0.006L B and fine grid with minimum grid spacing of ∆z = 0.0019L B were set up for benchmarking with the nominal mesh just described. Based on the surface forces computed using the three different grid setups, the thrust coefficient (C T ) was calculated using equation (2). The results from all three meshes are plotted in figure 3(b) for comparison. Qualitatively all three meshes agree well with each other, as the basic shape is retained in the coarsest mesh. Quantitatively, between the nominal and fine mesh, difference in the peak value was found to be 1.5%, while the difference in grid density was 50%. It is thus concluded that the nominal mesh was sufficient for our computation, and therefore the subsequent analyses were conducted with the nominal mesh

Hydrodynamics of baseline kinematics
The pressure forces acting on the body of a steady swimming dolphin and the same dolphin model in a neutral position being towed at the same speed are summarized in table 2. For both cases, simulations are continued for long enough for asymptotic steadystate values to emerge. The coefficient of lift, C L , is calculated in both cases using equation (2). The steadystate C T and C L are shown for the towing case, and the cycle-average coefficients of thrust and lift,C T andC L after steady state has been reached are shown for the forward swimming dolphin. Without oscillation, in the towing case, both the body and the pair of flukes produce net drag. This drag is primarily constituted of the form drag due to the pressure gradient along the streamwise direction. In the forward swimming case, the oscillation kinematics is found to greatly reduce the net drag experienced by the body, and the fluke becomes fully thrust producing. It is also found that the dolphin's morphology is moderately lift-producing, as shown by the net positive lift on the body of the towed dolphin. The body and tail flukes all contribute positively to lift generation when being towed, making up 63.6% and 36.4% of the total form lift. With oscillating kinematics, the dolphin's lift production multiplies, by 16 folds for the flukes and 6 folds for the body. The flukes generate a significant portion (61%) of the total lift on the forward swimming dolphin, demonstrating the flukes as powerful propulsive and control surface.
Interesting points to note so far in summary include, in the forward swimming case, (1) the reduction of body drag, (2) the production of forward thrust on the flukes and (3) lift generation of the fluke in forward swimming, and, in the towed case, (4) the lift generation of a static dolphin body. In the following analyses, we focus on the propulsive performance (1)-(2) with comparisons with the towing case as necessary, while lift generation (3)-(4) takes a minor role.
First, the flow around the towed dolphin is examined. As shown in figure 4(a), the forces have settled to a stable asymptotic value, with variation in C T between t * = 2 and t * = 3 less than 3 percents. Figures 4(b1) and (b2) visualizes the vortex structure around the dolphin using Q-criterion isosurfaces. It can be seen from figure 4 that vortices start to form near the head of the dolphin, agreeing well with early development of turbulence observed by Riedeberger and Rist [20]. When the flow encounters the pectoral fins, it rolls around its leading edge forming a shear layer that covers most of the pectoral fin's upper and lower surfaces, and is shed at the trailing edges. Near the pectoral fin roots, vortex streaks can be seen to form in figure 4. Similar vortex dynamics can be found on the dorsal fin. Vortices on the dorsal fin separated and were advected downstream, and minimal shear layer and vortices were seen in the narrowing body posterior to the dorsal fin. Figure 5 shows the wake structure and vortex dynamics of the oscillating dolphin through a combination of Q-criterion isosurfaces and vorticity contours. Similar to the towed dolphin, the oscillating dolphin's dorsal fin produces vortices passively due to the incoming free-stream flow and body shear   layer. Thus the dorsal-fin vortex (labeled DFV in figure 5) has a horizontal cross section ( figure 5(a2)) consisting of a drag-type pair of counter rotating vortices. Contrasting to the lack of vortex formation on the posterior body of the towed dolphin as observed in figure 4, around the oscillating dolphin, rotational flow was found to form in the posterior body, as depicted in figures 5(a1) and (b1), and were well attached on the body until reaching the flukes. At the flukes, the flow is greatly altered, as demonstrated by the isosurfaces around the flukes having markedly higher overall Q-criterion values than those around the body, revealing the flukes as playing a major role in dolphin's oscillatory swimming. As the flow continues downstream past the flukes, vortex rings are formed as shown in figures 5(a1) and (b1). As shown in the lateral view in figure 5(b1), a bifurcation of the wake structure is observed in the plane of oscillation, similar to fish swimming [35]. In the wake, two main vortex rings R D and R U are found to form, associated with the downstroke and the upstroke respectively. During the transition from upstroke to downstroke, another vortex ring, R T , is formed. Figures 5(b2) and (b3) shows the vortices in the cross sections of R T and R U , shown on slice cuts made at locations shown in figure 5(b1). The cross section of each vortex ring consists of two counter-rotating vortices that create a backward  facing jet in the middle, indicating the forming process of R T and R U as important for thrust production. By comparing the vorticity magnitude of R D and R U and the flow velocity in the middle of either vortex ring represented by the vectors, it can be seen that the downstroke produces a stronger vortex ring and associated jet, indicating it as the more effective at producing thrust. The jet production in the wake of the oscillating dolphin due to the shedding of R U and R D can also be seen in the time-averaged flow velocity, shown in figure 6(a1). The red and yellow region indicating high-momentum flow in the back of the dolphin is the jet produced by the oscillating dolphin. The far-wake jet in figure 6 can be seen diverging in three different directions corresponding to the production of R T , R U and R D , with the lower-most jet, corresponding to R D produced during downstroke, being the strongest. This provides supporting evidence for the slight downstroke bias of dolphin's propulsion. The blue region above the dolphin behind the dorsal fin shown in figure 6(a) indicates the flow being slowed down below the free-stream velocity, a sign of drag production. The time-average flow velocity contour for the non-oscillating, towed dolphin, shown in figure 6(a2) provides a contrast to the dolphin in oscillation, to highlight its benefits in drag reduction. The static dolphin exhibits a much more prominent and extensive slow-down of the flow in the region above the posterior body compared to the oscillating dolphin. The prevention of the slow-down and the minimization of its area of influence in the posterior body by the propelling dolphin shows the effect oscillation has to accelerate the flow around the posterior body, which was hypothesized to lead to boundary-layer attachment and drag reduction by Fish [14]. The vector field in figure 6 shows that associated with the lower-most jet, the flow is not only accelerated backward, but also downward. The downward-pointing velocity vectors are overall stronger than the upward-pointing ones, indicating a time-averaged net positive lift force generated. Figure 6(b) displaying the y-velocity contour immediately behind the towed dolphin without oscillation shows a net downward velocity in the wake. This indicates that dolphin body shape's effect in directing the flow downward, creating a down-wash like an asymmetric airfoil, leading to lift generation as seen in table 2. The downward velocity vectors shown in figure 6(a1) shows oscillation's enhancing the shapeinduced down-wash, leading to its increased lift force seen in table 2.
To examine the instantaneous forces acting on the oscillating dolphin and its power consumption, the time histories of C T , C PW and C L over a representative propulsive cycle are shown in figure 7. The coefficient of power, C PW is calculated using equation (2). Due to the high Reynolds number, short-period oscillations are observed in the force time history. However, an underlying long-period trend can still be observed. For both the body and the pair of flukes, C T and C PW follow similar trends. At the beginning of the cycle, both thrust production and power consumption are low, corresponding to the flukes' stroke transition as seen in figure 2(e). During the downstroke, the thrust force rapidly grows, and correspondingly, power consumption also grows. At the bottom of the downstroke, during stroke reversal, the thrust and power both drop to low values again. During the upstroke, the flukes maintain high thrust similar to those during downstroke, while the power consumption is lower than that of the downstroke. Corresponding to this, in figure 2(d), the peduncle is seen maintaining higher flexion angles during the upstroke than during the downstroke. As a result, the vertical force experienced by the flukes is reduced in the upstroke, as shown in figure 7(c). The fluke flexion is playing an important role in reducing power consumption while maintaining high thrust production, by reducing crossstream force components. The body experiences net drag over an entire cycle, as shown in table 2, and the instantaneous value oscillates rapidly between drag and thrust.
For lift generation, the body is found to oscillate rapidly between positive and negative, while the flukes undergo a more gradual process: during the downstroke, the lift gradually rises and falls as stroke reversal approaches, all the while staying mostly positive, and after the stroke reversal, the lift continues to fall below zero, and rises again as the next stroke reversal is approached.  figure 7(a), the body C T shows a local maximum. Correspondingly, the surface pressure distribution shows that in the posterior portion of the body, a thrust-beneficial pressure gradient is seen whereby high pressure is concentrated on the dorsal surface, pushing the dolphin, and suction force is concentrated on the ventral side of the posterior body, pulling the dolphin forward. Pressure gradient now subsides on the flukes, as it cambers and transitions from down-to upstroke. At t/T = 0.75, during upstroke, the flukes exhibit a strong pressure gradient, producing high thrust, and the body experiences moderate pressure gradient, and with little body deformation, the pressure gradient acts to create force in the y-direction, generating a net lift. As the dolphin finishes an upstroke, and transitions to downstroke, at t/T = 1.00, the pressure gradient on the flukes once again subsides, and the high-pressure region on the body moves back to the dorsal side, and the ventral side of the body has predominantly negative pressure.
The pressure distribution in figure 8 shows that though the posterior body occasionally aids with thrust production due to body deformation, the flexion of the peduncle allows the flukes to better utilize its pressure gradient to produce thrust. Figure 7 also shows that though the body experiences large amount of force, it mostly results in oscillation in vertical forces, whereas the flukes dominate over the body in producing thrust force. Next, the fluid dynamics of the flukes' thrust production is examined through an analysis of vortex forming and shedding over an entire propulsive cycle, as shown sequentially in figure 9, with the key vortices labeled. Parts (a1)-(a3) of figure 9 depict the flukes during the downstroke, and parts (b1)-(b3) depict the flukes during upstroke. The vortex forming and shedding during the downstroke (figures 9(a1)-(a3)) and the upstroke (figures 9(b1)-(b3)) are shown in the dorsal perspective and ventral perspective views, respectively, to ensure visibility of active vortex structures.
Leading-edge vortices (LEVs) are found to form and shed periodically during a flapping cycle. Its initiation occurs at the beginning of a stroke: as seen in figures 9(a1) and (b1), after a stroke reversal, coherent LEVs are fully formed. During flapping, periodic forming and shedding of multiple LEVs is observed during a stroke (half cycle). The coexistence of multiple LEVs is depicted in figures 9(a2) and (b2). Multiple Vortex forming and shedding is hypothesized as the mechanism behind the oscillation of thrust from the flukes observed in figure 7. The LEVs shed during an active flap are found to form the left and right bounds of the vortex rings R U and R D , as shown in figures 9(a2)-(a3) and (b2)-(b3). In the stroke reversal stage (figures 9(a3) and (b3)), the LEV from both the left and right leading edges of the fluke connect in the middle, and the shedding of this joint vortex encloses the vortex ring, and coherent R D and R U are formed. The process of R T formation during the stroke transition stage can be seen from figures 9(a1) and (b3). At the top of the upstroke, vortex tubes are shed from the shear layer on the flukes, and during the stroke reversal, the chord-wise flexion and the pitching of the flukes around the peduncle causes a strong vortex to form on the trailing edge, which then encloses the vortex loop to form R T .
The process and effect of multi-LEV forming and shedding can be demonstrated by an analysis of the fluid vorticity and pressure on a slice cut, shown in figure 10. Due to the good symmetry between the left and right flukes in their vortex formation and wake structure shown in figure 9, the analysis isolates and focuses on the left fluke, and the x-vorticity and pressure contour of the flow around it are shown at discrete sequential time instances in figure 10. At t/T = 0.15, peak thrust is reached in figure 7. Correspondingly, a coherent LEV is formed and fully attached to the left fluke in figure 10(a1). The highspeed swirling of the fluid associated with the LEV causes the fluid pressure to drop, and hence, a suction region is seen attached to the leading edge of the fluke depicted in figure 10(b1), leading to high thrust production. At the first dip in fluke thrust (t/T = 0.21), the first LEV core is seen to begin to detach ( figure 10(a2)). Correspondingly the core of the LEVinduced suction region is seen to separate from the leading edge of the fluke in figure 10(b2), though still exerting a suction force on the forward-facing upper surface of the fluke. Hence, the flukes only experience a small dip from peak thrust, but still maintain relative high thrust. At a later instant, t/T = 0.27, the first LEV has been fully separated, and advected further downstream, while staying relatively close to the fluke's upper surface, and a second LEV has started to form, as shown in figure 10 The combined effect of the second LEV and the detached first LEV help extend the suction region across the upper surface of the fluke, and recover a portion of the peak fluke thrust, as shown by the second small peak in C T in figure 7(a). The detached LEV creating a low-pressure region to still aid thrust production is in good agreement with Eldredge and Jones finding of flapping plates [36]. In all the instants depicted in figure 10, the fluke cross-section is shown at relatively large angles of incidence relative to the incoming flow, and hence the LEV production of suction region on the upper surface of the flukes creates desirable propulsive force, in agreement with Wang et al [29].

Effect of varying FR
By varying the flexion kinematics, according to the description provided in section 2.1, the performance changes that occur due to deviation from original kinematics can be isolated. Besides thrust force production, the stroke/cycle-average propulsive efficiency (η), calculated according to equation (3), is also compared among the cases. The propulsive efficiency is relevant when the fluke is thrust producing (C T > 0), when on average no net thrust is produced in a time period, the fluke also has no propulsive efficiency (η = 0). The cycle-averaged thrust production and propulsive efficiency, summarized in figure 11 shows that both high thrust and high propulsive efficiency are achieved near FR = 1. The peak cycle-average thrust production, maxC T is found when the flexion kinematics is scaled down to 67% of the original, 11.2% above the value reached with FR = 1, showing that swimming with a slightly stiffer peduncle is more conducive to the production of thrust force, whereas increasing the FR, i.e. swimming with a more flexible peduncle can cause a rapid loss of thrust force. The propulsive efficiency trend, on the other hand, favors a slightly enhanced (<+33%) FR, i.e. dolphin swimming with more flexion at the peduncle. With minimal difference (<1%) inη between FR = 1.00 and FR = 1.33, the peak is expected to be found with FR ∈ (1.00, 1.33). However, further increase of FR leads to rapid decline of propulsive efficiency. The upstroke is especially sensitive to FR change. With FR = 1.67, the flukes lose thrust completely and start experiencing drag during the upstroke, as shown in figure 11(a), and hence the upstroke-average propulsive efficiency is 0, as shown in figure 11 During the downstroke, the performance variation due to changes in kinematics is comparatively small. The peak value of downstroke-average thrust is found in the dolphin with the original kinematics, instead of at a reduced flexion of FR = 0.67. Though scaling up of flexion kinematics leads to diminishing downstroke-average thrust production and cycleaverage efficiency, the downstroke-average propulsive efficiency continues to rise with increasing FR. Next, the flow physics leading to variations in propulsive performance in the upstroke, the downstroke, and the whole cycle will be discussed.
The cases isolated for comparison here are 1. the case with the highest cycle-average thrust production (FR = 0.67), 2. the baseline (FR = 1.0) and 3. (FR = 1.33) the case with the low cycle-average thrust but a high cycle-average efficiency. By comparing case FR = 0.67 with the other two, hydrodynamics of high thrust production can be discovered, and comparing case 3 with the other two, hydrodynamics behind efficiency variations can be discussed.
Qualitatively, the time histories of the cases follow similar trends, with the thrust and power consumption peaking in the middle of a stroke and subsiding near reversals of strokes, as shown in figure 12. The differences in performance occur both near the peaks during a stroke (t/T = 0.1 to 0.35 and t/T = 0.5 to 0.8 in figure 12) and in the reversal stage (t/T = 0.4 to 0.5, and t/T = 0.95 to 1.05 in figure 12), when the thrust is the lowest. The gap in instantaneous C T is larger during the upstroke than during the downstroke, agreeing with the steeper slope in the upstroke curve than in the downstroke curve, to the right of FR = 0.67 in figure 11(a). The gaps in instantaneous C PW are similarly larger in the upstroke than in the downstroke. The large differences in instantaneous C T and C PW during the upstroke compound to result in large variations of upstroke propulsive efficiency, as shown in figure 11. The instantaneous C L curves show that the vertical force oscillates very little above and below C L = 0 for FR = 1.33. With decreasing FR, peak |C L | increases, indicating higher transverse forces experienced. High non-axial force can lead to excessive power consumption, and therefore possible reduction of propulsive efficiency.
To examine the flow behaviors that resulted in performance shown in figure 12(a) in a general manner, the near-fluke flow is visualized by Q-criterion isosurfaces shown in figure 13. The time instants depicted here correspond to moments of large gap in instantaneous C T between cases with different FRs. At t/T = 0.13 shown in figures 13(a1)-(a3), the dolphin's tail has finished the transition from upstroke, and moving downwards. R U formed during the upstroke is still visible for FR = 0.67 ( figure 13(a1)), but not quite for FR = 1.00 (figure 13(a2)) and FR = 1.33 (figure 13(a3)), indicating a stronger upstroke jet formed in FR = 0.67 resulting in the higher instantaneous and upstroke-averaged C T values observed in figures 11(a) and 12(a). The downstroke LEV (LEV D ) formation at this instant is also more thorough, covering the entire leading edges of the flukes in FR = 0.67 than in FR = 1.00 and FR = 1.33, likely due to differences in angles of attack relative to local flow. Prominent downstroke trailingedge vortices (TEV D ) are also formed at this instant for all three cases. The combination of LEV D and TEV D can result in the formation of coherent vortex ring generating a thrust jet.
During the upstroke, at an instant when a large difference in instantaneous C T value is observed among the three cases (figures 13(b1)-(b3)), the primary upstroke LEV (LEV U ) begins to separate from the leading edge and starts rolling over the ventral surfaces of the flukes, in FR = 0.67 and FR = 1.00, while for FR = 1.33, the LEV U does not separate to two parts. Both the primary and secondary LEV U 's are the most coherent for the thrust-optimal case of FR = 0.67, and the connection between them is stronger for FR = 0.67 than for FR = 1.00; both features are responsible for higher thrust. The flapping of the flukes at this reduced FR of 67% of the original amplitudes led to stronger formations of LEV U groups, which when shed, formed stronger R U as shown in figure 13(a1).
A slice cut is made at the quarter-span location as shown in figure 14, the same as in figure 10 to examine the interior structures of the leading-edge and trailing-edge vortices whose isosurfaces are shown in figure 13. By overlapping vectors indicating local flow velocity on the ω x , the cores of the LEV D and TEV D of FR = 0.67 can be visualized to cause more swirling than for FR = 1.00, and a stronger and more directed thrust jet is created in the narrow channel between the vortex cores for FR = 0.67 than for FR = 1.00. The cause of this at instant t/T = 0.13 can be attributed to the combining of the downstroke trailing-edge vortex and the shed upstroke leadingedge vortex. The down-stream advection of LEV U cores can be seen in all three cases, whose relative strengths differ in a manner consistent with relative thrust production. For all three cases, the shed LEV U is adjacent to the newly formed TEV D . Even though connection between LEV U and TEV D at the edges can be seen in figure 14(a2), for both FR = 1.00 and FR = 1.33 the cores of these vortices are apart from each other, with the LEV U core for FR = 1.33 comparatively much weaker, associated with lower thrust production during the upstroke. However, the reduced flexion (FR = 0.67) acquires a much stronger and more coherent combined TEV D core, due to the shallower angle of incidence of the flukes relative to the incoming flow resulting in the overlapping and absorption of LEV U by TEV D . The result of this combination of TEV D and LEV U is a stronger trail-edge vortex that leads to stronger LEV-TEV jet and thus higher instantaneous thrust.
During the upstroke, the comparison of figures 14(b1)-(b3) highlights the advantage of a shallower angle of incidence during flapping in producing thrust force. While the fluke's angle of incidence varies little among the cases at t/T = 0.13, the difference is more observable at t/T = 0.69. For all three cases, the overall vortex structures are similar, with LEV U formed on the leading-edge and advected over the ventral surface, and TEV U formed on the dorsal surface, advected back and shed at the trailing edge. The vector field indicates stronger LEV U and TEV U circulations for the case of FR = 0.67, and this is confirmed by a circulation calculation shown in figures 15(a3) and (b3). The vector field also indicates stronger jet formation in the case with FR = 0.67, leading to high forward propulsive force.
The cross sections of the LEV and TEV cores are depicted in figures 15(a1), (a2), (b1) and (b2) for cases with FR = 0.67 and FR = 1.33. The instant chosen for circulation calculation is one during the upstroke with C T gaps between these cases, the same as in figures 13(b) and 14(b). The fluke is seen exhibiting different angles of incidence, and as a result, the vortex coherence and strength differ. With FR = 0.67, though the LEV core is seen relatively far from the fluke surface, but the vortex formation is coherent. Whereas, for FR = 1.33, the vortex core is situated in close proximity to the fluke, but the vortex is incoherent; towards the fluke root, the leading-edge vortex core has diminished. The relative strengths of the LEV cores can be compared through circulation calculation. As shown in figure 15(a3), the LEV strength on average grows with decreasing FR, especially towards the middle of the leading edge (quarter-span location). Towards the tip, the LEV circulation diminishes similarly across cases with different FR. For the TEV, a similar trend of diminishing coherence with increasing FR can be observed by comparing figures 15(b1) and (b2). Similar to LEV, the TEV circulation trend for each case shows a trend of increasing magnitude with increasing x/L B , moving towards the mid-span location of the fluke trailing edge. The overall TEV strengths of the FR = 0.67 and FR = 1.00 cases are similarly high, enhanced by their interactions with previously shed vortices. The enhanced circulations of LEV and TEV combined in cases with reduced FR lead to stronger jet production seen in figure 14.
The enhanced LEV circulation in the case with FR = 0.67 at t/T = 0.69 also leads to increased  leading-edge ventral-surface suction. As shown in figure 16, the (blue) suction region a little upstream of the fluke becomes more prominent as FR is decreased, at the same instant as depicted in figure 14(b), due to the high-speed swirling of the fluid in vortical regions with high circulation. In FR = 0.67, this suction region is strong and covers the entire chord of the fluke. On the dorsal surface, a little downstream of the fluke, high-pressure region is found as a result of the large surface of the fluke pushing into the flow In FR = 1.00, the high-strength suction region is reduced to the leading edge, though relatively high suction intensity can still be seen covering the ventral side of the chord. In FR = 1.33, the suction strength at the leading edge becomes much lower than the other two cases, and the rest of the chord suffers a similar reduction of suction. The combination of highpressure region and suction region created beneficial pressure gradients that led to the strong overall force production. Though the propulsive force benefited from this, given the high C T at the depicted instant for FR = 0.67 as shown in figure 12(a), it also means that the upward-flapping fluke encounters stronger resistance, leading to higher power consumption as shown in the upstroke region of figure 12(b), and reduced propulsive efficiency. In the case with FR = 1.33, the  upper surface of the fluke is not dominated by any heightened pressure, and therefore the fluke does not need to overcome excessive resistance to move, leading to a lower power consumption and enhanced upstroke-average propulsive efficiency.
As flukes experience a force from the flow due to a pressure gradient, the flukes also exert a reaction force on the flow, and hence accelerating it. Whereas a stream-wise acceleration of the flow can correspond to forward thrust (C T ) production, a vertical acceleration of the flow corresponds to C L production. Across the chord length of the flukes, the flow accelerated steadily in the stream-wise direction in cases with FR = 0.67, 1.00 and 1.33, as shown in figure 17(b1), indicating steady thrust force production. At the roots of the flukes, in varying FR cases, the normalized net stream-wise flow speed, (w − w ∞ )/w ∞ , exhibited similar magnitudes, and the difference among them grew larger toward the tips. As the FR was decreased, the stream-wise speed of the flow experienced a larger acceleration due across the flukes' chord. For FR = 0.67, 1.00 and 1.33, the increase of (w − w ∞ )/w ∞ from the flukes' root to the tips was 0.051, 0.042 and 0.026, respectively. As the acceleration across the flukes' chord increased with decreasing FR, so did the flukes' C T . Along the vertical direction, the trend of flow speed over different chord-wise locations differs from that of the stream-wise flow speed. As the flukes flapped upward, the vertical flow speed also increased from the root to the tips of the flukes. Unlike the stream-wise flow speed, the vertical flow speed then quickly decreased aft the fluke tips, and crossed zero, as the flow in this region originated from the downstroke. The variation of overall amplitude of the vertical flow speed across cases with different FR was similar to that of the amplitude of the streamwise flow speed: enhanced flexion of a higher FR led to lower amplitudes of vertical flow speed. The lowered vertical-flow speed meant a lower proportion of hydrodynamic power that is diverted from thrust production, and hence a recovery of η despite lowered C T .

Conclusion
In this study, the propulsive performance and hydrodynamics of dolphin-inspired swimming kinematics was studied. The overall shape of the dolphin body also acts like an asymmetric airfoil that deflects incoming flow downwards to create a moderate downwash, leading to lift production without oscillation. Flow past a static dolphin produced shear layers on the thick anterior portion of the body, and the tapering posterior region saw minimal vortex formation compared to the oscillating dolphin. Dolphin-like oscillatory kinematics reconstructed from video recording of real dolphin was found to pre-accelerate the flow around the posterior portion of the body, resulting in body drag reduction and production of thrust by the posterior portion of the body.
Dolphin swimming is found to produce a multiring vortex structure, with the vortex rings produced during the upstroke, downstroke and stroke transitions. These vortex rings help dolphin's forward propulsion by producing backward facing jets. The downstroke produces stronger jets on average, leading to higher downstroke thrust and net lift over a cycle of dolphin-like oscillation. These results correspond well with the empirical observations by Fish et al [13,14,25] of a staggered array of vortex rings in the wake of a swimming dolphin, with visualization of the vortex rings and the propulsive jets done using PIV by tracing the position of microbubbles as the dolphin swam along a curtain of bubbles [25].
The results of the present study attributed the majority of thrust production to the oscillation of the flukes, with the posterior body supporting the thrust production. The forces acting on the flukes were a lot greater than most areas of the body. The flukes produce thrust by generating leading-edge vortices, which develops a leading-edge suction region that pulls the flukes forward and produces thrust. The shedding of the leading-edge vortices and trailingedge vortices lead to the thrust-generating vortex rings.
Due to the presence of the ball vertebra near the anterior insertion of the flukes on the peduncle, the increased bending of the peduncle can effectively camber the flukes [24]. Curvature measurements show that the peduncle joint experiences significant flexion, and flukes themselves experience high cambering. Dolphins in different maneuvers exhibit varying amounts of fluke flexion. Varying the flexion of the peduncle and flukes, we produced alternative dolphin-inspired oscillatory kinematic patterns to study the effect of fluke flexion on propulsive performance. Fluke flexion scaled down from the original benchmark was found to lead to 11.2% higher cycle-average thrust production, by generating stronger leading-edge vortices and stronger trailing-edge vortices that can generate stronger jet, and stronger leading-edge suction. Increased flexion can help with propulsive efficiency by reducing the production of non-axial forces; the orientation of the flukes allows it to better cut through the flow, and the leading-edge vortex is better attached to the fluke's surface, leading to, though weaker vorticity and thrust jets, a larger proportion of the force being directed forward, and thus enhancing the propulsive efficiency.

Data availability statement
All data that support the findings of this study are included within the article (and any supplementary files).