The flight of a butterfly, Pieris melete, was observed in the take-off phase and was analyzed theoretically from aerodynamic and kinetic viewpoints. A vortex method, which was recently developed by the present authors, was used in this analysis.

During the downstroke, the butterfly generates mainly a vertical force. The acceleration of the butterfly’s body during the first half of the downstroke is especially large, and this acceleration is mainly caused by a large unsteady pressure drag acting on the wings. This large unsteady pressure drag is generated by the vortices shed into the flow from the outer edges of each wing of a pair; it is increased by the interference effect between a pair of wings when the opening angle is small. This force can be estimated by the previous quasi-steady analysis when the force coefficient is changed to 4. In addition to the unsteady pressure drag, an aerodynamic force due to added mass is generated and this is also increased by the interference effect between a pair of wings.

During the upstroke the butterfly generates mainly a horizontal force. The change of direction of the forces during the down- and upstrokes is controlled by variation in the inclination of the stroke plane. The moment, which is created by the aerodynamic force acting on the wings and by abdominal motion, changes the thoracic angle, that is the inclination of the stroke plane.

Keywords:
butterfly, Pieris melete, take-off flight, vortex method, ‘near fling’

A butterfly has low-aspect-ratio wings, which are not suitable for cruising flight. Much attention has been paid to the flight mechanism of such a low-aspect-ratio wing (Betts and Wootton, 1988; Dudley, 1990). A much larger force coefficient than the quasi-steady-state value was obtained by observations of flight in the field (Dudley, 1991). A butterfly uses a ‘peel’ mechanism. The unsteady vortices and the interference effect play important roles in this mechanism (Ellington, 1984a; Kingsolver, 1985; Brodsky, 1991). Brackenbury (1991) carried out careful observation of take-off and climbing flight in butterflies and showed that the hindwings and abdomen act to increase the interference effect between a pair of wings. Bocharova-Messner and Aksyuk (1981) suggested that the jet force due to a ‘tunnel’ between a pair of wings is important in the flight of a butterfly.

The characteristics of the three-dimensional ‘fling’ mechanism, which is equivalent to the ‘peel’ mechanism using rigid wings, were investigated. The interference effect between a pair of wings was made clear quantitatively and a numerical method of calculation was developed to analyze the variation of the pressure distribution on the wing (Sunada et al. 1993). In this paper, this numerical calculation technique is applied to a study of the take-off flight of the butterfly Pieris melete.

The take-off phase of the flight of the butterfly Pieris melete was observed and filmed with a high-speed video camera during the first period of beating motion. The geometrical characteristics of the butterfly are shown in Fig. 1 and Table 1.

Table 1.

Geometric characteristics of the butterfly

Geometric characteristics of the butterfly
Geometric characteristics of the butterfly
View Large
Fig. 1.
Schematic configuration of Pieris melete.
View largeDownload slide

Schematic configuration of Pieris melete.

Fig. 1.
Schematic configuration of Pieris melete.
View largeDownload slide

Schematic configuration of Pieris melete.

Because the geometrical relationship between the fore- and hindwings is changed during flight, the planform shapes observed in the photographs were averaged as shown in Fig. 1. In flight, the butterfly bends its body at a point G. This point G is in accordance with the centre of gravity of the body when the body is stretched. The body is divided into two parts, the thorax (including the head) and the abdomen, at the point G, as shown in Fig. 1. The recorded movements of the point G in the inertial frame are shown in Fig. 2. When the body bends, there is a small difference between the point G and the true centre of gravity of the body. Referring to Fig. 18, the position of the true centre of gravity is given as follows:

Fig. 2.
Time variation of the position of point G.
View largeDownload slide

Time variation of the position of point G.

Fig. 2.
Time variation of the position of point G.
View largeDownload slide

Time variation of the position of point G.

The thoracic angle ΘT and the abdominal angle ΘA are defined as the angles from the horizontal, as shown in Fig. 3A. The time variations of ΘT and ΘA are shown in Fig. 3B. This figure shows that the thorax is directed horizontally during the downstroke, and that it is vertical during the upstroke. A similar variation of the thoracic angle was observed qualitatively by Brackenbury (1991). The physical explanation of this variation will be given later in this paper. Near the beginning of downstroke (t<0.2T) and near the end of the upstroke (t>0.8T); the abdominal angle cannot be observed because the wing conceals the abdomen. The algebraic functions indicated by solid lines are estimated from discrete values of ΘT (circles) and ΘA (squares). In Fig. 3C the angular accelerations obtained by differentiating these functions are indicated by the solid line and the broken line ⁠.
Fig. 3.
Thoracic angle and abdominal angle. (A) Definition of ΘT and ΘA; (B) time variation of ΘT and ΘA; (C) time variation of Θ··T and Θ··A.
View largeDownload slide

Thoracic angle and abdominal angle. (A) Definition of ΘT and ΘA; (B) time variation of ΘT and ΘA; (C) time variation of Θ··T and Θ··A⁠.

Fig. 3.
Thoracic angle and abdominal angle. (A) Definition of ΘT and ΘA; (B) time variation of ΘT and ΘA; (C) time variation of Θ··T and Θ··A.
View largeDownload slide

Thoracic angle and abdominal angle. (A) Definition of ΘT and ΘA; (B) time variation of ΘT and ΘA; (C) time variation of Θ··T and Θ··A⁠.

The wing motions are described by two vectors of the fore- and hindwings which are represented by vectors and respectively, as shown in Fig. 1. The flapping angle β and the lead-lag angle ζ of and/or are defined by equations A4 and A5, when and/or is coincident with the x″ -axis. The coordinate systems used in this paper are defined in Appendix A and Abbreviations. These flapping angles are shown by circles and squares in Fig. 4A. It is observed that the phase of βf is always ahead of that of βh; that is, the flapping motion of is always ahead of that of ⁠.This means that the butterfly uses a ‘peel’ mechanism (Ellington, 1984a; Kingsolver, 1985). The flapping motion, where the forewings are always ahead of the hindwings, causes a ‘tunnel’ between the right and left wings near the beginning of the downstroke (Bocharova-Messner and Aksyuk, 1981). The ‘tunnel’, however, is not formed near the beginning of the upstroke, because the flapping angles, βf and βh, do not reach −π/2. This is different from the observation of Bocharova-Messner and Aksyuk (1981). The lead-lag angles for (circles) and (squares) are shown in Fig. 5. It is observed in Fig. 5 that the lead-lag angle of the hindwing ζh is roughly constant, ζh=−0.6, and that the lead-lag angle of the forewing ζf during the downstroke is larger than that of the upstroke. This means that the geometrical relationship between the fore- and hindwings is changed slightly during the flight. The butterfly moves the fore- and hindwings almost as one wing, so that the two are referred to as a ‘whole’ wing in this paper. The difference in the phase of the flapping angles, βf and βh, causes a twist of the ‘whole’ wing. This twist angle δ is defined as the angle between the vector and the flapping axis, as shown in Fig. 1. The definition of the vector is given in Appendix B. The feathering angle at x=0.7xtip is, then, defined as θfh=π/2− δ. The time variation of the feathering angle θfh and its analytical expression θ0.7(t) are shown in Fig. 6. It is observed that the feathering angle is almost π/2 and is larger during the upstroke than during the downstroke. It is assumed here that the wing is twisted linearly from root to tip. Thus, the feathering angle of a wing element at x is expressed as:

Fig. 4.
Flapping angle. (A) Flapping angle βf, βh and (t); (B) angular velocity, β·; (C) angular acceleration, β··.
View largeDownload slide

Flapping angle. (A) Flapping angle βf, βh and (t); (B) angular velocity, β·⁠; (C) angular acceleration, β··⁠.

Fig. 4.
Flapping angle. (A) Flapping angle βf, βh and (t); (B) angular velocity, β·; (C) angular acceleration, β··.
View largeDownload slide

Flapping angle. (A) Flapping angle βf, βh and (t); (B) angular velocity, β·⁠; (C) angular acceleration, β··⁠.

Fig. 5.
Lead-lag angle.
View largeDownload slide

Lead-lag angle.

Fig. 5.
Lead-lag angle.
View largeDownload slide

Lead-lag angle.

Fig. 6.
Feathering angle. θfh, observed values; θ0.7(t), analytical expression.
View largeDownload slide

Feathering angle. θfh, observed values; θ0.7(t), analytical expression.

Fig. 6.
Feathering angle. θfh, observed values; θ0.7(t), analytical expression.
View largeDownload slide

Feathering angle. θfh, observed values; θ0.7(t), analytical expression.

In order to calculate the aerodynamic and inertial forces acting on the wing, the motion of the ‘whole’ wing is determined as follows. The analytical expressions of three angles defined in Appendix A, β, ζ and θ, which determine the x,y,z coordinate system, are assumed to be equal to the measured values βf, 0 and θ, respectively. By differentiating these functions of time, the flapping angular velocity and the feathering angular velocity are obtained. These values near t=0 are extrapolated so that and become 0. By differentiating and again, the flapping angular acceleration and the feathering angular acceleration are obtained. The values of and are shown in Fig. 4B,C, respectively. The broken lines in these figures indicate the extrapolated part.
The geometrical angle of attack of the wing element at x, ⁠, is defined as the angle between the chord line and the inflow velocity at (x,y = 0). It is expressed as:
where
In the above equations, Hij is the i, j component of matrix H as defined in Appendix B. The calculated results of the geometrical angle of attack at three span positions, x= 0.25xtip, x=0.5xtip and x=0.75xtip, are shown in Fig. 7. It is observed that the geometrical angle of attack is larger than 1rad≈60° during most of the stroke. This means that drag force (the force parallel to the inflow velocity) is more dominant than lift force (the force perpendicular to the inflow velocity). Therefore, the numerical calculation method developed by the present authors (Sunada et al. 1993) is applied to the analysis of the flight of a real butterfly.
Fig. 7.
Geometrical angle of attack at three span positions. θ∧(0.25xtip), values at x=0.25xtip; θ∧(0.5xtip), values at x=0.5xtip; θ∧(0.75xtip), values at x=0.75xtip.
View largeDownload slide

Geometrical angle of attack at three span positions. θ∧(0.25xtip)⁠, values at x=0.25xtip; θ∧(0.5xtip)⁠, values at x=0.5xtip; θ∧(0.75xtip)⁠, values at x=0.75xtip.

Fig. 7.
Geometrical angle of attack at three span positions. θ∧(0.25xtip), values at x=0.25xtip; θ∧(0.5xtip), values at x=0.5xtip; θ∧(0.75xtip), values at x=0.75xtip.
View largeDownload slide

Geometrical angle of attack at three span positions. θ∧(0.25xtip)⁠, values at x=0.25xtip; θ∧(0.5xtip)⁠, values at x=0.5xtip; θ∧(0.75xtip)⁠, values at x=0.75xtip.

Experiments were performed to measure the fluid-dynamic forces and moment acting on a rotating plate in front of a mirror in a water tank. The shape of the plate was similar to that of a real butterfly wing as shown in Fig. 8 and Table 2. A small weight was attached to the plate in order to generate an initial rotating moment. The initial opening angle was set at 0. The non-dimensional distance d between the rotation axis of the plate and the mirror was selected to be d=0.02, 0.06, 0.09, 0.19 and ∞. The normal and tangential forces, FN and FS, and the moment around the rotation axis, Mf, acting on the plate during ‘near fling’ were measured. Details of the experimental apparatus and the procedure are given in Sunada et al. (1993).

Table 2.

Geometric characteristics of a test plate and a weight attached to the plate

Geometric characteristics of a test plate and a weight attached to the plate
Geometric characteristics of a test plate and a weight attached to the plate
View Large
Fig. 8.
Plate shape for fluid-dynamic tests.
View largeDownload slide

Plate shape for fluid-dynamic tests.

Fig. 8.
Plate shape for fluid-dynamic tests.
View largeDownload slide

Plate shape for fluid-dynamic tests.

Fig. 9A–C shows the time histories of opening angle α, angular velocity and angular acceleration ⁠. In Fig. 9B,C, the results of two particular cases are given; the results of other cases are located within these lines. The angular acceleration at t=0 becomes smaller as the distance d becomes smaller. The angular velocity at is common in all cases. The Reynolds number at is given by:

Fig. 9.
Motion of the plate. (A) Time history of opening angle α; (B) time history of angular velocity α·; (C) time history of angular acceleration α··.
View largeDownload slide

Motion of the plate. (A) Time history of opening angle α; (B) time history of angular velocity α·⁠; (C) time history of angular acceleration α··⁠.

Fig. 9.
Motion of the plate. (A) Time history of opening angle α; (B) time history of angular velocity α·; (C) time history of angular acceleration α··.
View largeDownload slide

Motion of the plate. (A) Time history of opening angle α; (B) time history of angular velocity α·⁠; (C) time history of angular acceleration α··⁠.

The above Reynolds number is almost equal to that of the flight of a real butterfly, ⁠. The normal and tangential forces are shown by the shaded area in Fig. 10A for all cases. The relationship between the moment around the rotation axis Mf and the opening angle α is common for all cases, and this relationship is shown by the bold line in Fig. 10B. The position of the centre of pressure Mf/FNxtip is shown in Fig. 10C. The values of FN are distributed in the shaded area of Fig. 10A. Hence, the average value of FN at each opening angle is used to determine the position of the centre of pressure. The position of the centre of pressure is located near 0.5xtip at the beginning of the motion and moves slightly towards the wing tip as the opening angle increases. A typical position of the centre of pressure may be represented by:
Fig. 10.
Results of fluid-dynamic tests. (A) Normal and tangential forces, FN and FS. (B) Moment around the rotation axis Mf. Hatched area or ––––––, experimental results; ––––––, results calculated using the vortex method (d=0.02); – ––, results calculated using the simple method (d=0.02). (C) Position of centre of pressure obtained experimentally.
View largeDownload slide

Results of fluid-dynamic tests. (A) Normal and tangential forces, FN and FS. (B) Moment around the rotation axis Mf. Hatched area or ––––––, experimental results; ––––––, results calculated using the vortex method (d=0.02); – ––, results calculated using the simple method (d=0.02). (C) Position of centre of pressure obtained experimentally.

Fig. 10.
Results of fluid-dynamic tests. (A) Normal and tangential forces, FN and FS. (B) Moment around the rotation axis Mf. Hatched area or ––––––, experimental results; ––––––, results calculated using the vortex method (d=0.02); – ––, results calculated using the simple method (d=0.02). (C) Position of centre of pressure obtained experimentally.
View largeDownload slide

Results of fluid-dynamic tests. (A) Normal and tangential forces, FN and FS. (B) Moment around the rotation axis Mf. Hatched area or ––––––, experimental results; ––––––, results calculated using the vortex method (d=0.02); – ––, results calculated using the simple method (d=0.02). (C) Position of centre of pressure obtained experimentally.

Fig. 11.
Configuration of the numerical calculation model.
View largeDownload slide

Configuration of the numerical calculation model.

Fig. 11.
Configuration of the numerical calculation model.
View largeDownload slide

Configuration of the numerical calculation model.

The solid lines in Fig. 10A,B show the results of the vortex method (Sunada et al. 1993) for a test case with d = 0.02. The configuration of the mathematical model is shown in Fig.11, where a plate shape similar to the real butterfly wing is represented by a simple form in order to simplify the calculation. The bound vortex is settled on the plate and the wake of the vortex sheet is on the surface swept by line DE. The total values of the normal force and the moment around the rotation axis, (FN)V and (Mf)V, consist of those created by dynamic pressure, (FN1)V and (Mf1)V, and those created by impulsive pressure, (FN2)V and (Mf2)V. It is observed in Fig.10A,B that the total values obtained by the present vortex method are in good agreement with the measurements of normal force and moment obtained experimentally. Using the relationship between the tangential force and the normal force due to dynamic pressure shown in Fig. 10A, it is assumed that the tangential force is given as:
The calculated tangential force is compared with the measured results in Fig.10A. Again, good agreement is obtained. Therefore, it is verified that by using the present calculation method it is possible to predict the normal and tangential forces and the moment around a rotation axis acting on a plate with a shape similar to that of a real butterfly wing.

The variation of shape factor VF, which is proportional to added mass (Sunada et al. 1993), with distance d is calculated by the vortex method and is shown in Fig. 12A. This figure shows the results with α =0. When the distance d becomes smaller, the shape factor becomes larger. The relationship between VF/VF(d = ∞) and distance d is independent of plate shape. It should be noted that VF(d = ∞) is a shape factor that has no effect on the interference between the right and left plates. The position of the centre of pressure is almost independent of distance d. This means that the variation of VM with distance d is similar to that of VF. The ratios between measured and calculated shape factors VF and VM are shown in Fig. 12B. Because these ratios are almost 1 for all values of d, the effectiveness of the present calculation method is verified for the real butterfly wing shape.

Fig. 12.
Shape factor of plate C, (α=0). (A) Shape factor VF, which is proportional to added mass, and position of centre of pressure versus non-dimensional distance d; (B) ratios between measured and calculated shape factors.
View largeDownload slide

Shape factor of plate C, (α=0). (A) Shape factor VF, which is proportional to added mass, and position of centre of pressure versus non-dimensional distance d; (B) ratios between measured and calculated shape factors.

Fig. 12.
Shape factor of plate C, (α=0). (A) Shape factor VF, which is proportional to added mass, and position of centre of pressure versus non-dimensional distance d; (B) ratios between measured and calculated shape factors.
View largeDownload slide

Shape factor of plate C, (α=0). (A) Shape factor VF, which is proportional to added mass, and position of centre of pressure versus non-dimensional distance d; (B) ratios between measured and calculated shape factors.

Simple method

Since the vortex method presented above requires a large amount of computational time, a simpler method is examined in this section. This simple method has been used to calculate fluid-dynamic forces acting on a beating wing (Weis-Fogh, 1973; Ellington, 1984b; Sunada et al. 1993). The method uses a plate element parallel to the rotation axis, as shown in Fig. 8.

The fluid-dynamic normal force acting on the plate element is composed of two kinds of forces. One is proportional to the second power of the inflow velocity, (FN1)S, and the other is proportional to the acceleration of the inflow, (FN2)S. Therefore:
For the present fluid-dynamic tests, the two values, Vz and ⁠, are given by:
In the above equations, the correction factors, k1 and k2, are introduced to improve the accuracy of the simple method, as explained later. The quantity VF(d,α)/VF(d = ∞) is also introduced in order to include the interference effect between right and left plates on added mass. By integrating the force and the moment around the rotation axis acting on the plate element along the x-axis, the normal force and the moment around the rotation axis acting on a whole plate are obtained with the following equations:
The value of k2 is determined with the following equation:
The quasi-steady force coefficients on the plate element (Sunada et al. 1993) are given by:
The value of k1 is defined as the mean of the above two values, and Therefore:
As in equation 7 of the vortex method, the tangential force (FS)S is given as:
Using equations 8–14, the calculation is made for the case d=0.02 of the fluid-dynamic tests. The results are shown by dashed lines in Fig. 10A,B. The results obtained by this simple method are in good agreement with the experimental results. Therefore, this simple method is shown to be a convenient and practical tool to use to estimate the aerodynamic forces acting on real butterfly wings.

Acceleration of the centre of gravity of the body in the take-off phase

Assuming that the left and right wings beat symmetrically, the following equation of motion of a butterfly is derived from Appendix B:
The first term in the right-hand side of this equation shows the gravitational force acting on the butterfly, and the second term shows the inertial force acting on a pair of ‘whole’ wings. The value of mtot is measured and the value of inertial force can be estimated by using the data of the motion analysis alone. The third term is aerodynamic force acting on a pair of ‘whole’ wings. By using equation 15, the aerodynamic forces estimated by the motion analysis, (Fair,X)M and (Fair,Z)M, are obtained:
These aerodynamic forces are also obtained theoretically by the two methods of numerical calculation, the simple method and the vortex method, as follows. The distance between the flapping axes of right and left ‘whole’ wings is 0.1xtip and the value of d is, then, 0.05 in the following calculation.

Simple method

The acceleration of the inflow is given by:
The forces acting on a wing element, dFN1 and dFN, are obtained by using equations 4, 8, 11, 13, 14 and 17. Referring to the relationships derived in Appendix A, the aerodynamic forces in the X and Z directions are calculated as follows:

Vortex method

In the vortex method, it is assumed that the ‘whole’ wing is planar and the feathering angle θ is θ0.7(t) for all positions of x. Three components of the inflow velocity due to wing and body motions on a point (x, y, 0) of a ‘whole’ wing are given by the following equations:
Adding the above inflow velocity, the velocity induced by bound vortices on the other wing and separated vortices from both the wings, the boundary condition, that is the z component of the velocity is 0, is satisfied. These numerically obtained forces, (FN)V and (FS)V, are transferred to the (X, Y, Z) axes by using the relationship derived in Appendix A, and the aerodynamic forces are expressed as follows:
where

Fig. 13A,B shows a comparison of the acceleration of the centre of gravity obtained experimentally with the calculated accelerations obtained by using the two numerical methods. The experimental values, and ⁠, are obtained by differentiating equation 1 as follows:

Fig. 13.
Acceleration of the centre of gravity of the butterfly’s body. (A) The −ZE direction; (B) the −XE direction.
View largeDownload slide

Acceleration of the centre of gravity of the butterfly’s body. (A) The −ZE direction; (B) the −XE direction.

Fig. 13.
Acceleration of the centre of gravity of the butterfly’s body. (A) The −ZE direction; (B) the −XE direction.
View largeDownload slide

Acceleration of the centre of gravity of the butterfly’s body. (A) The −ZE direction; (B) the −XE direction.

The abdominal angle ΘA was not measured during the periods 0<t<0.2T and 0.8T<t<T, as stated before. During these periods, and are assumed to be and ⁠, respectively. The calculated results, and are obtained from equation 15 by using (Fair,X)V, (Fair,Z)V, (Fair,X)S and (Fair,Z)S. It is shown in Fig. 13A,B that the results of the two numerical methods agree with the experimental results except for the first half of the upstroke. The agreement means that with the purely theoretical methods it is possible to get an understanding of the flight of a butterfly. It should be noted that the present theoretical models do not include any ‘tunnel’ effect proposed by Bocharova-Messner and Aksyuk (1981). It is also indicated that the difference between the results from the experiment and the vortex method is smaller than that between the results from the experiment and the simple method. Fig.14A,B shows −2(FI,X), −2(FI,Z), −2(Fair,X)M, −2(Fair,Z)M, −2(Fair,X)V, −2(Fair,Z)V and their components given by equations 20 and 21. The value of −2(Fair,Z)V agrees fairly well with that of −2(Fair,Z)M, and the value of −2(Fair,X)V agrees reasonably well with that of −2(Fair,X)M. These agreements on aerodynamic forces result in the agreements on the acceleration of the centre of gravity of the butterfly body. It is demonstrated in Fig. 14A,B that the aerodynamic forces due to dynamic pressure, −2(FN1,Z)V and −2(FN1,X)V, are dominant. The aerodynamic forces due to impulsive pressure, −2(FN2,Z)V and −2(FN2,X)V, are also large during the first half of the downstroke. The inertial forces, −2(FI,Z) and −2(FI,X), are much smaller than the aerodynamic forces. The upward vertical force is generated during the first half of the downstroke. It was pointed out by Sunada et al. (1993) that the aerodynamic forces are increased by the interference between a pair of wings when the opening angle is small. The butterfly utilizes the large forces for generating the large acceleration during the first half of the downstroke. The large acceleration may make it easier for the butterfly to escape from a predator. During the other period of the stroke, forward force is generated. This change of direction of force is due to the variation of the inclination of the stroke plane, as discussed later.
Fig. 14.
Each component of force acting on the butterfly’s wings. (A) The −Z direction; (B) the −X direction.
View largeDownload slide

Each component of force acting on the butterfly’s wings. (A) The −Z direction; (B) the −X direction.

Fig. 14.
Each component of force acting on the butterfly’s wings. (A) The −Z direction; (B) the −X direction.
View largeDownload slide

Each component of force acting on the butterfly’s wings. (A) The −Z direction; (B) the −X direction.

Thoracic angle

The equation of motion of a body rotating around the point G is given by Appendix C as follows:
In the above equation, all the terms except 2Mair are estimated from the data of the motion analysis. The value of 2Mair can be calculated from equation 23. Each component of this equation is shown in Fig. 15. As stated before, the abdominal angle is observable during the period 0.2T<t<0.8T. Therefore, the moment due to the aerodynamic force acting on the wings 2Mair is calculated only during this period. It is shown in Fig. 15 that the calculated moment 2Mair always raises the thorax. This is because the forewings are shifted forwards during the downstroke and backwards during the upstroke, as shown in Fig.5. The rise of the thorax is suppressed by the moment transmitted from the abdomen to the thorax, ⁠. Hence, the butterfly utilizes the motion of the abdomen to change the thoracic angle. The angle between the thoracic angle ΘT and the inclination of the stroke plane, ΘS, is constant, as shown in equation A3. Therefore, the above mechanism which controls the thoracic angle also controls the inclination of the stroke plane, and it causes the change of direction of force acting on the wing in the inertial frame. If the stroke plane is kept constant in space, the upward vertical force generated during the downstroke is cancelled by the downward vertical force during the upstroke. This is because drag force, which is parallel to the stroke plane, is dominant, as stated before. Therefore, the variation of the inclination of the stroke plane during both the strokes is the key mechanism of butterfly flight.
Fig. 15.
Moment around the Y-axis.
View largeDownload slide

Moment around the Y-axis.

Fig. 15.
Moment around the Y-axis.
View largeDownload slide

Moment around the Y-axis.

Power

The inertial torques around the flapping axis QI and the gravitational torque QG are given by:
The aerodynamic torque around the flapping axis QA is obtained by three methods: calculating QA ≡ (QA)V = Mf by the vortex method; calculating QA≡(QA)S = Mf by the simple method; or obtaining QA ≡ (QA)M from an analysis of the body motion of a butterfly. QA ≡ (QA)M is obtained as follows. By adopting the same relationships as equations 20 and 21, the normal force obtained by using only the data from the motion of the butterfly is given by:
It is assumed from equation 6 that the centre of pressure is located at x=0.56xtip. The aerodynamic torque (QA)M is obtained, therefore, as follows:
Fig.16 shows each component of the torque, QI, QG, ( QA)M, ( QA)V and (QA)S. Three aerodynamic torques, (QA)M, (QA)V and (QA)S, obtained with the above three methods show reasonable agreement with each other. The inertial and gravitational torques, QI and QG, are much smaller than the aerodynamic torque, (QA)M, or (QA)V or (QA)S. The necessary power required to beat the wings is given by:
Fig. 16.
Torque around the flapping axis.
View largeDownload slide

Torque around the flapping axis.

Fig. 16.
Torque around the flapping axis.
View largeDownload slide

Torque around the flapping axis.

There is little difference among these three values of the necessary power obtained by use of the three methods. The necessary power per 1kg body mass for the flight is estimated to be about 30 W kg−1. This value is roughly equal to the average value obtained by Dudley (1991) for the forward flight of butterflies.

Conclusion

New numerical calculation methods have been applied successfully to an analysis of take-off flight of a butterfly. The theoretical prediction agrees well with the observed butterfly motion. The acceleration observed during take-off flight is mainly due to the pressure drag generated by the wings in ‘near fling’ motion. The butterfly utilizes the interference effect in ‘near fling’ motion to generate the large acceleration during the first half of the downstroke. This large acceleration can help the butterfly to take off suddenly. The necessary power for this flight is estimated to be about 30 W kg−1 bodymass. This value is roughly equal to the average value obtained by Dudley (1991) for the forward flight of butterflies.

The flight mechanism of the butterfly is determined to be as follows. The butterfly changes the direction of the aerodynamic forces by changing the inclination of the stroke plane. This variation in the inclination of the stroke plane is mainly due to the moment created by the aerodynamic force acting on the wings and to the moment generated by abdominal motion. The aerodynamic moment raises the thorax and this moment is suppressed by abdominal motion during the second half of the downstroke and the first half of the upstroke in the observed flight. The butterfly shifts the forewings forward during the downstroke. This shift increases the wing area, that is the aerodynamic force, and generates the aerodynamic moment which raises the thorax.

Appendix A

The geometrical relationships among the various coordinate systems used in this paper, which are explained in Fig. 17, are discussed in this Appendix. An orthogonal Cartesian coordinate system (XE, YE, ZE) is fixed with respect to the earth. The ZE-axis is vertical and positive downwards. The coordinate system (X, Y, Z), the origin of which is located at the point G, is transferred in parallel from the coordinate system (XE, YE, ZE). Then,

Fig. 17.
Relationship between the coordinate systems.
View largeDownload slide

Relationship between the coordinate systems.

Fig. 17.
Relationship between the coordinate systems.
View largeDownload slide

Relationship between the coordinate systems.

Fig. 18.
Forces and moments acting on the thorax, abdomen and wing element. (A) On the thorax and abdomen; (B) on a wing element.
View largeDownload slide

Forces and moments acting on the thorax, abdomen and wing element. (A) On the thorax and abdomen; (B) on a wing element.

Fig. 18.
Forces and moments acting on the thorax, abdomen and wing element. (A) On the thorax and abdomen; (B) on a wing element.
View largeDownload slide

Forces and moments acting on the thorax, abdomen and wing element. (A) On the thorax and abdomen; (B) on a wing element.

The coordinate system (X, Y, Z) is coincident with the coordinate system (XE, YE, ZE) at the beginning of beating motion.
The coordinate system (XS, YS, ZS) is defined by rotating the coordinate system (X, Y, Z) about the Y-axis by ΘS as follows:
The angle ΘS is given by:
where Λ is the angle between the flapping axis and the thorax axis, as shown in Fig. 1, and is nearly approximated by Λ = 0.5. The flapping axis is defined by two points, the roots of the fore- and hindwings.
The coordinate system (XS, YS, ZS) is rotated about XS-axis by − β and XS, YS and ZS are then replaced by y′, x′ and −z′, respectively. The following equation is obtained:
The coordinate system (x″, y″, z″) is given by a rotation of, about the z′-axis as follows:
The coordinate system (x, y, z) is finally obtained by a rotation (θ−π/2) about the x″-axis:
The relationship between the coordinate systems (x, y, z) and (XE, YE, ZE), therefore, is given by:

Appendix B

Let (x, y, z) be a wing-fixed coordinate system in this Appendix as follows. The x-axis is perpendicular to the flapping axis and crosses the line BD. The flapping angle β of the x-axis is defined by equation A.4 and is set to be βf. The vector is determined to be perpendicular to the x-axis, as shown in Fig. 1. The y-axis is defined as being parallel to ⁠. The orthogonal Cartesian coordinate system (x, y, z) is thus determined.

A wing element on the right ‘whole’ wing shown in Fig. 18B is considered in order to obtain the equation of motion of the ‘whole’ wing. The centre of gravity of this wing element has the coordinate (x, yG(x), 0) with respect to the (x, y, z) axes.

The coordinate (XE, YE, ZE) of the centre of gravity of the wing element with respect to the earth-fixed axis is given by:
where
It should be noted that the above equation is derived from equation A7 in Appendix A by using the relationship that R is a unit matrix. By differentiating equation B1, the acceleration of the centre of gravity of the wing element is given by:
The external force acting on the wing element is composed of the aerodynamic force Δ Fair, the shear force transmitted from the neighbouring wing elements Δ S and the gravitational force. This external force generates the acceleration of the centre of gravity of the wing element as follows:
Integrating equation B4 on the right ‘whole’ wing:
where wing mass of unit area κ is assumed to be constant over the ‘whole’ wing. A similar relationship is obtained for the left ‘whole’ wing. Thus, the equation for both right and left wings is given by:
The position of the centre of gravity of the body G is defined as in equation A1. The external force acting on the body is composed of the gravitational force, the shear force transmitted from the right and left wings and the aerodynamic force. This aerodynamic force FB, which is generated by the body alone, is small in comparison with other terms. Therefore:
Substituting equation B7 into equation B6, the following equation of motion of the butterfly is obtained:

Appendix C

The matrix H in equation 4 is defined as follows:
Then, each component of the above matrix is given as follows:

Appendix D

Assuming that feathering angle θ of the wing element is π/2, the equation of angular motion about the Y-axis for a wing element is obtained as follows (Refer to Fig. 18):
where
Integrating equation D1 on the right ‘whole’ wing, the equation of rotating motion of the right ‘whole’ wing is given by:
where
The equations of angular motion around the point G are given for the thorax and the abdomen as follows:
where
Using equations D3 and D6, the following equation of rotating motion around the Y- axis is obtained for the butterfly:
     
  • c

    chord length (m)

    •  
  • d

    non-dimensional half-distance between the rotation (flapping) axes of right and left plates (wings), d0/xtip

    •  
  • d0

    half-distance between the rotation (flapping) axes of right and left plates (wings) (m)

    •  
  • dm

    mass of wing element (kg)

    •  
  • Fair

    aerodynamic force generated by a right ‘whole’ wing (N)

    •  
  • FB

    aerodynamic force generated by body (N)

    •  
  • FN

    normal component of aerodynamic force acting on a right ‘whole’ wing (N)

    •  
  • FS

    tangential component of aerodynamic force acting on a right ‘whole’ wing (N)

    •  
  • FAT

    force transmitted from thorax to abdomen (N)

    •  
  • FWT

    force transmitted from thorax to wing (N)

    •  
  • g

    acceleration due to gravity, 9.81 ms−2

    •  
  • H

    transformation matrix as defined by equation B2

    •  
  • Hij

    i, j component of matrix H

    •  
  • H−1

    inverse matrix of H

    •  

  • i, j component of matrix

    •  
  • IT

    moment of inertia of thorax around Y-axis (kgm2)

    •  
  • IA

    moment of inertia of abdomen around Y-axis (kgm2)

    •  
  • IY

    moment of inertia of a right ‘whole’ wing around Y-axis (kg m2)

    •  
  • Iy

    moment of inertia of a right ‘whole’ wing around y-axis (kgm2)

    •  
  • k1, k2

    correction factors in the simple method

    •  

  • quasi-steady force coefficients at

    •  
  • lT

    length of thorax (m)

    •  
  • lT,G

    distance between the centre of gravity of the thorax and point G as shown in Fig. 18A (m)

    •  
  • lA

    length of abdomen (m)

    •  
  • lA,G

    distance between the centre of gravity of the abdomen and point G as shown in Fig. 18A (m)

    •  
  • Mair

    aerodynamic moment acting on a right ‘whole’ wing around Y-axis (Nm)

    •  
  • MT,I

    inertial moment acting on thorax around Y-axis (Nm)

    •  
  • Mf

    aerodynamic moment around rotation axis (Nm)

    •  
  • MAT

    moment transmitted from thorax to abdomen around Y-axis (Nm)

    •  
  • MA,I

    inertial moment acting on abdomen around Y-axis (Nm)

    •  
  • MWT

    moment transmitted from thorax to a right ‘whole’ wing around Y-axis (Nm)

    •  
  • MW,I

    inertial moment acting on a right ‘whole’ wing around Y-axis (Nm)

    •  
  • mA

    abdominal mass (kg)

    •  
  • mB

    body mass (kg), mA + mT

    •  
  • mT

    thorax mass (kg)

    •  
  • mtot

    total mass of a butterfly (kg) mB + mW

    •  
  • mW

    wing mass (kg)

    •  
  • p

    transformation matrix between (XS, YS, ZS) and (X, Y, Z)

    •  
  • Pn

    necessary power (W)

    •  
  • Q

    transformation matrix between (XS, YS, ZS) and (x′, y′, z′)

    •  
  • QA

    aerodynamic torque around y-axis (Nm)

    •  
  • QG

    gravitational torque around y-axis (Nm)

    •  
  • QI

    inertial torque around y-axis (Nm)

    •  
  • R

    transformation matrix between (x′, y′, z′) and (x″, y″, z″)

    •  
  • Re

    Reynolds number

    •  
  • S

    transformation matrix between (x, y, z) and (x″, y″, z″)

    •  
  • T

    period of one beating cycle (s)

    •  
  • t

    time (s)

    •  
  • VF

    shape factor that is proportional to added mass (m4)

    •  
  • VM

    shape factor that is proportional to added moment of inertia (m5)

    •  
  • Vx, Vy, Vz

    x, y and z components of inflow velocity on the wing due to wing and body motion (ms−1)

    •  
  • (X, Y, Z)

    coordinate system which is parallel to the coordinate system (XE, YE, ZE) and the origin of which is located at the centre of gravity of the butterfly body

    •  
  • (XE, YE, ZE)

    earth-fixed coordinate system

    •  
  • (XG, YG, ZG)

    position of the point G

    •  
  • (x, y, z)

    wing-fixed coordinate system

    •  
  • yG

    chordwise position of centre of gravity of wing element (m)

    •  
  • α

    opening angle (rad)

    •  
  • β

    flapping angle (rad)

    •  
  • Δ

    value for wing element

    •  
  • ΔMS

    moment transmitted from neighbouring wing element to a wing element (Nm)

    •  
  • ΔS

    shear force transmitted from neighbouring wing elements to a wing element (N)

    •  
  • δ

    twist angle (rad)

    •  
  • ζ

    lead-lag angle (rad)

    •  
  • ΘA

    abdominal angle (rad)

    •  
  • ΘS

    flapping axis angle (rad), ΘT+ Λ

    •  
  • ΘT

    thoracic angle (rad)

    •  
  • θ

    feathering angle (rad)

    •  
  • θ0.7

    feathering angle defined by (rad)

    •  

  • geometrical angle of attack defined by equation 2 (rad)

    •  
  • κ

    wing mass of unit area (kg m−2)

    •  
  • Λ

    angle between flapping axis and thorax axis (rad)

    •  
  • ν

    kinematic viscosity (m2 s−1)

    •  
  • ξ

    angle defined in equation C2 (rad)

    •  
  • ρ

    density (kg m−3)

    Subscript or superscript

    •  
  • a

    air

    •  
  • f

    fluid or value defined by

    •  
  • h

    value defined by

    •  
  • I

    inertial component

    •  
  • ( )M

    value calculated from data of motion analysis alone

    •  
  • N

    normal force

    •  
  • S

    tangential force

    •  
  • ( )S

    value calculated by the simple method

    •  
  • tip

    wing tip

    •  
  • ( )V

    value calculated by the vortex method

    •  
  • w

    water

    •  
  • X, Y, Z

    X, Y, Z component

    •  
  • 1

    dynamic pressure

    •  
  • 2

    impulsive pressure

Betts
,
C. R.
 and
Wootton
,
R. J.
 (
1988
).
Wing shape and flight behavior in butterflies (Lepidoptera: Papilionoidea and Hesperioidea): A preliminary analysis
.
J. exp. Biol.
138
,
271
–
288
.
https://doi.org/10.1242/jeb.138.1.271
Google Scholar
Crossref
Search ADS
 
Bocharova-Messner
,
O. M.
 and
Aksyuk
,
T. S.
 (
1981
).
Production of a tunnel by the wings during flight in butterflies (Lepidoptera, Rhopalocera)
.
Dokl. Akad. Nauk SSSR
260
,
1490
–
1493
(in Russian).
Google Scholar
 
Brackenbury
,
J. H.
 (
1991
).
Kinematics of take-off and climbing flight in butterflies
.
J. Zool., Lond.
224
,
251
–
270
.
https://doi.org/10.1111/j.1469-7998.1991.tb04803.x
Google Scholar
Crossref
Search ADS
 
Brodsky
,
A. K.
 (
1991
).
Vortex formation in the tethered flight of the peacock butterfly Inachis io L. (Lepidoptera, Nymphalidae) and some aspects of insect flight evolution
.
J. exp. Biol.
161
,
77
–
95
.
https://doi.org/10.1242/jeb.161.1.77
Google Scholar
Crossref
Search ADS
 
Dudley
,
R.
 (
1990
).
Biomechanics of flight in neotropical butterflies: morphometrics and kinematics
.
J. exp. Biol.
150
,
37
–
53
.
https://doi.org/10.1242/jeb.150.1.37
Google Scholar
Crossref
Search ADS
 
Dudley
,
R.
 (
1991
).
Biomechanics of flight in neotropical butterflies: aerodynamics and mechanical power requirements
.
J. exp. Biol.
159
,
335
–
357
.
https://doi.org/10.1242/jeb.159.1.335
Google Scholar
Crossref
Search ADS
 
Ellington
,
C. P.
 (
1984a
).
The aerodynamics of flapping animal flight
.
Am. Zool.
24
,
95
–
105
.
https://doi.org/10.1093/icb/24.1.95
Google Scholar
Crossref
Search ADS
 
Ellington
,
C. P.
 (
1984b
).
The aerodynamics of hovering insect flight. II. Morphological parameters
.
Phil. Trans. R. Soc. Lond. B
305
,
17
–
40
.
https://doi.org/10.1098/rstb.1984.0050
Google Scholar
Crossref
Search ADS
 
Kingsolver
,
T. G.
 (
1985
).
Butterfly engineering
.
Scient. Am.
253
,
90
–
97
.
https://doi.org/10.1038/scientificamerican0885-106
Google Scholar
Crossref
Search ADS
 
Sunada
,
S.
,
Kawachi
,
K.
,
Watanabe
,
I.
 and
Azuma
,
A.
 (
1993
).
Fundamental analysis of three-dimensional ‘near fling’
.
J. exp. Biol.
183
,
217
–
248
.
https://doi.org/10.1242/jeb.183.1.217
Google Scholar
Crossref
Search ADS
 
Weis-Fogh
,
T.
 (
1973
).
Quick estimates of flight fitness in hovering animals, including novel mechanisms for lift production
.
J. exp. Biol.
59
,
169
–
230
.
https://doi.org/10.1242/jeb.59.1.169
Google Scholar
Crossref
Search ADS
 
© 1993 by Company of Biologists
1993