Propulsion performance of a skeleton-strengthened fin

Skeleton-reinforced bio-membranes are ubiquitous and play critical roles in many biological functions. For example, the flexibility of the wings of many insects is determined primarily by the architecture of the embedded veins. By experimentally measuring and numerically characterizing the wing structure of 16 different species of insects, Combes and Daniel concluded that an anisotropic flexibility is achieved through the architecture of these composite structures (Combes and Daniel,2003a; Combes and Daniel,2003b). In all the species they studied, the spanwise bending stiffness was found to be at least 1–2 orders of magnitude higher than the chordwise bending stiffness (see also Newman and Wootton, 1986; Wootton, 1992). It was attributed mostly to the existence of clustered or thickened veins near the leading edge. This suggests that, through specific arrangement of supporting veins, these insect wings allow certain passive deformations while discouraging others. Structures with similar architecture are also found in the fins of aquatic animals. Many fishes possess caudal, dorsal, pectoral or ventral fins that are stiffened by fin rays embedded in a thin layer of collagenous membrane (Lauder,1989; Lauder and Drucker,2004; Fish and Lauder,2006; Lauder and Madden,2007). Overall, there are over 25 000 species of fish with these ray-strengthened fins. According to morphological studies (Harder, 1971; Kardong, 1998), a fin ray contains a central bundle of collagen surrounded by small segmented bony elements called hemitrichs. These elements are paired and resemble a bimetallic strip with two elongated bony elements separated by the central collagen core. Along the length of a hemitrich, and at the ends, there exist short ligaments and elastic fibers. In addition, the basal end of each ray attaches to four separate muscles. This architecture enables a fish to control the motion of each individual ray and thus achieve three-dimensional fin deformations. This property is expected to have fundamental effects on the hydrodynamic performance of the fins in generating propulsion force during cruising, bursting and maneuvering.

The exact function of the anisotropic flexibility of skeleton-strengthened membranes in insect flying and fish swimming remains to be fully understood. Despite a limited number of investigations(Cheng et al., 1998; Pedley and Hill, 1999),existing studies in fish swimming and insect flying are predominantly based upon the assumption that the wings and fins are rigid(Liu and Kawachi, 1998; Ramamurti and Sandberg, 2002; Ramamurti and Sandberg, 2007; Wang, 2000; Pullin and Wang, 2004) or with prescribed deformations (Liu and Bose,1997; Ramamurti et al.,2005; Mittal et al.,2006). It has always been speculated that the structural flexibility of these composite bio-structures, with their stiffness determined by the coupling of the embedded skeletons and the surrounding materials,attributes to the highly efficient force generation through flapping motions. Evidence of the beneficial effects of structural flexibility in the generation of lift and propulsion forces comes from different sources, in particular the analytical, numerical and experimental investigations of flapping foils, as discussed below.

Unlike conventional aero- and hydro-foils, which utilize steady fluid dynamics for force generation, flapping foils are modeled on fish fins and insect wings so that they produce aerodynamic or hydrodynamic lift, propulsion or maneuvering forces using unsteady foil oscillations(Triantafyllou et al., 1991; Tuncer and Platzer, 1996). Through various investigations, it has been shown that structural flexibility, per se, increases the propulsion efficiency in many cases. For example, by experimentally studying the performance of a semi-flexible flapping foil, Yamamoto et al. reported up to a 27% increase in propulsion efficiency compared with a rigid foil(Yamamoto et al., 1995). Similar performance enhancement was observed by Heathcote et al. in their investigation of a thin steel plate undergoing periodic heaving motion in still water (Heathcote et al.,2004). By theoretically examining the oscillation of a two-dimensional flexible plate, Katz and Weihs were able to obtain a 20%increase in propulsion efficiency (Katz and Weihs, 1978). The performance of a flexible foil, by contrast,has been found to depend not only on its elasticity but also on its inertia(Pederzani and Haj-Hariri,2006). This conclusion was confirmed in our recent work(Zhu, 2007). In that study, we investigated the three-dimensional fluid–structure interaction of a flapping foil. A number of extreme, yet representative, cases were examined,including chordwise vs spanwise flexibility and fully-coupled (the foil deformation is driven mostly by the fluid) vs partially decoupled (the foil deformation is driven mostly by its own inertia)interactions. By using a boundary-element model coupled with a two-dimensional plate model, we were able to show that both the structural flexibility and the inertia of a flapping foil had significant effects on its capability for thrust generation. Among all the numerical tests we have conducted, two cases– one with chordwise flexibility and fluid-driven deformation and the other with spanwise flexibility and inertia-driven deformation – yield higher propulsion performance than rigid foils. The former case shows increased efficiency (by 20%, consistent with previous studies) and the latter shows increased thrust (by almost twofold in some cases) without any significant loss in efficiency.

Structural flexibility has also been found to effectively reduce fluid drags. Studying the two-dimensional soap film flow around an elastic glass fiber as an example, Alben et al. characterized the drag reduction effect of structural flexibility, illustrating the transition from a rigid-body drag scaling law of U² to U^4/3 (Uis the speed of incoming flow) when the fiber flexibility exceeds a threshold(Alben et al., 2002; Alben et al., 2004).

Despite the abovementioned efforts to characterize dynamics of flexible foils or fibers, the hydrodynamic significance of three-dimensional anisotropic flexibility is still not illustrated. In particular, comprehensive studies are required to understand the coupled fluid–structure interactions of bio-membranes with morphological relevance, e.g. those with mechanical properties similar to skeleton-reinforced fins or wings.

We herein present a numerical investigation by coupling fluid dynamics and structural mechanics to examine the propulsion performance of a skeleton-reinforced fin that is geometrically and structurally similar to the caudal fin of a fish. We study cases in which the unsteady motions of the rays are individually controlled by imposing sway–yaw motions at their basal ends. The rest of the ray may deform under the combined effect of its own inertia, the hydrodynamic load and constraint from the membrane. Our approach includes an Euler–Bernoulli beam model for the strengthening rays, as well as a boundary-integral equation representation for the surrounding fluid. To account for large three-dimensional deformations, the beam model is fully nonlinear and allows stretching, bending and twisting. Via systematic simulations, dynamics of a flexible fin and that of a rigid one are compared with each other using the performance metrics, including the capacity and efficiency of the fin to generate thrust and the reduction of transverse force. In addition, we investigate both the homocercal mode, in which the fin deformation is symmetric across the span, and the heterocercal mode, in which the deformation of the dorsal side of the fin is different from that of the ventral side.

The rest of the paper will be organized as follows. In the next section we define the geometry and internal structure of the fin, together with the mathematical formulations and the numerical methods to solve the coupled fluid–structure interaction problem. Following that, numerical results,including the thrust force, the transverse force and the propulsion efficiency of the flexible fin in comparison with a rigid fin, will be presented. Also included will be comparisons between the homocercal and the heterocercal modes. Finally, conclusions will be drawn.

As a model of the real ray-strengthened fin(Fig. 1A), we consider a geometrically and structurally simplified fin as shown in Fig. 1C. This trapezoid fin has chord length c and leading edge span. In a Cartesian coordinate system (x, y, z), the center plane of the fin at its mean position locates within the x–z plane and is symmetric with respect to the x–y plane. This fin is reinforced by nine embedded rays, numbered 1 to 9. The basal ends of the rays are distributed evenly on the leading edge, as are the angles they form with respect to the x-axis, ranging from 30° (Ray 1) to –30°(Ray 9). The thickness of the fin, d, is assumed to be uniform except for tapering at the leading and the trailing edges (note that the value of d is only relevant in the hydrodynamic problem). To put it in perspective, this configuration is similar to the homocercal caudal fin of species such as bluegill sunfish (Lepomis macrochirus).

Despite the fact that the real fin rays are non-uniform, anisotropic and have curvatures actively controlled by offsets of tendons linked to their bases (see Fig. 1A,B), for simplicity in the present investigation we model the rays as uniform Euler–Bernoulli beams with circular cross-sections with passive deformations only. By this definition, these rays are able to sustain stretching, bending and twisting loads. We further idealize the rest of the fin as a membrane that can sustain stretching/compression but not bending and replace it by distributions of spring-dampers between neighboring rays. Coupled with this beam–spring structure, we model dynamics of the surrounding fluid by using potential flow approximation.

We consider a three-dimensional flapping motion of the fin in combination with a constant forward speed, U, in the –x direction. The flapping motion is actuated by imposing sway–yaw motions on each ray. By assumption, the rays share the same sway motion y(t)=y₀sin(ωt) following the motion of the caudal peduncle. The yaw motion with respect to the z-axis of Ray i isθ(t)=θ_isin(ωt–ψ_i), i=1,..., 9. The Strouhal number, a non-dimensional parameter characterizing the unsteady fluid-dynamic effects, is defined to be S_t=2y₀f/U, where f=1/T=ω/2ο is the frequency of the motion and T is the period.

In the following, we elaborate the mathematical formulations of the ray motion and the fluid dynamics.

We apply a three-dimensional nonlinear formulation and model each ray as an Euler–Bernoulli beam that could be stretched, bent and twisted(Tjavaras et al., 1998). Based upon the Euler–Lagrangian dual-coordinate approach and the robust Euler parameters, this method achieves fully-nonlinear simulation of arbitrary ray deformations. A finite difference algorithm, the modified box method, is then employed to solve these nonlinear equations (see Appendix A for details).

The flow around the fin is assumed to be irrotational except for zero-thickness sheets of vorticity shed from the trailing edge so that it can be described by a flow potential Φ(x,t), where x=(x, y, z) is the position vector. The near-body flow field as well as the hydrodynamic loading on the fin are then predicted using a boundary-element approach based upon the boundary-integral framework by segmenting the fin surface S_b into N_bpanels S_bj, j=1,..., N_b(Appendix B).

As mentioned before, the structural function of the membrane in which the rays are embedded is simplified as linear spring-dampers spanning between neighboring rays. Specifically, in the following calculations, we chose the spring constant to be 20 N m^–2, and the damping constant to be 2 Ns m^–2 (both are measured per unit length). With such a constraint, in the following simulations the maximum variation of the surface area of the fin is around 5–10%.

We employ an iteration method to solve the fully coupled fluid–structure interaction problem. Below we list the primary steps of this algorithm.

1. At time t, we start from an initial guess,

   \(S_{\mathrm{b}}^{(0)}\)
   ⁠, of the configuration of the fin.

2. The hydrodynamic problem is solved based on

   \(S_{\mathrm{b}}^{(0)}\)
   ⁠.

3. After the hydrodynamic pressure, p, within panel j is evaluated using Eqn 1, we assume that the hydrodynamic force acting on this panel, pnΔS_bi(ΔS_bi is the area of the panel), is transferred equally to four grid points that are closest to the panel on the two neighboring rays (two points on each ray). Summing up the forces on all the panels, the hydrodynamic loads, F_h, on the nine rays are thus determined.

4. The Bernoulli–Euler beam equations(EqnA9) are numerically integrated to obtain the updated positions of the rays. Based upon these, the updated configuration of the fin, S_b, is determined via linear interpolation.

5. If S is not sufficiently close to

   \(S_{\mathrm{b}}^{(0)}\)
   ⁠, we set
   \(S_{\mathrm{b}}^{(0)}=S_{\mathrm{b}}\)
   and repeat Steps 2–4.

To eliminate the well-known added mass induced instability, we apply an implicit added mass scheme. In this approach, the added mass effects are subtracted from both sides of Eqn A9 during time integration. As long as convergence of the iteration is achieved, the exact values of the added mass will not affect the results. A detailed description of this method is provided in Connell and Yue(Connell and Yue, 2007).

The validity and accuracy of both the Euler–Bernoulli beam solver and the boundary integral solver have been extensively tested and documented through numerous convergence studies and comparisons with experiments or other numerical/theoretical predictions(Tjavaras et al., 1998; Zhu et al., 2002; Zhu et al., 2006; Zhu, 2007). These are omitted here for brevity. Instead, we will concentrate primarily on the convergence of the algorithm with respect to the time step, as it reflects the accuracy of our iteration scheme discussed above, which is the core of the interaction solver.

Convergence tests are performed by examining the propulsive force generated by a fin with leading-edge span s₀=0.02 m, chord length c=0.1 m, and thickness d=0.002 m. The rays are structurally identical except for their lengths. Each has a cross-sectional diameter of 1 mm, a Young's modulus, E, of 0.4 GPa, and a hysteretic (material)damping coefficient, D, of 0.4 sGPa. The density of the rays is chosen to be the same as that of the water, i.e. 10³ kg m^–3.

Kinematically, we set the forward speed to be U=0.1 ms^–1. The sway amplitude, y₀=0.5c. The Strouhal number, S_t=0.2. We also assume that the nine rays shear the same yaw motion, with θ_0i=θ₀=10⁰ andψ _i=90° (i=1,..., 9).

Table 1 shows the mean thrust coefficient, F̄_t/½ρU²c²,as a function of the number of body elements, N_b, and the time step Δt. Linear convergence is achieved with respect to both parameters.

Lauder (Lauder, 2000)studied kinematics of caudal fins of three species: the leopard shark(Triakis semifasciata), the sturgeon (Acipenser transmontanus) and the bluegill sunfish (Lepomis macrochirus). Two of these fishes, the shark and the sturgeon, possess caudal fins with dorso-ventral asymmetry (heterocercal fins). The bluegill sunfish, on the other hand, has caudal fins that are geometrically symmetric across the span(homocercal fins). Surprisingly, his observations illustrated that, despite the geometric difference between homocercal and heterocercal fins, their motions are nevertheless heterocercal in many cases. Specifically, the dorsal side of the fin was found to be ahead of the ventral side. In such a heterocercal mode, a small vertical force is generated, which has to be balanced by forces generated by the anterior part of the fish body. The heterocercal motion was attributed to the embedded muscle through which the kinematics of the rays can be controlled. Similar heterocercal modes in homocercal fins have been found on chub mackerel (Scomber japonicus)(Gibb et al., 1999). On the other hand, in another experiment of a bluegill sunfish, it was observed that the caudal fin performed mostly homocercal motions(Tytell, 2006). The implication is that a fish with a homocercal caudal fin may be able to switch between different swimming modes.

By specifying the motion of each individual reinforcing ray, in our model we are able to achieve controlled anisotropic deformations, both homocercal and heterocercal, of the fin. We assume that the basal ends of the rays share the same sway motion following the undulation of the anterior part of the fish body. The yaw motions, i.e. the rotation of the rays around their basal ends,can be varied individually. This imitates the muscle action that controls the ray kinematics (Lauder and Drucker,2004). In the following, we will present our numerical results of the force generation by a rigid fin and a flexible fin with identical geometries. The flexible fin has the same geometric and structural properties as the one applied for the convergence tests. In addition, we assume that in all the simulations the forward speed, U, is 0.1 m s^–1, the sway amplitude, y₀, is 0.05 m,and the phase angle between sway and yaw, ψ_i, is 90°(i=1,..., 9). The focus of our numerical simulation will be the comparison between force generation by a rigid fin and that by a flexible fin. We will also study the differences and similarities between homocercal and heterocercal modes.

Numerically, 1000 boundary elements are applied on the fin surface. On each ray, the number of grids is N_r=32. Finally, a time step ofΔ t=T/128 is chosen.

In the homocercal motion, we assume that the yaw amplitudes of all the rays are the same, i.e. θ_0i=θ₀ (i=1,...,9). A typical sequence of the deformation of a flexible fin over ¼period is displayed in Fig. 2. At t=0, when the fin is at its mean sway position (y=0) and maximum yaw angle (θ=10°), we observe that there is no pronounced spanwise (dorso-ventral) bending. As t increases towards T/8 and T/4, accompanying the increasing sway, significant dorso-ventral bending is observed, creating a bow-shaped trailing edge so that the central part of the fin bends into the flow. If we plot the trajectories of the trailing ends of Ray 5 (at the center) and Ray 1 (at the upper edge), it is clear that the upper and lower edges of the fin undergo smaller and lagged-behind undulations compared with the center part(Fig. 3). In Fig. 2, we also observe that the dorso-ventral bending reverses its direction near the upper and lower edges. There are consequently two inflection points.

The shape of the trailing edge during flapping is attributed to the combined effect of the anisotropic flexibility of the fin as determined by the embedded rays and the hydrodynamic loads on its surface. Firstly, the central rays are shorter in length than those closer to the upper and lower edges and thus sustain less bending during the flapping, leading to a bow-shaped trailing edge. This effect is reinforced by the pressure distribution on the fin surface during its flapping motion. As shown in Fig. 4, on the fin surface we see areas of large negative pressure near the upper and lower edges, caused by large flow velocity near the upper and lower edges, resembling flow around acute corners. The suction force created by this negative pressure creates larger than normal bending on the rays close to the two edges (e.g. Ray 2 and Ray 8). By contrast, the rays located at the upper and lower edges (Ray 1 and Ray 9) are only attached to the membrane from one side (i.e. these rays are`semi-free'), so the hydrodynamic loading on them is less than that on their neighboring rays. This mechanism reduces the bending of the edge rays in comparison with their neighboring rays and causes the inflection points in the displacement curve of the trailing edge.

In reality, vorticity shedding from the upper and lower edges, which is not included in our current model, may affect the near-body flow field as well as the hydrodynamic loading on the fin. The strength and behavior of these edge vortices depend on kinematic parameters such as Reynolds number (Re),frequency, sway amplitude and angle of attack (e.g. Triantafyllou et al., 2004; Birch et al., 2004). These vortices are expected to cause quantitative (but not qualitative) changes in the predicted pressure distribution on the fin surface. Indeed, on delta wings, the concentrated negative pressure near the edges has been observed experimentally (e.g. Furman and Breitsamter, 2006).

The bow-shaped trailing edge of the fin, as well as the existence of inflection points near the upper and lower edges, is qualitatively consistent with a recent report about kinematics of the caudal fin of a bluegill sunfish(Tytell, 2006) (see fig. 3 in that paper).

Figs 5, 6, 7 summarize the mean thrust force, the propulsion efficiency and the amplitude of the transverse force generated by both a rigid fin and a flexible fin within the parameter range of 0.1<S_t<0.4 and–10°<θ₀<40°. From Fig. 5 we see that the flexible fin generates a larger thrust when–10°<θ₀<∼5°. Otherwise, the rigid fin creates a larger thrust. With regard to the propulsion efficiency, the flexible fin outperforms the rigid fin through the whole range(Fig. 6). The maximum propulsion efficiencies achieved by the rigid fin and the flexible fin are found to be 0.48 and 0.54, respectively. More importantly, the flexible fin achieves good performance within a much larger parameter range than the rigid fin. Indeed, at certain places (e.g. S_t=0.3,θ ₀=–10°), almost a 100% increase in propulsion efficiency is observed (also refer to Fig. 8). Thus, we conclude that by using a flexible fin, we can greatly reduce the dependence of its propulsion performance on kinematic parameters.

Another advantage of the flexible fin is that it significantly reduces the transverse force (Fig. 7). According to our simulations, an estimated 10–40% (case dependent)reduction in the transverse force coefficient is recorded. The implication is that less power is required to actuate the sway motion and is therefore consistent with the increased efficiency.

Both the increase in propulsion efficiency and the decrease in transverse force indicate performance enhancement through three-dimensional structural flexibility. These characteristics are more clearly illustrated in Fig. 8, where cross-sections of Figs 5, 6, 7 are plotted.

In our simulations, the heterocercal mode, i.e. a mode with dorso-ventral asymmetry, is achieved by varying the amplitudes (but not the frequency and phase) of yaw motions of the rays. Specifically, we assume that the ray at the upper edge on the dorsal side (Ray 1) performs no yaw motion(θ₀₁=0°). For the rest of the rays, the yaw amplitude increases linearly with the largest yaw motion occurring at the lower edge(Ray 9), i.e. θ₀₁=i/9×θ₀. Hereafter, we take a closer look at a flexible fin.

As shown in Fig. 9, the ray at the upper edge (Ray 1) leads both the central ray (Ray 5) and the ray at the lower edge (Ray 9). This trend is more clearly shown in Fig. 10, where the trajectories of the trailing ends of these three rays are plotted. Such behavior is consistent with experimental observations of bluegill sunfish as reported by Lauder (Lauder,2000).

Fig. 11 displays the capacity of force generation by this flexible fin in a heterocercal mode within the range 0.1<S_t<0.4 and 10°<θ₀<50°. As we see, compared with the homocercal mode, the heterocercal mode yields a slightly smaller propulsion efficiency. For example, the maximum propulsion efficiency achieved in the heterocercal case is 0.51, in comparison with 0.54 in the homocercal case.

An important advantage of the heterocercal motion over the homocercal motion may be that it further reduces the transverse force(Fig. 11C). For example, by comparing the two motions at fixed Strouhal numbers but varying yaw angles, we find that the minimum transverse force is 5–10% smaller in the heterocercal mode. This is attributed to the diminished spanwise coherence of the near-body flow in the heterocercal mode.

Fig. 11D shows the existence of a mean lifting force as generated by the heterocercal motion of the fin. Although this force is 1–2 orders of magnitude smaller than the thrust, it may be an important source of maneuvering force in the vertical direction.

To further investigate the fluid–structure interaction problem involved in the sway/yaw motion of the fin, and to compare the homocercal mode with the heterocercal mode, we numerically visualized both the near-body flow field and the vorticity distribution in the wake.

As illustrated in previous studies, an important phenomenon concerning flexible objects in flow is that these objects often bend their surfaces to match the incoming stream. This has been confirmed in various experiments,including, for example, the bending of a glass fiber by a soap film flow(Alben et al., 2004) and the meandering motion of a dead fish between incoming vortices(Liao et al., 2003). To confirm if this is the case in our problem, in Fig. 12 we checked the inplane streamlines within two horizontal planes, z=0 and z=0.3c, around a rigid fin and a flexible fin undergoing sway motions. It was found that at the mid-span (z=0), the fin is bent just a little towards the incoming flow. This results from the fact that the central ray (Ray 5) is the shortest and least deformable. On the plane z=0.3c, however, two different tendencies are observed. The front part of the fin tends to align with the incoming flow, echoing features of the two-dimensional fluid-driven foil deformations(Zhu, 2007). The hind part of the fin, however, undergoes much smaller reorientation. This is caused by the three-dimensional effect, in particular interactions with rays closer to the center via the membrane.

Fig. 13A,B shows the iso-surfaces of vorticity in the wake behind a rigid fin and a flexible fin undergoing a homocercal flapping motion, respectively. It is evident(especially from the insets of Fig. 13A,B) that the wake consists of a sequence of vortex rings symmetric in the dorsal/ventral direction (in two-dimensional experimental visualizations or numerical simulations, these vortex rings appear as pairs of vortices). The wakes of the rigid fin and the flexible fin resemble each other except for a subtle difference at the connection between neighboring rings. Specifically, in the flexible case, the figure shows elongated forklike structures between neighboring rings (Fig. 13B). These structures resemble the hairpin structures proposed by Tytell (Tytell, 2006) based upon analysis of PIV data obtained from the near-body flow around a bluegill sunfish. The hairpins were attributed to the cupping of the fin during flapping, as is also shown in our simulations (see Fig. 2).

Fig. 13C shows an iso-surface of vorticity when the flexible fin performs a heterocercal oscillation (the wake behind a rigid fin does not show significant difference,and is thus omitted here). It is seen that the dorso-ventral symmetry of the wake is broken, resulting in a chain of tilted vortex rings. This is consistent with predictions based upon PIV pictures obtained from experimental visualizations of the wake behind a bluegill sunfish in heterocercal mode(Lauder, 2000). The typical size of these vortex rings is comparable to the height of the fin, which confirms another characteristic predicted via PIV data (the smaller ring at the end of the chain is generated during the start-up period, in which the sway/yaw amplitudes of the fin are increased smoothly from zero to the requested values) (see also Nauen and Lauder, 2002). Since these vortex rings are tilted backwards, they induce a downwash flow in the wake and create the lifting force as observed above. This downwash is evident in Fig. 14B, where the flow field within the y=0 plane is shown. By plotting the y-component of the vorticity and the inplane velocity vectors, it is clear that each vortex pair (as a cross-section of a vortex ring) induces between them a jet flow pointing slightly downwards. Such a flow field is very similar to the DPIV image of the inplane flow behind a chub marckerel in Nauen and Lauder (Nauen and Lauder, 2002) (see fig. 4 in that paper). For comparison, in Fig. 14A we also plot the inplane flow field behind a fin undergoing homocercal motion. The figure shows counter-rotating vortex pairs with outgoing jets induced between them.

By using a fluid–structure interaction model combining a nonlinear Euler–Bernoulli beam model and a boundary-integral hydrodynamics model,we have investigated the force generation by a flexible fin made of a thin membrane reinforced by embedded rays. This structure resembles the ray-strengthened fins of many fishes. Structurally, the distribution of reinforcing rays creates a anisotropic flexibility. Two different locomotion modes, the homocercal (dorso-ventral symmetric) motion and the heterocercal(dorso-ventral asymmetric) motion, are studied.

Our results show that compared to rigid fins with the same geometry, a fin with anisotropic flexibility yields higher propulsion efficiency. The structural flexibility also reduces the transverse force, another beneficial effect in fish swimming. Further reduction of transverse force is observed by applying the heterocercal mode.

An interesting phenomenon is that in terms of propulsion efficiency, the fluid–structure coupling greatly reduces the sensitivity of the fin performance to kinematic parameters such as the Strouhal number and the amplitude of yaw motion. Even in cases when these parameters are not optimized, the flexible fin is able to function well as a propeller. This property may have significant effects on fish locomotion. The biological implication is that it may not be critical for a fish to control the motion of its caudal fin accurately to achieve high propulsion performance. In biomimetic applications, this suggests that in the future development of mechanical devices imitating the swimming mechanism of fishes, the controlling system could be significantly simplified provided that the structural properties of the propeller are designed to mimic the ray-reinforced fins.

In addition to the locomotion modes examined in our work, there may exist other combinations of fin structures and kinematics, leading to different fin deformations in the swimming process. For example, by observing the fin actuation of a dace (Leuciscus leuciscus), Bainbridge concluded that the maximum dorso-ventral curvature occurred when the caudal fin is close to its maximum sway displacements (Bainbridge,1963). His observations also show that the upper and lower edges of the fin usually precede the center part in their lateral motions. This contradicts what we found in the homocercal mode when the upper/lower edges are behind the central part. Similarly, although the fin kinematics predicted by our model appear to resemble the experimental observations of Tytell(Tytell, 2006), subtle differences, such as the relationship between the center part and the upper/lower edges in terms of amplitudes of deflection and phases, suggest more complicated fin kinematics in his experiment than in our simulations. In species with heterocercal tails, such as spiny dogfish (Squalus acanthias), a sophisticated ring-within-a-ring vortex structure has been observed (Wilga and Lauder,2004). This was attributed to the asymmetry in vortex formation at the upper and the lower lobes. A more sophisticated ray model, with actively controlled curvatures, is also required. To study cases with closer biological relevance, and to fully understand the benefits of different fin architectures and locomotion modes, interdisciplinary investigations including fish morphology and physiology, fish behavior, control, as well as fluid–structure interactions are required.

The perspective gained from our work and previous studies makes it appropriate to envisage future development of biomimetic propellers imitating flexible flapping fins made of skeleton-reinforced membranes. Compared with the conventional rigid foil design, these novel biomimetic propellers possess the following advantages.

1. Enhanced performance: as illustrated in our study, a flexible foil with three-dimensional deformability yields better performance than a rigid one,manifested in the increased efficiency, reduced sensitivity to kinematic parameters and diminished force in the transverse direction.

2. Controllability: with strengthening skeletons, the structural properties of the composite membrane can be adjusted in different operation modes. For example, our results indicate that a more rigid fin may generate larger thrusts in a certain parameter range. Therefore, in burst swimming, it may be desirable to increase the stiffness of the supporting rays.

3. Structural strength and lightness: the biomimetic composite membrane provides a light structure with high strength.

4. Deployability: the skeleton-reinforced membrane resembles deployable structures such as cable roof buildings. These structures can be easily folded and unfolded, making them an ideal design as a portable device.

It is necessary to point out the limitations of the current model. First,our work concerns a geometrically and structurally simplified caudal fin. Effects of the fish body and attachments such as dorsal, pectoral or anal fins are not considered. As illustrated by Wilga and Lauder(Wilga and Lauder, 2002), in leopard sharks (Triakis semifasciata) and bamboo sharks(Chiloscyllium punctatum), the inclination angle of the body plays a pivotal role in balancing the torque created by the tail. In addition, as shown by Tytell (Tytell,2006), vortices shed from the dorsal and the anal fins of a bluegill sunfish may contribute significantly to the near-body flow. Recent experiments with a mechanical fish have shown that the fish body itself might also generate a sequence of weaker vortex rings in the wake(Brucker and Bleckmann, 2007). The dynamic effects of the vortices generated upstream of the caudal fin depend not only on their strength but also on their phase with respect to the caudal fin (Zhu et al., 2002). This, in turn, is determined by the exact morphology of the fish as well as its kinematics. Another effect not considered in the present study is vorticity generation from the fin surface in locations other than the trailing edge, in particular the upper and lower edges. This vorticity generation will not only change the local dynamics (e.g. the pressure distribution) around the upper/lower edges but will also affect the overall performance of the fin in propulsion. To fully account for this effect, a Navier–Stokes solver capable of studying three-dimensional fluid–structure interactions is required [for a possible method to achieve this, see Connell(Connell, 2006)]. A fully viscous study will not only illustrate effects of vortex generation from places other than the trailing edge but will also provide more accurate predictions of the flow field and the hydrodynamic forces, especially in cases when the Reynolds number is low. These Navier–Stokes simulations,however, are computationally expensive. In that respect, a potential flow method like ours provides an alternative way to conduct systematic tests over a wide range of parameters. To improve the accuracy of our current method,investigations are underway to develop methods to study vortex generation from the leading edge area or the upper/lower edges of the fin. The knowledge gained from these studies will pave the way for future simulations with more comprehensive numerical tools and experiments using mechanical devices. These studies may eventually contribute to the creation of biomimetic apparatus that will be installed on AUVs (autonomous underwater vehicles) or other vehicles(see Tangorra et al.,2007).

Finally, we would like to mention that, in addition to insect wings and fish fins, there exist a large number of bio-structures in nature with similar structural characteristics, which can thus be categorized as skeleton-reinforced membranes. A typical example is found in the membrane of the erythrocyte (red blood cell), which features a composite structure including a lipid bilayer and a protein skeleton consisting primarily of actin and spectrin (Mohandas and Evans,1994). The combined visco-elasticity of the bilayer and the skeleton, as well as their dynamic inter-connectivity, impart remarkable stability and deformability essential for the cell functionality when it circulates around the body and squeezes through capillaries half its own diameter. Similar structures are found in the biopolymer membranes within the shells of mollusks such as the red abalone (Haliotis rufescens) and the nautilus (Nautilus pompilius). Remarkable mechanical properties of these structures (e.g. structural strength, stability, durability and deployability) suggest biomimetic applications and justify concentrated research efforts to understand the detailed structural mechanics and fluid–structure interactions of these bio-structures.

Overall, 13 boundary conditions are required at the two ends of each ray. At the basal end, we apply: (1) β₀=0; (2) the velocities V_τ, V_n and V_b, as determined by the prescribed sway motion; (3) M_τ=0; and (4) the three-dimensional orientation of the ray as determined by the prescribed yaw motion. In practice, condition 4 is enforced by adding a stiff rotational spring between the basal end of the ray and a reference direction representing the yaw motion. This leads to moments M_n and M_b at the basal end, whose values increase as the ray deviates from the reference direction. Six additional boundary conditions are applied at the trailing end, including: (1)

An iterative algorithm is employed to solve this system of implicit equations (Tjavaras et al.,1998). The method possesses second-order accuracy in space. The accuracy of time integration, on the other hand, is first order. The purpose of selecting a lower-order time integration scheme rather than the more conventional box method, which possesses second-order accuracy, is to providing better numerical stability.

The flow around the fin is assumed to be irrotational except for zero-thickness sheets of vorticity shed from the trailing edge. The flow field is then easily described by a velocity potential Φ(x,t)(the position vector x≡(x, y, z) represents an arbitrary point in space). We further express Φ as the linear superposition of two parts, Φ=ϕ_b+ϕ_w, where ϕ_brepresents the contribution from the fin body and ϕ_w represents the contribution from the wake. The problem is thus decomposed into two interdependent ones: the boundary value problem for ϕ_b and the evolution of the wake that determines ϕ_w.

To determine ϕ_w, we model the wake as a distribution of dipoles on a sheet shed from the sharp trailing edge of the fin. With the continuous generation of vorticity from the trailing edge, this wake sheet elongates over time. The strength of newly shed wake is determined by the Kutta condition, i.e. the strength of the newly shed wake ϕ_wequals the instantaneous difference of body influence potentialϕ _b between the upper and lower surfaces near the edge so that the local singularity is eliminated. The rest of the wake, treated as a material surface, is carried downstream by the collective effect of the upcoming flow and the self-induced velocity, and its strength is kept unchangeable due to the absence of dissipation. In order to stabilize the evolution of the discretely distributed vorticity in the wake, it is necessary to introduce a desingularization algorithm, in which we assign a finite core radius δ to each (body or wake) panel of dipole distribution(Krasny, 1986). It has been demonstrated that the hydrodynamic forces on the body are not sensitive to the value of δ (Zhu et al.,2002). The desingularization also eliminates singularity in the wake vorticity distribution, allowing us to visualize the near-body flow field.

A low-order boundary-element algorithm is employed to find the unknown body influence potential ϕ_b on S_b by solving Eqn B2(Katz and Plotkin, 1991). To achieve this, the fin surface S_b is discretized into N_b panels. The low-order panel method refers to the fact that in this algorithm within each panel S_bj(j=1,..., N_b), the body potentialϕ _b and its normal derivative are assumed to be constants. By collocating at the centroids of the panels, Eqn B2 is transformed into a system of linear equations and can be numerically solved. Afterwards, a forward Euler algorithm is applied for time integrations. As an illustration,the distribution of boundary elements on the fin surface and the wake is demonstrated in Fig. B1. The flow potential Φ is then determined by summing up the body influence potential ϕ_b and the wake influence potentialϕ _w.

LIST OF SYMBOLS AND ABBREVIATIONS

 
• c

  chord length of the foil

 
• C

  rotation matrix

 
• d

  thickness of the fin

 
• D

  material damping coefficient of the rays

 
• E

  Young's modulus of the rays

 
• EI

  bending stiffness of the rays

 
• F

  hydrodynamic force on the fin

 
• F_h

  hydrodynamic force on rays

 
• F_r^(r)

  first-harmonic transverse force

 
• F_t(F̄_t)

  (mean) thrust force

 
• GJ

  torsional stiffness of the rays

 
• (i, j, k)

  unit vectors of the global coordinate system

 
• N_b

  number of boundary elements on the fin surface

 
• N_r

  number of grids on a ray

 
• p

  hydrodynamic pressure

 
• P(P̄)

  (mean) power spent

 
• s

  distance from a point on a ray to the basal end measured along the unstretched ray

 
• s₀

  span at the leading edge

 
• S_t

  Strouhal number

 
• T

  period

 
• T=(T_τ, T_n, T_b)

  internal force inside a ray

 
• U

  forward speed

 
• V=(V_τ, V_n, V_b)

  translational velocity of a ray

 
• (x, y, z)

  global coordinate system

 
• y₀

  sway amplitude

 
• β0, β1, β2, β3

  Euler parameters

 
• ϵ

  strain in a ray

 
• Δt

  time step

 
• η

  propulsion efficiency

 
• θ0i

  yaw amplitude of the i-th ray

 
• (τ, n, b)

  tangential, normal and bi-normal unit vectors in the Lagrangian coordinate system

 
• ϕb

  potential due to the fin

 
• ϕw

  potential due to the wake

 
• Φ

  flow potential

 
• ψi

  phase between sway and yaw of the i-th ray

 
• θ

  frequency

 
• ω=(ωτ, ωn, ωb)

  rotational velocity of a ray

 
• Ω=(Ωτ, Ωn, Ωb)

  Darboux vector measuring material torsion and curvatures of a ray