Hydrodynamic characteristics of a lunate shape oscillating propulsor

Ocean Engineering 26 (1999) 519–529

Hydrodynamic characteristics of a lunate shape oscillating propulsor Pengfei Liua, Neil Boseb,* a

MMC Engineering and Research, P.O. Box 314, Burin Bay Arm, Newfoundland, Canada, AOE 1GO b Ocean Engineering Research Centre, Memorial University of Newfoundland St. John’s, Newfoundland, Canada, A1B 3X5 Received 18 August 1997; accepted 13 November 1997

Abstract Research was conducted to study the hydrodynamic efficiency of a foil with aft-swept wing tips. A potential flow based time domain panel method was formulated to predict the performance of a lunate and rectangular foil in large amplitude, unsteady motion. Skin drag was approximated and boundary layer growth and separation were also estimated. Hydrodynamic efficiency was evaluated in terms of propulsive efficiency and thrust coefficient of the foil. Results are presented for a lunate shaped planform and for a rectangular foil. Predictions show that the lunate shaped planform has a substantially higher propulsive efficiency (13% higher) than the rectangular foil under heavy load conditions when the feathering parameter is zero, throughout a range of reduced frequencies (0.2 to 1.8). Under a medium load condition, however, the rectangular foil gave a higher propulsive efficiency at reduced frequencies less than 0.5 and the same efficiency value at a reduced frequency of 1.8. For a practical range of reduced frequencies between 0.5 and 1.0, the lunate tail gave higher propulsive efficiency. The lunate planform gave a lower thrust coefficient at a heavy load and higher thrust at a medium load condition than the rectangular planform for all reduced frequencies.  1998 Elsevier Science Ltd. All rights reserved. Keywords: Hydrodynamic efficiency; Panel method; Time domain; Lunate planform

* Corresponding author. 0029-8018/99/$—see front matter  1998 Elsevier Science Ltd. All rights reserved. PII: S 0 0 2 9 - 8 0 1 8 ( 9 8 ) 0 0 0 1 8 - 3

520

P. Liu, N. Bose / Ocean Engineering 26 (1999) 519–529

1. Introduction Lunate-shaped wings and flukes are found in many fast moving animals such as birds and marine mammals. Numerous previous studies have been done on the propulsive performance of oscillating foils. Work on the propulsive performance of lunate planforms was presented by Newman and Wu (1974). Using a 3-D unsteady lifting surface theory, Chopra and Kambe (1977) studied a set of five planforms. A lunate-tail planform was reported to give the best overall propulsive performance (Chopra and Kambe, 1977). Aerodynamic efficiency of a number of wings with swept tips were studied (Van Dam, 1985) using a steady panel method VSAERO. Higher efficiency was reported for swept-tip planforms at both high and low angles of attack (Van Dam, 1985). The characteristics of tip-swept planforms in steady flow were also used to analyze the propulsive performance of marine swimmers (Van Dam, 1986). A parametric analysis for swept wings with sheared tips in steady flow has also been presented (Fremaux et al., 1990). Study of naturally occurring planforms for propulsive performance includes work by Bose and Lien (1989); Liu and Bose (1993). Recently an unsteady panel method has been used to calculate the propulsive performance of oscillating wings with nonuniform spanwise heave amplitude (spanwise flexibility or flapping oscillating wings): examples are Liu and Bose (1996); Smith (1996); Vest and Katz (1996); Liu and Bose (1997). In the present study, a time domain panel code was used to predict the propulsive performance for a rectangular and a lunate planform. Results and comparisons are presented for the propulsive efficiency and thrust coefficient with regard to the reduced frequencies and feathering parameters. The question was, when large amplitude theory is applied to the problem of a non-zero thickness foil, which gives control of the instantaneous hydrodynamic angle of attack to avoid a non-practical prediction due to severe separation and stall, how well would the lunate-shaped foil outperform the rectangular foil?

2. Motion parameters and fluid dynamics model Many motion parameters of an oscillating foil originated from studies of unsteady aerodynamics, namely, wing flutter and response to gusts (e.g. Theodorsen, 1935; Von Ka´rma´n and Sears, 1938). Major motion parameters for an oscillating foil are the heave amplitude h, pitch amplitude ␣, oscillating frequency ␻, foil forward speed V, and phase angle ⌽ between pitch and heave. The thrust T and propulsive efficiency ␩, the ratio of output thrust power to input power are a function of foil loading. Another loading related motion parameter, reduced frequency, k ⫽ ␻Cr/V, where Cr is the rootchord, was often used as the abscissa to present results of thrust T or thrust coefficient Ct and propulsive efficiency ␩. The larger the value of k, the higher the load on a foil. For wings with varied chord length, the reduced frequency k varies across span and an advance ratio

P. Liu, N. Bose / Ocean Engineering 26 (1999) 519–529

J⫽

521

␲V ␲ ␪␲ ⫽ , ⫽ ␻h kh0 ␣

was introduced (Bose and Lien, 1989). This advance ratio J provides a reasonable basis for comparison for planforms with nonuniform chord length. A feathering parameter, ␪ ⫽ ␣V/␻h, is also used in evaluation studies of the performance of an oscillating foil. This feathering parameter is the ratio of the pitch amplitude ␣ to the angular amplitude of the kinematic velocity of the foil,

␣␯ ⫽ tan−1

冉 冊

␻h ␻h ⬇ , V V

under a small amplitude assumption. The hydrodynamic angle of attack at a feathering parameter ␪ ⫽ 1 is zero if the phase angle between pitch and heave (pitch leading heave) ⌽ ⫽ 90°. At ␪ ⫽ 0, the hydrodynamic angle of attack is maximum and hence the foil has the highest loading at a given oscillating frequency, ␻. The small amplitude assumption worked well in the past for wings in flutter and gust applications. For oscillating foils, whose primary function is to produce sufficient thrust, a large amplitude of heave is often required. In large amplitude motion, the term ␻h/V is no longer small, especially when the oscillating frequency ␻ is big. Therefore, a large amplitude feathering parameter ⌰ was used (Liu, 1996). This large amplitude feathering parameter has a time dependent value, expressed as

␪instant ⫽

␣ sin ␻t . ␻h sin ␻t tan−1 V

冉

冊

An instantaneous hydrodynamic angle of attack ␣instant ⫽ ␣␯ ⫺ ␣, was also used. The value of the instantaneous angle of attack of a foil is negative when ⌰instant is greater than 1, in which case the foil produces negative thrust. A very large pitch angle will produce a large ⌰ with a value greater than one (⌰ > 1). The fluid flow modelling is based on the classic low order steady panel method (Katz and Plotkin, 1991; Maskew, 1987). This low order boundary element method was extended to a time domain for unsteady flow and to a case in which a flexibility in both chordwise and spanwise directions was present. A detailed formulation of rigid foils was described by Katz and Plotkin (1991); a surface and wake panel arrangement and doublet collocation algorithm were illustrated clearly by Maskew (1987). A brief description on the formulation of the method and its solution are given below. Panel methods are also referred to as boundary integral methods or boundary element methods. They were initiated in the early 1960s for aerospace applications. Perhaps the best known example, the Douglas Neumann program, was developed by Hess and Smith (Hess and Smith, 1964, 1967). For the numerical model in the present study, the governing equation is formed from Green’s third identity (Katz and Plotkin, 1991):

522

P. Liu, N. Bose / Ocean Engineering 26 (1999) 519–529

1 兰 4␲

冕

␴

foil

冕

1 1 兰 dS ⫹ r 4␲

foil

␮

冉冊

1 ∂ 1 dS ⫹ 兰 ∂n r 4␲

冕

␮

wake

冉冊

∂ 1 dS ⫽ 0, ∂n r

(1)

with an appropriate kinematic boundary condition on the surface of the flexible foil at each time step which accounts for both the sinusoidal motion and the flexibility of the foil. The boundary conditions are: →

␴⫽→ n ·V ⬁

(2)

where → n is the normal vector of the foil’s surface, and therefore →

→ ⵜ (⌽ ⫹ ⌽⬁)·n ⫽0

(3)

on the body boundary. The influence coefficient matrices for Eq. (1) were obtained using Newman’s algorithms (Newman, 1986). Morino’s steady Kutta condition was applied at the trailing edge with the consideration of the effect of a non-zero thickness trailing edge on the prediction of the velocity profile on the surface (see Maskew, 1987). From the combination of the doublet and source coefficients, along with the boundary conditions and the Kutta condition, a set of linear equations was obtained and solved for each time step of the motion of the foil. In many computational jobs, coefficient matrices that are produced from integral-L.P.D.E (linear partial differential equations) are symmetrical and sparse. The coefficient matrices containing the coefficients of the equation obtained in this study, however, are dense, asymmetrical matrices. In a computer program OSFBEM (oscillating foil boundary element method) for the present work, the Bi-CGSTAB (BiConjugate gradient stabilized) method was used (Freund et al., 1991) to invert these matrices. The Bi-CGSTAB method was not only efficient but also stable for the linear equation systems tested. The matrix solver calculates the doublet velocity potential strength at each panel at each time step. A quadratic curve fit for the velocity potential in time history was done and this was used to calculate the derivatives of the velocity potential with respect to time, which is the term ∂⌽/∂t in the unsteady Bernoulli equation that is used to evaluate the hydrodynamic forces. The pressure distribution around each section was obtained by using the unsteady Bernoulli equation. The pressure distribution was then used to find the lift coefficient cl and thrust coefficient ct of the oscillating foil at each time step. A 2-D boundary layer method (e.g., Moran, 1984) was used to estimate friction coefficients for this 3-D foil in unsteady motion. Each strip section on the 3-D planform was assumed to be a 2-D section (for details, see Liu, 1996). The boundary layer growth was assumed to start at a stagnation point at the leading edge. In the laminar flow region, cf, the coefficient of skin friction was calculated using the Thwaites method; at the transition point, Michel’s method was applied. In turbulent flow, Cf was obtained by Head’s method (Moran, 1984). At each time step, the total skin friction coefficient Cf was calculated as Cf ⫽

1 Afoil

冘 N

(cf)n(Apanel)n,

n⫽1

(4)

P. Liu, N. Bose / Ocean Engineering 26 (1999) 519–529

523

where N was the total number of panels, (cf)n was the local skin friction coefficient at each panel and Apanel was the surface area. This value of Cf was then used to modify the total thrust and input power calculated at each time step, similarly to that described in Liu and Bose (1993). 3. Results and discussions A convergence study was done to ensure stability of the results. This study covered the number of time steps, number of panels, panel spacing (cosine, uniform and log, spanwise and chordwise, a total of nine combinations), the size of the time steps and the effect of all the above on the stability and computing efficiency of the results. Detailed discussions and plots were presented by Liu (1996). Results showed that a good convergence rate and computing efficiency was obtained when the spanwise panel spacing was taken as uniform and the chordwise one was taken as cosine. While a smaller time step size gave a more accurate result, about 300 steps per period of the cycle was found to be acceptable. To examine the effect of a swept planform on the propulsive performance of a foil in a large amplitude sinusoidal motion, predictions were made for a man-made lunate planform to observe the differences between a rectangular and a lunate tail planform. This man-made lunate wing was taken from Chopra and Kambe (1977): the B2 planform. In their work, the offsets of the foil were calculated from

再 冉 冊冎

c ⫽ Cr 1 ⫺

y s

2

冉冊

⫹ Ctip

y s

2

(5)

for the local chord length and

冉冊

xl ⫽ K

y s

2

(6)

for the leading edge offsets. When K is taken as 5/4, s as 3, Cr ⫽ 1.0 and Ctip ⫽ 0.25 (the tip chord length), the foil has an aspect ratio of AR ⫽ 6s(2Cr ⫹ Ctip) ⫽ 8 and an area of S ⫽ 23 s(2Cr ⫹ Ctip) ⫽ 4.5. Input data were taken for nine stations and recorded on the input worksheet; then they were curve fitted by a spline function in a geometry interpolation program (JobContr.FOR). The sectional offsets were taken as a NACA 0012 designation (with 41 stations) and were interpolated by the same program. The right graph in Fig. 1 shows a perspective view of the foil with cosine and uniform spacing in chordwise and spanwise directions, respectively. The upper left graph shows a 3-D view of the right-most (starboard) tip section. The lower left graph shows the top view of this tip section. Calculations were done based on the following geometry and motion parameters: 앫 heave amplitude hheave ⫽ 1.0 m; hence, the factor h0 ⫽ hheave/Cr ⫽ 1.0; 앫 forward swimming velocity Vswim ⫽ 1.0 m/s; then, the velocity factor V*swim ⫽ 1.0/s;

524

P. Liu, N. Bose / Ocean Engineering 26 (1999) 519–529

Fig. 1. 3-D geometry of a man-made lunate planform with an aspect ratio of 8. This 3-D surface panel is illustrated by the perspective view.

앫 앫 앫 앫 앫

pitching axis x*pitch ⫽ xpitch/Cr ⫽ 1.0; small amplitude feathering parameter ␪ ⫽ 0.0, 0.4; large amplitude feathering parameter ⌰ ⫽ 0.0, 0.5093; angular velocity ␻ ⫽ 0.2 苲 1.8 to obtain k ⫽ 0.2 苲 1.8; pitching amplitude ␣0 ⫽ 0.08 苲 0.72 rad.

It is noted that the value of ⌰ corresponding to ␪ ⫽ 0.4 is not uniform for each curve by the present method because of the non-linear effect of the large amplitude oscillation. This gives for ␻ ⫽ 0.2, 0.6, 1.0, 1.4 and 1.8, values of ⌰ of 0.4093, 0.4441, 0.5093, 0.5892, 0.6769, respectively. For this man-made lunate foil, two comparisons were done and are presented as follows. 1. The present method versus previous ones. The present method takes the large amplitudes and the sectional thickness into account. The computations were done in the time domain. The previous methods are based on small amplitude lifting surface theory (zero sectional thickness). 2. Propulsive efficiency and thrust from a rectangular oscillating wing versus those from the man-made lunate planform predicted by the present method. Fig. 2 shows the propulsive efficiency predicted by a small amplitude lifting surface theory (Chopra and Kambe, 1977) and by the present method. In the calculation, this lunate tail propulsor was set to have a NACA 0012 section. The present method gives a lower efficiency prediction, especially at a large frequency, because of large amplitude considerations and the inclusion of skin friction in the present method. In addition to the large difference in the feathering parameters which exists between the two approaches, for a small amplitude feathering parameter ␪ ⫽ 0.4, when the reduced frequency is k ⫽ 1.8, the pitch amplitude ␣0 is 0.72 rad whereas for the large amplitude theory the pitch amplitude is 0.426 rad. The error for the pitch amplitude estimation itself from the small amplitude assumption is about 50%. There is also

P. Liu, N. Bose / Ocean Engineering 26 (1999) 519–529

Fig. 2.

525

Efficiency ␩ of a man-made lunate planform with an aspect ratio of 8.

about the same amount of error in heave amplitude. Therefore, the discrepancy in efficiency at larger heave and pitch amplitudes is significant. For a small feathering parameter, for example, ␪ ⫽ 0.0, only the error of heave amplitude exists; hence, the difference in efficiency is smaller. In Fig. 3, the difference in predictions of thrust coefficient by the two methods is substantial: lifting surface theory predicted a much larger thrust. The slope of the thrust curve obtained by the lifting surface theory at a higher pitch amplitude (the same as a higher reduced frequency k at a fixed feathering parameter ␪), is steep. In the lifting surface theory, the steepness is due to the linear effect of the lift coefficient Cl slope and the increasing thrust portion (the x-component of the normal force), as ␻ increases. Such a large difference in the thrust coefficient is not surprising, because the angle of attack of the small amplitude theory is almost twice as large as the large amplitude theory one.

Fig. 3.

Thrust of a man-made lunate planform with an aspect ratio of 8.

526

P. Liu, N. Bose / Ocean Engineering 26 (1999) 519–529

Table 1 Comparison of the hydrodynamic angles of attack at ␻t = ␲/2 and at the small amplitude feathering parameter ␪ = 0.4 k ␣␯l ␣␯p ␣al ␣ap

0.200 11.46° 11.31° ⫺11.46° ⫺11.61°

0.600 34.38° 30.96° 11.46° 8.04°

1.000 57.30° 45.00° 34.38° 22.08°

1.400 80.21° 54.46° 57.29° 31.54°

1.800 103.1° 60.95° 80.18° 38.03°

Table 1 gives the angles of attack of the foil calculated using the small amplitude and large amplitude theories at each reduced frequency k. The differences are the main reason for these large differences on the thrust shown in Fig. 4. That is, the larger the angle of attack, the higher the thrust, in a linear relation for a potential based theory. In this table, the pitch amplitude was set at 0.4 (22.92°). The normalized heave amplitude was h0 ⫽ 1.0, and the normalized swim velocity V0 ⫽ 1.0/s, based on the rootchord Cr. Therefore, in this case the oscillating frequency ␻ equals the reduced frequency k. Variables ␣␯l and ␣␯p are the angular kinematic velocity of the pitch axis for the lifting surface theory (small amplitude assumption) and the panel method, respectively. The values of ␣al are the effective angle of attack for the lifting surface theory due to the interpolation of the small amplitude theory, and ␣ap stands for the effective angle of attack for the present method. There is a big difference between these values at the same oscillating parameters at high reduced frequencies. These angles of attack were transient values at ␻t ⫽ ␲/2 (about the position of maximum angle of attack) and at a phase angle of 90° (heave lagging pitch). The amplitude of the instantaneous angle of attack varies with the oscillation of the foil. It can be seen that the angle of attack increases at high reduced frequencies greater than about k ⫽ 1.0 to values at which separation and vortex shedding would occur

Fig. 4.

Efficiency ␩ of a man-made lunate planform and a rectangular oscillating foil.

P. Liu, N. Bose / Ocean Engineering 26 (1999) 519–529

527

Table 2 Comparison of the hydrodynamic angles of attack at ␪ = 0.0 and at ␻t = ␲/2 k ␣␯l ␣␯p

0.200 0.2 rad 0.1974 rad

0.600 0.6 rad 0.5404 rad

1.000 1.0 rad 0.7854 rad

1.400 1.4 rad 0.9506 rad

1.800 1.8 rad 1.064 rad

on the oscillating foil. For attached flow, the range of reduced frequency would be less than about 1.0 or a value given of about k ⫽ 2␻Cr/V. To quantitatively examine the difference in thrust between these two methods, the angle of attack for the lifting surface and panel method were tabulated at different reduced frequencies for the small amplitude feathering parameter values of ␪ ⫽ 0.0 and 0.4 (Tables 2 and 3). At ␪ ⫽ 0.0 and thus a pitch amplitude ␣0 ⫽ 0.0, the angular kinematic velocity of the pitching axis is ␣␯l ⫽ ␻h0/V0 ⫽ ␻ by the small amplitude theory and

␣␯p ⫽ tan−1

冋 册

␻h0 ⫽ tan−1[␻] V0

by the present method. In this case, the angles of attack, ␣␯l and ␣␯p, at ␻t ⫽ ␲/2 equal the angular kinematic velocities of the pitching axis, because h0 and V0 were assumed to be 1.0. It can be seen from Table 2 that, at the reduced frequency k ⫽ 1.0, the angle of attack for the small amplitude assumption is ␣␯l ⫽ 1.0 rad and this angle of attack corresponds to that at a reduced frequency of about ␻ ⫽ 1.8 using the present method. The thrust values, from these two methods at the corresponding values of the angle of attack are very close (see Fig. 2). When the feathering parameter has a larger value at ␪ ⫽ 0.4, the change in the angle of attack is nonlinear for present method. It is noted that by considering a shift of the oscillating frequency ␻ (or the reduced frequency k) in terms of the equivalence of the angle of attack, thrusts predicted by the present method and the linearized theory were in general agreement. To observe the effect of a lunate tail planform shape on efficiency ␩ and thrust coefficient Ct, results from a rectangular oscillating propulsor were obtained. This wing had a span of 6 m and a chord of 0.75 m. The forward swimming velocity Vswim was set as 0.75 m/s in order to have all the same motion parameters as used by Chopra and Kambe (1977).

Table 3 Comparison of the hydrodynamic angles of attack at ␪ = 0.4 and at ␻t = /2 k ␣␯l ␣␯p

0.200 0.12 rad 0.1174 rad

0.600 0.36 rad 0.3004 rad

1.000 0.60 rad 0.3854 rad

1.400 0.84 rad 0.3906 rad

1.800 1.08 rad 0.3437 rad

528

P. Liu, N. Bose / Ocean Engineering 26 (1999) 519–529

Figs. 4 and 5 show that the man-made lunate planform had a significant efficiency and a small thrust increase over those from the rectangular planform, at a very small feathering parameter (⌰ ⫽ ␪ ⫽ 0.0), throughout the reduced frequency range. At ⌰ ⫽ 0.0, the instantaneous angle of attack of the foil reaches a maximum and the foil is under a heavy load condition. The lunate tail planform had a better performance. At a medium feathering parameter ⌰ ⫽ 0.4, the instantaneous angle of attack is smaller (so that the thrust is smaller), the skin friction drag became dominant over the thrust, and the required pitch moment portion of the input power also increased relatively (hence the efficiency was decreased). This made the efficiency and thrust less than those from a rectangular foil over most of the reduced frequency range studied. Hydrodynamics in unsteady conditions are important and hydrodynamic efficiency cannot be inferred from calculations in steady conditions (e.g. Van Dam, 1986). Though the performance of swept planforms outdoes that of rectangular planforms in steady conditions (Van Dam, 1985), rectangular planforms were found to give a better performance (greater efficiency and thrust coefficient) when the performance was calculated using unsteady methods and they were relatively lightly loaded. The large amplitude assumptions used here gave a more complete picture of propulsive performance from an oscillating wing than small amplitude thin wing theory. Unsteady lifting surface theory may still be a useful numerical method for the analysis of unsteady foils, but for large amplitude pitch and heave motion modifications are required for accurate results. At high loadings (small feathering parameter) the lunate planform swept wing was found to have a significantly higher propulsive efficiency and a small thrust increase over a rectangular foil undergoing similar oscillatory motions. This conclusion follows (as do the other points) from a large amplitude theory which assumes attached flow around the wing section. It is recognized that as the instantaneous effective angle of attack of an oscillating foil increases, separation occurs on the foil surface and large vortices can be shed from both leading and trailing edges. Under these conditions further significant augmentation of the thrust can occur under favourable

Fig. 5.

Thrust coefficient Ct of a man-made lunate planform and a rectangular oscillating foil.

P. Liu, N. Bose / Ocean Engineering 26 (1999) 519–529

529

combinations of the oscillatory parameters. No attempt was made in this study to assess the influence of vortex shedding from the foil at these large angles of attack. Acknowledgements The authors thank the Natural Sciences and Engineering Research Council Canada for financial support. References Bose, N., Lien, J., 1989. Propulsion of a fin whale (Balaenoptera physalus): why the fin whale is a fast swimmer. Proc. R. Soc. London B237, 175–200. Chopra, M.G., Kambe, T.K., 1977. Hydromechanics of lunate-tail swimming propulsion, part 2. J. Fluid Mech. 79 (Part1), 46–69. Fremaux, C.M., Vijgen, P.M.H.W., van Dam, C.P., 1990. Parametric analysis of swept-wing geometries with sheared wing tips. AIAA Paper 90-3196, 13 pp. Freund, R., Golub, G., Nachtigal, N., 1991. Iterative solution of linear systems. Acta Numerica 1991, 57–100. Hess, J.L., Smith, A.M.O., 1964. Calculation of nonlifting potential flow about arbitrary three dimensional bodies. Journal of Ship Research 8 (2), 22–44. Hess, J.L., Smith, A.M.O., 1967. Calculation of potential flow about arbitrary bodies. In: Ku¨chenmann D. (Ed.), Progress in Aeronautical Sciences, 8. Pergamon Press, pp. 1–139. Katz, J., Plotkin, A., 1991. Low-speed Aerodynamics. McGraw-Hill, New York, 1991, 632 pp. Liu, P., 1996. A time-domain panel method for oscillating propulsors with both chordwise and spanwise flexibility. Ph.D. thesis, Memorial University of Newfoundland, St. John’s, Newfoundland, Canada, 226 pp. Liu, P., Bose, N., 1993. Propulsive performance of three naturally-occurring oscillating propeller planforms. Ocean Engineering 20 (1), 57–75. Liu, P., Bose, N., 1996. A time domain panel method for wings of arbitrary, flexible planform under large amplitude, sinusoidal motion. In: Proceedings of the 2nd International Conference on Hydrodynamics, Hong Kong, A.A. Balkema, Rotterdam, Brookfield, 1996, pp. 113–118. Liu, P., Bose, N., 1997. Propulsive performance from oscillating propulsors with spanwise flexibility. Proc. R. Soc. London 453, 1963, 1763-1770. Maskew, B., 1987. Program VSAERO theory document, NASA CR-4023, September, 100 pp. Moran, J., 1984. Theoretical and Computational Aerodynamics. John Wiley and Sons, New York, 464 pp. Newman, J.N., Wu, T.Y., 1974. Hydromechanical aspect of fish swimming. In: Wu, T.Y., Brokaw, J., Brennen, C. (Eds.), Swimming and Flying in Nature, II. Plenum Press, pp. 615–634. Newman, J.N., 1986. Distributions of sources and normal dipoles over a quadrilateral panel. Journal of Engineering Mathematics 20, 113–126. Smith, M.J.C., 1996. Simulating moth wing aerodynamics: towards the development of flapping-wing technology. AIAA Journal 34 (7), 1348–1355. Theodorsen, T., 1935. General theory of aerodynamic instability and the mechanism of flutter. NACA Report No. 496, Washington DC, pp. 413–433. Van Dam, C.P., 1985. Swept wing tips shapes for low-speed airplanes. In: Aerospace Technology Conference and Exposition, Long Beach, CA, October, 10 pp. Van Dam, C.P., 1986. Drag-reduction characteristics of aft-swept wing tips. AIAA Paper 86-1824, pp. 327–338. Vest, M.S., Katz, J., 1996. Unsteady aerodynamic model of flapping wings. AIAA Journal 34 (7), 1433–1447. Von Ka´rma´n, T., Sears, W.R., 1938. Airfoil theory for non-uniform motion. J. Aeronaut. Sci. 5 (10), 379–390.