SUMMARY
We investigate the thrust generation capacity of a thin foil consisting of a membrane strengthened by embedded rays that is geometrically, structurally and kinematically similar to pectoral fins of bony fishes during liftbased labriform locomotion. Our numerical model includes a fully nonlinear Euler–Bernoulli beam model of the skeleton and a boundaryelement model of the surrounding flow field. The fin undergoes a dorso–ventral flapping activated by rotations of the rays. Both the trailing edge vortices (TEV) and the leading edge vortices (LEV) are accounted for and modeled as shear layers. The thrust generation and propulsion efficiency are examined and documented. Our results show that synchronization of rays is pivotal to the performance of the system. A primary factor that determines the performance of the fin is phase lags between the rays, which create variations of the effective angle of attack at the leading edge as well as shape changes throughout the fin surface. Structural flexibility of the rays leads to passive deformations of the fin, which can increase the thrust generation and the propulsion efficiency.
INTRODUCTION
Fish locomotion is actuated through rhythmic undulation of flexible bodies and/or unsteady flapping of bodyattached fins. Morphologically, these fins fall into two categories: median fins (e.g. dorsal fins, ventral fins and caudal fins) and paired fins (e.g. pectoral fins and pelvic fins). Paired fins are employed mostly in motion stabilization, maneuvering and labriform swimming (Webb, 1973; Blake, 1979; Vogel, 1994; Standen, 2008).
In most bony fishes, the attached fins possess a characteristic skeletonreinforced membrane structure – a soft (and thin) collagenous membrane strengthened by bony fin rays. The Young's modulus of the fin rays is much larger than that of the membrane so that the bending stiffness of the fin is determined mostly by the embedded rays (Lauder et al., 2006; Lauder and Madden, 2007). The nonuniform distribution of these rays imparts anisotropic structural flexibility so that certain deformations may be favored while others are restricted. The basal end of each ray attaches to four separate muscles. This architecture enables a fish to control the motion of each ray individually. In addition, the curvature along a ray can also be actively controlled. According to morphological studies (Harder, 1975; 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. A hemitrich is connected with short ligaments and elastic fibers so that bending moments can be implemented along the length of a ray (Alben et al., 2007).
Recently, ray fins, with their capacity of sophisticated shape control, have attracted much attention owing to their possible applications on bioinspired engineering. The structural lightness and robustness, controllability, versatility and deployability of these biostructures make them promising prototypes of novel propulsion systems. Indeed, experimental tests have been undertaken to study hydrodynamic performance of a mechanical replica of a ray fin (Tangorra et al., 2007).
Our study is motivated not only by the scientific need to understand the structure vs function of ray fins but also by the biomimetic applications of these biostructures on autonomous underwater vehicles (AUVs) or microaerial vehicles (MAVs). The high efficiency and versatility of this design are especially useful under complex conditions, e.g. in confined space, low speed and unsteady currents. Our primary objective is to identify essential characteristics of ray fins that contribute most to their locomotion performance. From the biomimetic point of view, it is convenient to create simple, easily manufactured devices possessing essential characters of fish fins rather than exact copies of nature. For this reason, instead of exactly duplicating the fin of a specific species (Ramamurti et al., 2002; Mittal et al., 2006), in the current investigation we consider idealized fins with simplified geometry, internal structure and kinematics.
In a previous work we have computationally examined the dynamics of a trapezoid fin supported by a number of nonlinear beams resembling the caudal fins of bony fishes (Zhu and Shoele, 2008). Numerically, the fluid–structure interaction problem was solved by coupling a large Reynolds number boundaryelement method (BEM) with a nonlinear Euler–Bernoulli beam model. At the base ends, the orientations of the rays are controlled individually, a pivotal characteristic that allows for different locomotion modes to be achieved through synchronized ray motions. For example, by synchronizing the rays a caudal fin can undergo both homocercal (with dorso–ventral symmetry) and heterocercal (with dorso–ventral asymmetry) motions. In both cases, passive flexibility of the rays greatly increases the propulsion efficiency of the fin. Moreover, it also reduces lateral force generation and the sensitivity of propulsion efficiency to kinematic parameters. Both of these features are beneficial to the performance of the fin.
We herein propose a numerical study of a fin that is geometrically, structurally and kinematically similar to a pectoral fin during labriform swimming [specifically, as discussed later, the parameters chosen in this study are close to data measured from striped surfperch (Drucker and Jensen, 1996a; Drucker and Jensen, 1996b; Drucker and Jensen, 1997)].
At low speed, labriform swimming using pectoral fins is a primary mode of locomotion (Webb, 1973). Two extreme cases of fin movements, dragbased (i.e. rowing) and liftbased (i.e. flapping), have been identified (Blake, 1983; Vogel, 1994), although it is known that pectoral fin motions are usually complicated and rarely exhibit pure dragbased or liftbased movements (Gibb et al., 1994). In dragbased propulsion, it is found that with rigid fins positive thrust is generated in half of the period (i.e. no thrust generation during the recovery phase) so that the duty factor is 50% (Blake, 1981; Lauder and Jayne, 1996). Recently Mittal et al. (Mittal et al., 2006) and Tangorra et al. (Tangorra et al., 2008) have demonstrated that with structural flexibility the duty factor can be increased to 100%. In liftbased propulsion positive thrust generation is obtained during the whole flapping period and the duty factor may approach 100% (Vogel, 1994; Walker and Westneat, 2002). It was also found that liftbased propulsion generated higher efficiency, while dragbased propulsion generated higher thrust. Based upon this, it was suggested that dragbased mode was employed mostly in low speed maneuvering, while liftbased mode was used in powerconserving cruising (Kato, 1999; Archer et al., 1979). Incidentally, the liftbased mode is much simpler kinematically than the dragbased mode and thus easier to materialize in biomimetic devices.
Our focus is upon the performance of a pectoral fin in liftbased thrust generation (Westneat, 1996; Walker and Westneat, 1997). The major aspects we are going to test in this numerical study are: (1) effects of leading edge vortices (LEV) to the propulsion performance of the fin; (2) significance of the synchronization of ray motions; and (3) effects of the structural flexibility of the rays. Towards this end, we apply a fully coupled fluid–structure interaction model in which the rays are modeled as Euler–Bernoulli beams and the fluid motion is solved as potential flow. Vorticity generations from both the trailing edge and the leading edge are modeled as shear layers and mathematically represented as distributions of dipoles on thin sheets.
The paper is organized as follows. We first specify the geometry, internal structure and kinematics of the idealized fin. The fluid–structure interaction model, including the BEM (with LEV incorporated) and the fully nonlinear structural model, is then briefly described. Using this model, we investigate the dynamics and nearbody flow field of the flapping fin, focusing on the effects of fin kinematics as well as structural flexibility. Finally, conclusions are drawn.
MATERIALS AND METHODS
Physical model of a paired fin and its kinematics
As shown in Fig. 1, we consider a skeletonreinforced fin, which is reminiscent of the pectoral fins of a fish (e.g. Drucker and Jensen, 1997; Thorsen and Westneat, 2005). The fin is supported by 15 rays, tagged Ray 1 to Ray 15. A Cartesian coordinate system (x, y, z) is defined, within which the base ends of the rays form a baseline, which lies on the xaxis with length x_{0}. The trailing ends outline the trailing edge of the fin. The base ends of the rays are pinned during the motion. The fin has a uniform thickness d except for the vicinity of the trailing edge (within a region around 4% of the chord length) where tapering is implemented. The purpose of this tapering on the profile of the fin is to create a sharp trailing edge required for implementation of the Kutta condition as discussed later. It does not affect the rays (i.e. they are still considered to be uniform beams). θ_{i}, the angle between the base end of Ray i and the xaxis, varies evenly from 90 deg. (Ray 1) to 30 deg. (Ray 15) and does not change with time during the flapping motion. The length of Ray i, L_{i}, is given as L_{1}{1–0.8[(i–1)/14]^{5/4}}, where L_{1} is the length of Ray 1 (the ray at the leading edge).
The fin undergoes a translational motion in the x direction with constant speed U. To achieve liftbased thrust generation, a dorso–ventral flapping is activated by rotating each individual ray around the baseline. The angle between Ray i and the x–z plane is denoted as α_{i} and prescribed as , where ω is frequency, is the phase of the ith ray in rolling and t is time. According to this depiction, during the flapping motion each ray rotates (with the baseline as the rotating axis) between α_{i}=0 andα _{0} (α_{0} is maximum angle of rolling of each ray in one stroke). During the downstroke period a ray goes fromα _{1}=0 to α_{0}. The recovery fromα _{0} to 0 occurs in the upstroke period.
We also define a Strouhal number S_{t}≡ωL_{1}/2πU, where U is the forward speed. In this study we focus on a range of S_{t} between 0.1 and 0.5. To imitate observed fin activation at moderate swimming speed (Gibb et al., 1994), we assume that between fin strokes there are periods of inactivity during which the fin is held in its vertical position (α_{i}=0). Thus, in our simulations we consider only one period of flapping motion. The motion starts from the vertical position (Fig. 2). It is then followed by a downstroke and an upstroke (also known as abduction and adduction, respectively). The motion starts from Ray 1 and ends at the moment it returns to its original position (i.e. the duration of simulation is T=2π/ω). Unless otherwise specified, for simplicity we assume that the phase varies linearly from to ; thus, is used as a metric of the phase difference between neighboring rays. Both the structural deformation of the rays and the phase difference between them can cause the bow shape in lateral view as observed in experiments (Drucker and Jensen, 1997).
In the following calculations, we use L_{1}=0.1 m, x_{0}=0.03 m, d=2 mm, α_{0}=145 deg. and U=0.5 m s^{–1}. Although we aim at examining the mechanical principles about this type of design in general rather than the performance of a particular species, the choice of these parameters has certain biological relevance. Specifically, most parameters in this study are chosen following the experimental measurements by Drucker and Jensen of striped surfperch (Embiotoca lateralis) (Drucker and Jensen, 1996a; Drucker and Jensen, 1996b). In highspeed labriform cruising, this fish achieves a speed of 0.48 m s^{–1}, around 80% of pectoral caudal gait transition speed. As shown in their experiments, the flapping amplitude of the fin tip is about twice the span of the fin, and the Strouhal number based upon the flapping period reported by Drucker and Jensen (Drucker and Jensen, 1997) is around 0.3–0.4. These are within the range of the fin kinematics considered in our present study.
Structural dynamics model
Despite the fact that rays have complicated internal structures and are usually tapered/branched near the trailing ends, for simplicity in our model the rays are depicted as cylindrical Euler–Bernoulli beams with uniform crosssections. These beams can be stretched (with tensile stiffness EA, where E is the Young's modulus and A is the crosssectional area) or bent (with bending stiffness EI). In this study we concentrate on the passive deformability of the rays, while the active deformability activated by the muscles/tendons connected to the rays (Alben et al., 2007) is not considered.
Numerically, a threedimensional nonlinear model is employed to simulate the dynamics of the rays (Zhu and Shoele, 2008). This method was originally developed to study the fully nonlinear dynamics of marine cables (Tjavaras et al., 1998), and is characterized by its capacity to simulate arbitrary configurations and large deformations owing to the implementation of a dual Euler–Lagrangian coordinate system and the robust Euler parameters (e.g. Junkins and Turner, 1986) instead of the more conventional Euler angles, which suffer from singularities.
To elaborate, by invoking the conservation of momentum we have (Tjavars, 1996): (1) where m is the mass per unit length, s is a Lagrangian marker measuring the distance from a point on the ray to the basal end along the unstretched ray, ϵ is the strain, V and ω are the translational and angular velocities, respectively, T is the internal force inside the ray, F_{c} is the constraining force from the membrane that controls the distance between neighboring rays and F_{h} is the hydrodynamic force on each ray.
Similarly, from conservation of angular momentum another set of equations is obtained as: (2) where M is the internal moment. Ω is the Darboux vector measuring the material torsion and curvatures of the ray. τ is the tangential vector along the ray. M and Ω are related by M_{τ}=GJΩ_{τ}, M_{n}=EIΩ_{n} and M_{b}=EIΩ_{b}. The subscripts τ, n and b represent components in the tangential, normal and binormal directions, respectively, and GJ represents the torsional stiffness. In the following simulations with flexible rays, unless otherwise specified we assume E=1 GPa (e.g. Lauder et al., 2006); a much larger E (100 GPa) is applied to approximate the behavior of a rigid ray. Because no twisting motion is considered, the value of GJ is irrelevant. To model the viscoelastic behavior of the rays, we further replace E by E+D∂/∂t, where the viscous coefficient D is chosen to be 0.4 sGPa. The diameters of the rays are estimated by using the argument in Alben et al. (Alben et al., 2007), i.e. during normal swimming the bending energy stored in the fin must be balanced by the kinetic energy of the flow. According to their estimation, for this to be true the EI of a ray, its length L and the swimming speed U must be related by EI/L<U^{2}L^{3}/12, based upon which we choose the radius of a ray to be 1 mm.
For the ith ray, its base end is pinned in space with the orientation determined by the angles θ_{i} andα _{i}. Its trailing end is treated as a free end with no external load. The resulting system of nonlinear equations is solved by an efficient finite difference algorithm, the modified box method (Tjavaras, 1996).
The hydrodynamic load on the rays is imposed via the attached membrane, which is modeled as arrays of linear springs connecting the neighboring rays. According to Lauder et al. (Lauder et al., 2006), the membrane between fin rays has a modulus of about 0.3–1.0 MPa. If its thickness is smaller but comparable with the effective diameter of the ray, the spring constant per unit length should be O(10^{3})Nm^{–2}; in this study we choose the spring constant to be 2000 N m^{–2}. The rays are assumed to have the same density as water. The mass of the membrane is not considered because it is much smaller than the added mass.
Fluid dynamics model
We consider cases with large Reynolds numbers (Re) so that the flow field can be described by defining a flow potentialΦ (x,t), where x=(x, y, z) is the position vector [as demonstrated by Anderson et al. (Anderson et al., 1998)], predictions from a potential flow model are in reasonable agreement with experimental measurements for a flapping foil when Re∼O(10^{4}). The vorticity shed from the trailing edge of the fin is physically modeled as a zerothickness shear layer and mathematically represented by a thin sheet of distributed dipoles. The strength of the newly shed vortices is determined by the Kutta condition at the trailing edge (Zhu et al., 2002).
A key characteristic of our current fluid dynamics model is the inclusion of vorticity shed from the leading edge area. Following the treatment of leading edge separations within the potential flow framework (Katz, 1980; Dong, 2007), we model the leading edge separation as thin shear layers (i.e. in the same fashion as the trailing edge separation). This depiction, per se, is consistent with experimental visualizations that show that vorticity shed from the leading edge area starts as thin vorticity layers before they rollup and form individual vortices (Taneda, 1977).
Katz studied the dynamics of twodimensional flow around a thin plate (Katz, 1980). In that model, LEV was depicted as a shear layer initiated from a prescribed separation point near the leading edge. Good comparisons with experiments were achieved over a large range of angle of attack. A similar model was applied to solve the flow around a thin plate, in which vorticity shed from both edges are represented as vortex sheets emanating from the sharp edges (Jones, 2003). Dong extended this method to three dimensions, unsteady flow and finite foil thickness (Dong, 2007). An algorithm was developed to predict the unsteady separation lines by using the adverse pressure gradient (Faltinsen and Pettersen, 1987). The accuracy of this method has been tested via comparisons with experimental measurements of a foil undergoing pitching and heaving motions (Read, 2000). It was demonstrated that by including this simplified LEV representation the accuracy of the potential flow based prediction is significantly improved. Specifically, it was found that in general the incorporation of LEVs greatly reduce the predicted thrust generation.
Taking advantage of the special geometry of the fin we study, in our model we assume that the locations of leading edge separation are fixed at sharp edges existing near the leading edge, so that although we model LEV in the same manner as the abovementioned studies, the complicated algorithm to locate separation lines is not required. Furthermore, by using two vortex sheets shed alternatively, one from the upper side and the other from the lower side of the leading edge, we assume that each vortex sheet remains on its side of the fin, so that the singularity associated with a vortex sheet passing through the fin is avoided. Otherwise to eliminate this singularity it requests special treatment (Dong, 2007).
As shown in Fig. 3, in our approach the LEV is modeled as two shear layers, one from the upper side and the other from the lower side of the leading edge, where sharp angles exist between neighboring panels. Similar to their trailing edge counterpart, the strengths of these shear layers are determined by invoking the Kutta condition. In addition, to prevent interceptions of these shear layers with the fin itself, which will induce singularity, at each time step the distances between the shear layers and the fin are carefully monitored. Whenever interception is likely to occur, the shear layer is pushed away from the fin so that singularity is prevented. Although our model allows continuous vorticity generation from both separation lines at the leading edge, it has been demonstrated in our results that at any given moment only one of them is associated with significant vorticity shedding.
Fluid–structure coupling
To solve the fluid and the structural problems simultaneously in a fully coupled algorithm, we employ the following iteration method: at each time step to start the iteration we apply an initial guess of fin deformation and solve the hydrodynamic problem through the BEM. Based upon this the hydrodynamic load is updated and applied to calculate the body deformation; if this recalculated deformation is not sufficiently close to the initial guess it is used as the initial guess in the next iteration step.
Detailed description and validation of this iteration method are provided in previous publications (Zhu, 2007; Zhu and Shoele, 2008).
Numerical issues
Numerically, the fluid dynamics problem is solved using a boundaryelement approach by discretizing the fin surface S_{b} into N_{b} panels S_{bj}, j=1,...., N_{b}, and the wake surfaces S_{w} into N_{w} panels, where N_{b} is the total number of panels on the fin surface, S_{bj} represents one panel and N_{w} is the number of panels on the wake. By satisfying the no flux condition at the centroids of the panels on the fin surface this boundary value problem is formulated as a system of linear equations and solved numerically (Katz and Plotkin, 1991; Zhu et al., 2002).
After obtaining Φ(x,t), the hydrodynamic pressure (p) is determined through Bernoulli's equation, and given as: (3) where ρ is the density of water (10^{3} kg m^{–3}). We assume that the hydrodynamic force acting on this panel, pΔS_{bj} (ΔS_{bj} 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).
Hydrodynamic forces
Integrating p over the fin surface, we obtain the overall hydrodynamic force on the fin (F) as: (4) where n is a unit normal vector pointing into the fin. The power (P) required to drive the fin is given as: (5) where V_{b} is the instantaneous velocity at the surface of the fin. The propulsion efficiency η is defined as: (6) where F̄_{T}=–F̄_{x} and P̄ represent the thrust force and power consumption averaged over one period T, respectively. To characterize the hydrodynamic performance of the fin, we herein define a mean thrust coefficient .
The dynamics of the fin is also characterized by generation of lift force F_{z} and lateral force F_{y}. Correspondingly, we will examine the lift coefficient and the lateral force coefficient , focusing upon the firstharmonic components and .
RESULTS
Using the fully coupled model we study cases with S_{t}=0.1∼0.5 and varying from 30 deg. (π/6) to 120 deg. (2π/3). To put this choice of phase lag in perspective, we note that, according to Webb (Webb, 1973), the phase lag between leading and trailing edges varies with swimming speed, ranging fromπ at 0.2 BL s^{–1} to 0.2π at higher velocities.
Flow visualization and fin deformation
We first study the nearbody flow field around the fin as well as the deformation of the fin itself during a single downstroke–upstroke period.
In Fig. 4 we plot the inplane streamlines within a plane parallel to the rotating axis (i.e. the baseline). This plane rotates with Ray 1 and its distance to the baseline remains as 0.75L_{1}. The fin is flexible with S_{t}=0.40 and deg. Four snapshots at different phases of motion are shown. At t=T/8 (i.e. at the beginning of the downstroke motion), a pair of counterrotating circulations are generated, one from the leading edge (marked as C1) and the other from the trailing edge (marked as C2). These circulations are then shed into the wake (see the snapshot at t=T/2). After the transition from downstroke to upstroke motions, a new pair of circulations (C3 and C4) appears (t=3T/4). Thus, at the end of the upstroke period (t=T), two pairs of vortices are left in the wake. This sequence of vortex shedding is qualitatively consistent with the depiction by Lauder and Drucker (Lauder and Drucker, 2004).
Fig. 5 displays structural deformations of the fin during a flapping period. For comparison, in the same figure we also draw the shapes of a fin with rigid rays. It is seen that effects of ray flexibility is most pronounced near the leading edge. This is attributed to the several facts. First, the rays in that vicinity are longer and thus more deformable than others. Second, the ray at the leading edge (Ray 1) is loaded from one side only so that its deformation is different from neighboring rays, which sustain hydrodynamic loads from both sides.
The deformability of the rays has profound effects upon the interaction between the fin and the surrounding flow field. Our simulations demonstrate that it may significantly change the effective angle of attack at the leading edge. For example, in Fig. 6 we plot the inplane streamlines of the flow field measured in a reference system fixed on the leading edge of the fin (i.e. this reference system follows both the forward and the lateral motions of the intersection point between this plane and the leading edge). During most of the cycle (except for the snapshot at t=T/2, when the transition from downstroke to upstroke occurs), the leading edge of the flexible fin is better aligned with the incoming flow than that of the rigid fin. This effect is more clearly shown by plotting the magnitude of the effective angle of attack at the tip of the fin (i.e. the angle between the leading edge and the incoming flow vector measured in the moving coordinate system) (Fig. 7). As illustrated in previous experiments using flapping foils (e.g. Anderson et al., 1998), such a reduction in effective angle of attack may mitigate generation of LEV.
Fig. 8 shows the threedimensional flow field visualized via dipole distribution within the vorticity layer as well as isosurfaces of vorticity behind a fin when only the TEV are accounted. It is seen that concentrations of dipoles, displayed in Fig. 8A as red and blue areas, correspond to threedimensional vortex rings in Fig. 8B. Within one flapping period, two connected vortex rings, one during upstroke and the other during downstroke, are generated, forming a figure of `8' shape in the wake. As shown in Fig. 9A,B, with LEV included during the downstroke period the LEV is generated primarily at the upper side, while in the upstroke period the vorticity sheds from the lower side. In Fig. 9C we plot the isosurfaces of the vorticity field created by the LEV only (i.e. influences of the TEV and the body itself are not considered). Owing to the complexity of the vorticity field, it is difficult to observe clearly defined vortex rings in the wake.
To summarize, during the downstroke motion two concentrations of vorticity, one from the trailing edge of the fin and the other from the upper side of the leading edge, are generated and shed into the wake. Two additional vorticity concentrations are shed during the upstroke period from the trailing edge and the down side of the leading edge, respectively.
Thrust generation and propulsion efficiency of the fin
In Fig. 10A,B we plot the mean C_{T} and the (defined as U× the mean thrust divided by the mean power expenditure) as functions of the S_{t} and , the phase lag of Ray 15 (the phase varies linearly from 0 when i=1 to when i=15). No LEV is included in the simulation and the rays are assumed to be rigid. As shown in the figures, for fixed , C_{T} increases with S_{t}, while η reaches its maximum value at an optimal S_{t}, which depends upon . However, when S_{t} is fixed, positive thrust and η are achieved only when is below a threshold value, which increases with S_{t}.
The inclusion of LEV (even though the rays are still rigid) significantly changes C_{T} (Fig. 10C) as well as η (Fig. 10D). Most notably, the mean thrust is greatly reduced. Within the range of parameters considered in this study (i.e. 0.1<S_{t}<0.5, 30 deg.< deg.), as predicted by the model without LEV the maximum C_{T} is around 1.3 whereas with LEV this value is merely 0.35. In both cases, however, predictions of the maximum values of η are around 0.75. The effect of LEV is most pronounced in small values of . Our simulations show that, with LEV included, positive thrust and efficiency can be obtained only when is sufficiently large (>∼30 deg.). Large values of also cause negative thrust. Indeed, large corresponds to significant pitching of the fin, which may be applied as a braking mechanism (Drucker and Lauder, 2003). Besides, at large when the thrust generation is close to zero, the discrepancy between predictions with and without LEV becomes insignificant. This is attributed to the fact that in these cases the fin motion resembles the feathering behavior so that the effective angle of attack at the leading edge is minimized and LEV is thus mitigated.
Hereafter, we include LEV in all the testing cases.
The effects of fin flexibility upon thrust generation and η are displayed in Fig. 10E,F. By comparing the mean C_{T} and the η in flexible cases with rigid cases, it is seen that flexibility has limited effect in low S_{t}. However, in higher S_{t} it is able to improve the performance of the flapping fin. To elaborate, the threedimensional anisotropic flexibility imparted by flexible rays can significantly increase the mean C_{T} (by 40% in some cases), accompanied by a slight increase in η (by about 10%).
Fig. 11A shows the time history of the C_{T} over one flapping period at S_{t}=0.4 and =40 deg. It is seen that although the thrust force remains positive over most of the period, there exists a short time interval with negative thrust during the transition between the downstroke and the upstroke motions so that the corresponding duty factor is clearly less than 100%. This is qualitatively consistent with the experimental study of the liftbased swimmer bird wrasse (Gomphosus varius) as reported by Walker and Westneat (Walker and Westneat, 2002). Employing the experimentally recorded fin kinematics of this fish (Walker and Westneat, 1997), Ramamurti et al. (Ramamurti et al., 2002) has carried out fullyviscous simulations using an unstructured finiteelement method, yielding a time history of thrust similar to our prediction (see fig. 6 in that report). Based upon descriptions in their studies, we estimate that the fin size is around 2.5–3 cm so that peak values of the C_{T} obtained in their simulations are around 0.5–0.7, which are also comparable with our results. In Fig. 11B we study a case with =60 deg. and find that the period with negative thrust is greatly reduced. In both cases (=40 deg. and =60 deg.), the structural flexibility increases the thrust and its effect is most pronounced during the upstroke period, which is consistent with the fin deformation as shown in Fig. 5.
As suggested in Fig. 11, the duty factor clearly depends upon the kinematic parameters as well as the fin flexibility. To further illustrate this, in Fig. 12 we plot the duty factor as a function of S_{t} and . The effect of flexibility to increase the duty factor is herein quantitatively confirmed. Within this range of parameters, our simulations record a maximum duty factor of 87%.
In Fig. 13 we plot the firstharmonic components of the lift coefficient and the lateral force coefficient as functions of the S_{t} and the phase lag . Similar to the case in thrusts generation, the influence of structural flexibility is also more pronounced in high S_{t}. Its specific effect depends upon the range of kinematic parameters under consideration. For example, at S_{t}=0.4 this flexibility decreases the lateral force when <45 deg. and increases it otherwise; the lift force is increased when 45 deg.<<105 deg.
We finally check the effect of different levels of ray flexibility, achieved by varying E while keeping the ratio D/E and other parameters unchanged. As shown in Fig. 14, peaks of both C_{T} and η occur when E is around 0.6 GPa. For sufficiently soft fins significant drop of thrust generation is observed, echoing findings from flexible foils (e.g. Katz and Weihs, 1978; Zhu, 2007).
DISCUSSION
By using a fully coupled model, we have numerically investigated the fluid–structure interaction mechanisms involved in the liftbased propulsion of an idealized pectoral fin. The fin is strengthened from inside by a skeleton of fin rays, which are modeled as Euler–Bernoulli beams. The fin motion is activated by rotations of each individual ray. A `winglike' dorso–ventral motion characteristic of the liftbased mechanism is achieved by synchronized ray motions.
Our results indicate that the hydrodynamic performance of such a system is determined by its frequency of motion (represented by the Strouhal number), phase difference among the rays, as well as structural properties of the rays. By varying the phases of the rays, the fin achieves a combination of roll and pitch motions. Based on this mechanism the effective angle of attack at the leading edge as well as the shape of the fin can be adjusted so that the thrust generation and the propulsion efficiency can be optimized. With properly chosen kinematic parameters, positive thrust is observed in most of the downstroke–upstroke period. High propulsion efficiency (with a peak value of ∼0.8 within the range of parameters considered in this study) is observed.
We studied the effect of passive deformation of the fin due to structural flexibility of the rays. Our results show that such passive deformations indeed increase the thrust generation as well as the propulsion efficiency [although the increase in efficiency is not as significant as in the case of caudal fins (Zhu and Shoele, 2008)]. By visualizing the nearbody flow field we find that a primary effect of structural flexibility is the reduction of the effective angle of attack, which suggests a mechanism for the change in performance.
Further investigation shows that performance of the system depends not only upon , the maximum phase lag between all the rays, but also upon the detailed distribution of the phases. For example, in Fig. 15 we plot the thrust coefficient and propulsion efficiency of the fin when are quadratic (rather than linear) functions of i. It is clear that the performance of the fin is considerably changed. This not only demonstrates the subtlety of liftbased propulsion but also brings in more parameters to be optimized for the design of future biomimetic devices.
As mentioned in the introduction, in dragbased swimming (rowing mode) with flexible fins, the duty factor of swimming reaches 1 (Mittal et al., 2006; Tangorra et al., 2008). However, the studies by Ramamurti et al. (Ramamurti et al., 2002) as well as our current simulation of liftbased swimmers (in flapping mode) predict duty factors smaller than 1. This may be attributed to the following factors: (1) these studies were carried out at different forward speeds. As indicated by Tangorra et al. (Tangorra et al., 2008), as forward speed increases (with other parameters fixed) the capacity of thrust generation decreases and this may eventually generated a lower than 1 duty factor. (2) It appears that structural flexibility has a larger effect on the performance of a rowing fin. Indeed in the rowing mode the duty factor of 1 is only achieved with a highly flexible fin. To quantitatively understand this effect, future investigations using fins with different levels of flexibility are required. (3) The sequences of vortex shedding (dorsal/ventral edge vortices in flapping mode vs dorsal/ventral edge vortices along with abduction/adduction tip vortices in rowing mode) are different between flapping and rowing modes. These differences may affect the thrust drop during the transition from downstroke (instroke) to upstroke (outstroke) motions and consequently the duty factor.
One of the key characteristics of the current study is the incorporation of a LEV model. As suggested by our results, two concentrations of vorticity, one generated during the downstroke motion and the other during the upstroke motion, are shed from the leading edge into the wake. LEV significantly affects the thrust generation and propulsion efficiency of the fin. The current model of LEV stems from the potential flow framework. At a result, although this method provides a convenient alternative to the Navier–Stokes equations, it is based upon simplified descriptions of the otherwise complicated mechanisms involved in vorticity shedding, vortex formation and vortex–body interactions. For example, as indicated in our recent Navier–Stokes study (Zhu and Peng, 2009), when the vortex shedding and the body motion are properly synchronized, the energy of LEV can be effectively recovered so that the performance of the system is enhanced. For accurate prediction and comprehensive understanding of mechanisms like this and their effects upon the dynamical performance of the fin, more sophisticated solvers (e.g. those based on Navier–Stokes equations) are necessary. Although these models are expected to be more accurate, they are computationally expensive so that it is critical to carefully choose the cases to be tested. The highly efficient potentialflowbased simulations discussed in this study can provide valuable guidance.
LIST OF ABBREVIATIONS
 A
 crosssectional area of a ray
 BEM
 boundaryelement method
 C_{T}(C̄_{T})
 (mean) thrust coefficient
 C_{L}[C^{(1)}_{L}]
 (firstharmonic) lift coefficient
 C_{y}[C^{(1)}_{y}]
 (firstharmonic) lateral force coefficient
 d
 thickness of the fin
 D
 material damping coefficient of the rays
 E
 Young's modulus of the rays
 EA
 tensile stiffness
 EI
 bending stiffness
 F
 hydrodynamic force on the fin
 F_{c}
 constraining force from the membrane that controls the distance between the neighboring rays
 F_{h}
 hydrodynamic force on a ray
 F_{T}(F̄_{T})
 (mean) thrust force
 F_{y}
 lateral force
 F_{z}
 lift force
 GJ
 torsional stiffness
 L
 length of ray
 LEV
 leading edge vortices
 L_{i}
 length of Ray i
 L_{1}
 length of Ray 1 (the leading edge ray)
 M
 mass per unit length
 M
 internal moment
 n
 unit normal vector pointing into the fin
 N_{b}
 number of panels on the fin surface
 N_{w}
 number of panels on the wake
 p
 hydrodynamic pressure
 P(P̄)
 (mean) power spent
 Re
 Reynolds number
 s
 distance from a point on a ray to the basal end measured along the unstretched ray
 S_{b}
 fin surface
 S_{t}
 Strouhal number
 S_{w}
 wake surface
 U
 forward speed
 t
 time
 T
 period (2π/ω)
 T
 internal force inside a ray
 V
 translational velocity of a ray
 V_{b}
 instantaneous velocity at the surface of the fin
 x_{0}
 length of baseline
 (x, y, z)
 global coordinate system
 α0
 maximum angle of rolling of each ray in one stroke
 αI
 angle between base end of Ray i and the x–z plane
 ϵ
 strain in a ray
 η
 propulsion efficiency
 θ_{i}
 angle between base end of Ray i and the xaxis
 ρ
 density of water
 (τ, n, b)
 tangential, normal and binormal unit vectors in the Lagragian coordinate system
 phase of the ith ray in rolling
 Φ
 flow potential
 ω
 frequency
 ω
 angular velocity of a ray
 Ω
 Darboux vector measuring material torsion and curvatures of a ray
FOOTNOTES

Computational supports from TeraGrid resources provided by the San Diego Supercomputer Center (SDSC) and the National Center for Supercomputing Applications (NCSA) are acknowledged.