SUMMARY
We examine numerically the performance of a thin foil reinforced by embedded rays resembling the caudal fins of many fishes. In our study, the supporting rays are depicted as nonlinear Euler–Bernoulli beams with threedimensional deformability. This structural model is then incorporated into a boundaryelement hydrodynamic model to achieve coupled fluid–structure interaction simulation. Kinematically, we incorporate both a homocercal mode with dorsoventral symmetry and a heterocercal mode with dorsoventral asymmetry. Using the homocercal mode, our results demonstrate that the anisotropic deformability of the rayreinforced fin significantly increases its capacity of force generation. This performance enhancement manifests as increased propulsion efficiency, reduced transverse force and reduced sensitivity to kinematic parameters. Further reduction in transverse force is observed by using the heterocercal mode. In the heterocercal model, the fin also generates a small lifting force, which may be important in vertical maneuvers. Via threedimensional flow visualization, a chain of vortex rings is observed in the wake. Detailed features of the wake, e.g. the orientation of the vortex rings in the heterocercal mode, agree with predictions based upon particle image velocimetry (PIV) measurements of flow around live fish.
INTRODUCTION
Skeletonreinforced biomembranes 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 raystrengthened 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 threedimensional 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 skeletonstrengthened 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 biostructures, 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 hydrofoils, 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 semiflexible 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 twodimensional 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 HajHariri, 2006). This conclusion was confirmed in our recent work (Zhu, 2007). In that study, we investigated the threedimensional fluid–structure interaction of a flapping foil. A number of extreme, yet representative, cases were examined, including chordwise vs spanwise flexibility and fullycoupled (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 boundaryelement model coupled with a twodimensional 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 fluiddriven deformation and the other with spanwise flexibility and inertiadriven 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 twodimensional 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 rigidbody drag scaling law of U^{2} to U^{4/3} (U is 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 threedimensional anisotropic flexibility is still not illustrated. In particular, comprehensive studies are required to understand the coupled fluid–structure interactions of biomembranes with morphological relevance, e.g. those with mechanical properties similar to skeletonreinforced fins or wings.
We herein present a numerical investigation by coupling fluid dynamics and structural mechanics to examine the propulsion performance of a skeletonreinforced 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 boundaryintegral equation representation for the surrounding fluid. To account for large threedimensional 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.
MATERIALS AND METHODS
Fin structure and kinematics
As a model of the real raystrengthened 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 xaxis, 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 nonuniform, 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 crosssections 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 springdampers between neighboring rays. Coupled with this beam–spring structure, we model dynamics of the surrounding fluid by using potential flow approximation.
We consider a threedimensional 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_{0}sin(ωt) following the motion of the caudal peduncle. The yaw motion with respect to the zaxis of Ray i isθ (t)=θ_{i}sin(ωt–ψ_{i}), i=1,..., 9. The Strouhal number, a nondimensional parameter characterizing the unsteady fluiddynamic effects, is defined to be S_{t}=2y_{0}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.
Mathematical formulations and numerical algorithms
We apply a threedimensional 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 dualcoordinate approach and the robust Euler parameters, this method achieves fullynonlinear 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 zerothickness 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 nearbody flow field as well as the hydrodynamic loading on the fin are then predicted using a boundaryelement approach based upon the boundaryintegral framework by segmenting the fin surface S_{b} into N_{b} panels S_{bj}, j=1,..., N_{b} (Appendix B).
By solving the potential flow around the body and determining the flow potential Φ, the hydrodynamic pressure p on the fin surface S_{b} is readily obtained through Bernoulli's equation. We have: (1) where ρ is the density of water (10^{3} kg m^{–3}). Integrating p over the fin surface, we obtain the hydrodynamic force as: (2)
The power required to drive the fin is given as: (3)
By definition, the thrust force, F_{t}, is the– x component of F. The propulsion efficiency, η, is defined as: (4) where F̄_{t} and P̄ represent the thrust force and power consumption averaged over one period T, respectively. The fin oscillation also creates a transverse force F_{r}, which is the y component of F. In the heterocercal mode, there also exists a nonzero lifting force, F_{z}, leading to a mean lifting coefficient F̄_{z}/½ρU^{2}c^{2}.
We thus characterize the performance of the fin by its mean thrust coefficient, F̄_{t}/½ρU^{2}c^{2}, the propulsion efficiency, η, the mean lifting coefficient, and the transverse force coefficient, F_{r}^{(1)}/½ρU^{2}c^{2}, where: (5)
As mentioned before, the structural function of the membrane in which the rays are embedded is simplified as linear springdampers 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.

At time t, we start from an initial guess, , of the configuration of the fin.

The hydrodynamic problem is solved based on .

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.

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.

If S is not sufficiently close to , we set and repeat Steps 2–4.
To eliminate the wellknown 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).
Convergence tests
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 leadingedge span s_{0}=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 crosssectional 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^{3} kg m^{–3}.
Kinematically, we set the forward speed to be U=0.1 ms^{–1}. The sway amplitude, y_{0}=0.5c. The Strouhal number, S_{t}=0.2. We also assume that the nine rays shear the same yaw motion, with θ_{0i}=θ_{0}=10^{0} andψ _{i}=90° (i=1,..., 9).
Table 1 shows the mean thrust coefficient, F̄_{t}/½ρU^{2}c^{2}, 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.
RESULTS
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 dorsoventral 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_{0}, 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.
Homocercal modes
In the homocercal motion, we assume that the yaw amplitudes of all the rays are the same, i.e. θ_{0i}=θ_{0} (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 (dorsoventral) bending. As t increases towards T/8 and T/4, accompanying the increasing sway, significant dorsoventral bending is observed, creating a bowshaped 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 laggedbehind undulations compared with the center part (Fig. 3). In Fig. 2, we also observe that the dorsoventral 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 bowshaped 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 `semifree'), 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 nearbody 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 bowshaped 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°<θ_{0}<40°. From Fig. 5 we see that the flexible fin generates a larger thrust when– 10°<θ_{0}<∼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,θ _{0}=–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 threedimensional structural flexibility. These characteristics are more clearly illustrated in Fig. 8, where crosssections of Figs 5, 6, 7 are plotted.
Heterocercal modes
In our simulations, the heterocercal mode, i.e. a mode with dorsoventral 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 (θ_{01}=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. θ_{01}=i/9×θ_{0}. 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°<θ_{0}<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 nearbody 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.
Visualization of threedimensional flow field
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 nearbody 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 midspan (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 twodimensional fluiddriven foil deformations (Zhu, 2007). The hind part of the fin, however, undergoes much smaller reorientation. This is caused by the threedimensional effect, in particular interactions with rays closer to the center via the membrane.
Fig. 13A,B shows the isosurfaces 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 twodimensional 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 nearbody 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 isosurface 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 dorsoventral 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 startup 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 ycomponent of the vorticity and the inplane velocity vectors, it is clear that each vortex pair (as a crosssection 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 counterrotating vortex pairs with outgoing jets induced between them.
DISCUSSION AND CONCLUSIONS
By using a fluid–structure interaction model combining a nonlinear Euler–Bernoulli beam model and a boundaryintegral 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 raystrengthened fins of many fishes. Structurally, the distribution of reinforcing rays creates a anisotropic flexibility. Two different locomotion modes, the homocercal (dorsoventral symmetric) motion and the heterocercal (dorsoventral 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 rayreinforced 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 dorsoventral 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 ringwithinaring 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 skeletonreinforced membranes. Compared with the conventional rigid foil design, these novel biomimetic propellers possess the following advantages.

Enhanced performance: as illustrated in our study, a flexible foil with threedimensional 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.

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.

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

Deployability: the skeletonreinforced 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 nearbody 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 threedimensional 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 biostructures in nature with similar structural characteristics, which can thus be categorized as skeletonreinforced 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 viscoelasticity of the bilayer and the skeleton, as well as their dynamic interconnectivity, 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 biostructures.
APPENDIX A. MODEL OF RAYS VIA NONLINEAR BEAM DYNAMICS
Our ray model is based upon a threedimensional nonlinear Euler–Bernoulli beam formulation reported by Tjavaras et al. (Tjavaras et al., 1998). This paradigm employs an Euler–Lagrangian dual coordinate system, including the abovedefined global coordinate system (x, y, z) and a local reference system with unit vectors τ, n and b in the tangential, normal and binormal directions of the ray, respectively. The two coordinate systems are related by: (A1) where i, j and k are unit vectors in the x, y and z directions, respectively. C is the orthogonal tensor representing the change of orientation from the global to the local reference frames. C is constructed by using the singularityfree method of Euler parameters based upon Euler's principal rotation theorem, which depicts an arbitrary reorientation as a single rotation by a principal angle about the principal unit vector of C, i.e. l. Correspondingly, four Euler parameters – β_{0},β _{1}, β_{2} and β_{3} – are defined as: (A2) The rotation matrix, C, thus becomes: (A3)
Dynamic equations of the ray movement are derived by considering the conservation of translational and angular momenta of an infinitesimal segment of the ray. We have: (A4) and (A5) where m is the mass per unit length, and s is the distance from this point to the basal end along the unstretched ray. V(s,t)=V_{τ}τ+V_{n}n+V_{b}b and ω(s,t)=ω_{τ}τ+ω_{n}n+ω_{b}b are the translational and angular velocities, respectively. T(s,t)=T_{τ}τ+T_{n}n+T_{b}b is the internal force inside the ray. M(s,t)=M_{τ}τ+M_{n}n+M_{b}b is the internal moment. ϵ(s,t) is the strain. Ω(s,t)=Ω_{τ}τ+Ω_{n}n+Ω_{b}b is the Darboux vector measuring the material torsion and curvatures of the ray. F_{c} is the constraining force from the membrane that controls the distance between neighboring rays. F_{h} is the hydrodynamic force to be defined later.
The compatibility relations, which state that the ray's configuration must be continuous in both space and time, should also be satisfied as: (A6)
The internal forces and moments are related to the strain ϵ and the Darboux vector Ω through the constitutive relations T_{τ}=EAϵ, M_{τ}=GJΩ_{τ}, M_{n}=EIΩ_{n}, M_{b}=EIΩ_{b}. By definition, A is the crosssectional area. EI and GJ are the bending and torsional stiffnesses, respectively. Our current simulations do not consider cases with ray twisting; therefore, the value of GJ is irrelevant. At this point, a hysteretic (material) damping can be incorporated by replacing the Young's modulus E with E+D∂/∂t, where D is the damping coefficient. In addition, the angular velocity ω is expressed in terms of the Euler parameters by: (A7)
The governing equations are closed by applying the spatial derivative of the Euler parameters. We have: (A8)
Finally, the system of equations is summarized into a vector form, i.e.: (A9) where Y=[ϵT_{n}T_{b}V_{τ}V_{n}V_{b}β_{0}β_{1}β_{2}β_{3}Ω_{τ}Ω_{n}Ω_{b}]^{T}. Detailed description of the matrices H and P is provided by Tjavaras et al. (Tjavaras et al., 1998).
Overall, 13 boundary conditions are required at the two ends of each ray. At the basal end, we apply: (1) β_{0}=0; (2) the velocities V_{τ}, V_{n} and V_{b}, as determined by the prescribed sway motion; (3) M_{τ}=0; and (4) the threedimensional 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) ; (2) the internal force T=0; and (3) the internal moments M_{n} and M_{b} disappear so thatΩ _{n}=Ω_{b}=0.
To solve Eqn A9 and determine the dynamics of a single ray, we numerically segment this ray into N_{r}–1, each with length Δs, by distributing N_{r} points s_{k} (k=1,..., N_{r}) along its unstretched length. A modified box method is then applied to integrate the system of equations from time step t_{i–1} to t_{i}=t_{i–1}+Δt. We have: (A10)
An iterative algorithm is employed to solve this system of implicit equations (Tjavaras et al., 1998). The method possesses secondorder accuracy in space. The accuracy of time integration, on the other hand, is first order. The purpose of selecting a lowerorder time integration scheme rather than the more conventional box method, which possesses secondorder accuracy, is to providing better numerical stability.
APPENDIX B. THE BOUNDARYINTEGRAL FORMULATION AND THE BOUNDARYELEMENT METHOD
The flow around the fin is assumed to be irrotational except for zerothickness 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 ϕ_{b} represents 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}.
The boundaryvalue problem for ϕ_{b} consists of the Laplace equation ^{2}ϕ_{b}=0 inside the fluid domain, the farfield condition ϕ_{b}→0 as x→∞, and the noflux condition on the fin surface S_{b}: (B1) where V_{b} is the instantaneous velocity at the surface of the fin, n is the unit normal vector pointing into the fin. The integral equation for ϕ_{b} is then obtained by invoking Green's theorem. At any point x on the fin surface we have: (B2)
Note that the term on the righthand side is known since the value of n·ϕ_{b} on S_{b} is given through Eqn B1. Solving the integral equation (B2), we obtain the unknown potential ϕ_{b} on S_{b}. The value of ϕ_{b} at any point x inside the fluid is then evaluated via Green's theorem, i.e.: (B3)
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 ϕ_{w} equals 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 selfinduced 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 nearbody flow field.
A loworder boundaryelement 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 loworder 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)}
 firstharmonic 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_{0}
 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_{0}
 sway amplitude
 β0, β1, β2, β3
 Euler parameters
 ϵ
 strain in a ray
 Δt
 time step
 η
 propulsion efficiency
 θ0i
 yaw amplitude of the ith ray
 (τ, n, b)
 tangential, normal and binormal 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 ith ray
 θ
 frequency
 ω=(ωτ, ωn, ωb)
 rotational velocity of a ray
 Ω=(Ωτ, Ωn, Ωb)
 Darboux vector measuring material torsion and curvatures of a ray
 © The Company of Biologists Limited 2008