Derivation of ODEs and Bifurcation Analysis of a Two-DOF Airfoil Subjected to Unsteady Incompressible Flow
An airfoil subjected to two-dimensional incompressible inviscid flow is considered. The airfoil is supported via a translational and a torsional springs. The aeroelastic integro-differential equations of motion for the airfoil are reformulated into a system of six first-order autonomous ordinary differential equations. These are the simplest and least number of ODEs that can present this aeroelastic system. The differential equations are then used for the bifurcation analysis of an airfoil with a structural nonlinearity in the pitch direction. Sample bifurcation diagrams showing both stable and unstable limit cycle oscillation are presented. The types of bifurcations are assessed by evaluating the Floquet multipliers. For a specific case, a period doubling route to chaos was detected, and mildly chaotic behavior in a narrow range of velocity was confirmed via the calculation of the Lyapunov exponents.
A great deal of qualitative information can be obtained about wing flutter by studying the aeroelasticity of a simple two degree-of-freedom-system (DOF). The system, shown in Figure 1, consists of a rigid two-dimensional airfoil supported by torsional and translational springs. The translational and torsional springs are representative of the wing's flexural and torsional stiffnesses, respectively. The airfoil is subjected to two-dimensional incompressible inviscid flow. This model has been used frequently by researchers to obtain a better understanding of flutter [1–6].
The model is an extreme simplification of a wing. Despite this simplification, the underlying dynamics of the model could still be quite complicated if a structural nonlinearity is taken into account. The nonlinear dynamic behavior of this model has been investigated by several researchers, and limit cycle and chaotic oscillations have been detected for velocities below the linear flutter speed [1–4]. Previously, in order to take the advantages of available nonlinear dynamics techniques for the analysis of this simplified aeroelastic system, the governing integro-differential equations were converted to eight first-order ordinary differential equations (ODE) . The reformulated equations are autonomous ODEs, that is, they do not depend explicitly on the independent variable . The bifurcation analysis was then performed for the airfoil containing a 0.5 degree freeplay nonlinearity in the pitch stiffness. Sedaghat et al.  also converted the governing equations into eight nonautonomous ODEs and used normal form theory (NFT) to predict and characterize limit cycle oscillations (LCO) of the system.
In the present analysis, the integro-differential equations are further simplified to a set of six autonomous first-order ODEs. The new set of equations are much simpler than the previously presented equations [5, 6]. A bifurcation analysis is also performed for the airfoil containing a cubic nonlinearity in the pitch stiffness. Both defect-controlled and collocation methods are employed for the analysis [7, 8].
2. Derivation of ODEs from Integro-Differential Equations
where and denote differentiation with respect to nondimensional time ; and are nonlinear functions representing the restoring force and moment in the plunge and pitch directions, respectively, normalized with respect to the linear stiffnesses (in the case of a linear system simply and ); and are the nondimensional aerodynamic force and moment defined as Other parameters are defined in the nomenclature.
For incompressible flow the aerodynamic force and moment may be obtained for any arbitrary motion of the airfoil , giving where It is convenient to convert these integro-differential equations into a set of ordinary differential equations. The conversion could simply be accomplished if the approximate formula for the Wagner function is considered, where and .
Previously , using this approximate formula, the aeroelastic equations (1) and (3) have been transformed into eight autonomous first-order ODEs by twice integration of the equations. Also, in  the aeroelastic equations are transformed into eight nonautonomous ODEs via introducing 4 auxiliary variables for the integral terms, which is equivalent to ten autonomous ODEs. The auxiliary variables in Refeence  are ,, and . In the present analysis, the equations are converted into six autonomous ODEs with a better selection of two auxiliary variables as The following two auxiliary differential equations could then be constructed:
By inspection, the differential equations (7) can be written in the following simple forms By substituting and from (6) into (3), the aerodynamics terms and become After some algebra the equations of motion (1) are reformulated into the following second-order differential equations: These two second ODEs are much simpler than the set of ODEs presented in the previous publications [5, 6]. An important point to note regarding these equations is that all the coefficients are functions of airfoil and air parameters, and they are independent of time. Therefore, the reformulated equations are nonlinear autonomous ODEs and they do not depend explicitly on the independent variable . The reformulated equations (10) are then transformed into six first-order autonomous ODEs, in order to be able to analyze them using well-known analytical or numerical methods.
The above set of ODEs are either analyzed using collocation method via AUTO software package  or integrated numerically using defect-controlled method. The defect-controlled method employed here is an 8th-order Runge-Kutta method , which is known to be a reliable tool for detecting chaotic oscillations. The approximate Floquet multipliers are obtained by applying a standard eigenvalue routine to the approximation to the linearized Poincaré map.
3. Bifurcation Analysis Taking into Account a Cubic Nonlinearity in Pitch Stiffness
The linear flutter velocity , is first obtained by setting , and in (10) and solving the resultant eigenvalue problem. Stability of the linear system depends on the six eigenvalues of the system. The linearized set of equations is solved for the following parameters; ,,,. Figure 2 shows the variation with nondimensional velocity, , of the real part of the eigenvalues obtained for this set of parameters. Eigenvalues with large negative real parts are not shown in the figure.
The nonlinear set of ODEs are used for bifurcation analysis of the system taking into account a cubic nonlinearity of the form in the pitch direction and linear stiffness in plunge. The fixed points of the system can be calculated by setting all speeds and accelerations to zero in (10). In addition, the auxiliary variables and are set to zero as they are functions of speeds and accelerations.
This yields For , the only real solution is . A bifurcation diagram showing both stable and unstable solutions is presented in Figure 3. As shown, the origin is a stable fixed point until , point 1. At this point a supercritical Hopf bifurcation occurs, destabilizing the origin, and giving rise to a stable limit cycle. This stable periodic solution becomes unstable through a pitchfork bifurcation at point 2, and it changes its direction at the limit or turning point 3. Decreasing , the periodic solution remains unstable until it reaches a second turning point (point 5) where it restabilizes. It then remains stable for velocities greater than the linear flutter speed.
Another unstable periodic branch starts at point 2 and changes its direction and becomes stable at a limit point 6. This branch becomes unstable after a bifurcation to an invariant torus at point 7, reaches the limit point 4, turns back to a stable solution via another bifurcation into a torus at point 8. Passing the limit point 9, this branch becomes unstable and returns to point 2 where it started.
As explained above, a pitchfork bifurcation and two bifurcations into torus were detected for this case. The variations of Floquet multipliers around these bifurcation points are shown in Figure 4. At point 2 of Figure 3, the stable branch becomes unstable and two other unstable branches start from this point. Also, Figure 4(a) shows that at this point a Floquet multiplier crosses the unit circle at +1. These are characteristic of a pitchfork bifurcation . The variation of Floquet multipliers about point 7 is also presented in Figure 4(b). In this case two complex conjugate multipliers cross the unit circle; this means that the periodic orbit looses its stability to the torus. The phenomenon of bifurcation into a torus is sometimes called Hopf bifurcation of periodic orbit, secondary Hopf bifurcation or generalized Hopf bifurcation, and this bifurcation can lead to chaos . Although chaos was not obtained for this set of airfoil parameters, changing the parameters can lead to chaos through the torus bifurcations.
Figure 5 represents the bifurcation diagram obtained using the defect-controlled method along with the stable results given by AUTO presented in Figure 3. The results of the defect-controlled method show the value of when for the particular initial conditions given by and ; while AUTO results show only the maximum value of the pitch angle, . Both methods predict the same magnitude of the stable limit cycle and the same supercritical Hopf bifurcation point. Other tests using the defect-controlled method for various initial conditions show that both methods are also in agreement with the location of the other bifurcation points. For regions where more than one stable solution exists, the defect-controlled scheme converges to one or another solution depending on the initial conditions.
A bifurcation diagram was also constructed using AUTO for the same airfoil but with and . The bifurcation diagram and an expanded view of a partial region are shown in Figure 6. Once again, the origin is stable for low velocities; then, at it becomes unstable, but in this case through a subcritical Hopf bifurcation. The unstable limit cycle, which starts at point 1, gains stability at the limit point 2. This branch of periodic solution looses its stability at point 3 through a pitchfork bifurcation, to finally become stable at the limit point 5. This branch is stable for speeds much higher than the linear flutter boundary. Two other unstable periodic solutions start at point 3 and become stable at the limit points 6 and 9. Both branches of stable limit cycle undergo cascades of period doubling bifurcations. At first, period doubling bifurcations occur at points 7 and 8. Stable period-two solutions start at these bifurcation points to undergo further period doubling bifurcations at points 10 and 11. These period doubling bifurcations continue and lead finally to chaotic oscillations. All the unstable branches of the period doubling cascade end at point 4, which is an intersection point with the main branch of periodic solution 1-2-3-4-5.
The period doubling cascade leading to chaos was also predicted in the defect-controlled results. To compare the results obtained using two different methods, both the stable results obtained using AUTO and the results obtained using the defect-controlled method are presented in Figure 7. It is clear that both methods are in excellent agreement in the amplitude of pitch motion and the starting point of the chaotic region. However, the results obtained via AUTO not only give a more complete bifurcation diagram but also verifies the route to chaos.
Lyapunov spectra were also calculated for the same airfoil and nonlinearity as Figure 6. Figure 8(a) shows time evolutions of the largest Lyapunov exponent at in the chaotic region, and in the periodic region of Figure 6. A zero value for the largest Lyapunov exponent at clearly indicates periodic motion; while, a small positive exponent at indicates mildly chaotic behavior of the airfoil motion at this airspeed. In this paper “mildly chaotic” refers to a motion with small positive Lyapunov exponent indicating small deviation from periodic motion. The variation of the largest Lyapunov exponent with respect to is also shown in Figure 8(b). It is clear that for a narrow range of speed, approximately, there is a positive Lyapunov exponent indicating mildly chaotic behavior. This is also in agreement with the bifurcation diagrams obtained via both AUTO and the defect-controlled method. In this paper "mildly chaotic" refers to a chaotic motion with a small positive Lyapunov exponent.
Assuming the approximate exponential expressions for the Wagner function, the aeroelastic equations of motion for a two-DOF airfoil are conveniently transformed into a set of six autonomous first-order ODEs; then, the equations are analyzed using defect-controlled method and AUTO software package. Bifurcations diagrams have been obtained for the nonlinear system showing both stable and unstable solutions such as fixed points or limit cycle oscillations. Interesting bifurcation points and a period doubling route to chaos are predicted via calculating Floquet multipliers and Lyapunov exponents. The diagrams show that the system can have several stable and unstable solutions for a set of parameters. One may not be able to detect all the possible solutions experimentally or using time marching methods, unless one tries large number of initial conditions.
|:||Nondimensional distance measured from airfoil mid-chord to elastic axis|
|:||Semichord of airfoil|
|:||Nonlinear structural restoring force|
|:||Plunge motion of the airfoil|
|:||Mass moment of inertia about elastic axis|
|:||Linear structural stiffness in heave|
|:||Linear structural stiffness in pitch|
|:||Aerodynamic lift force|
|:||Mass of airfoil per unit span|
|:||Nonlinear structural restoring moment|
|:||Aerodynamic pitching moment about elastic axis|
|:||Nondimensional aerodynamic force|
|:||Nondimensional aerodynamic moment|
|:||Nondimensional radius of gyration about elastic axis|
|:||Nondimensional free stream velocity,|
|:||Nondimensional linear flutter velocity|
|:||Free stream velocity|
|:||Nondimensional distance from elastic axis to center of mass|
|:||Pitch rotation of the airfoil|
|:||Viscous damping ratio in pitch|
|:||Viscous damping ratio in heave|
|:||Airfoil-air mass ratio,|
|:||Nondimensional heave displacement,|
|:||Uncoupled frequency in pitch,|
|:||Uncoupled frequency in heave,|
The authors would like to acknowledge the financial support of the Natural Sciences and Engineering Research Council of Canada.
B. H. K. Lee and J. Desrochers, “Flutter analysis of a two-dimensional airfoil containing structural nonlinearities,” Aeronautical Report LR-618 27833, National Research Council of Canada, 1987.View at: Google Scholar
T. O'Neil and T. W. Strganac, “Aeroelastic response of a rigid wing supported by nonlinear springs,” Journal of Aircraft, vol. 35, no. 4, pp. 616–622, 1998.View at: Google Scholar
H. Alighanbari and S. J. Price, “The post-Hopf-bifurcation response of an airfoil in incompressible two-dimensional flow,” Nonlinear Dynamics, vol. 10, no. 4, pp. 381–400, 1996.View at: Google Scholar
A. Sedaghat, J. E. Cooper, J. R. Wright, and A. Y. T. Leung, “Prediction of non-linear aeroelastic instabilities,” in Proceedings of the 22nd Congress of the International Council of the Aeronautical Sciences (ICAS '00), pp. 464.1–464.10, Harrogate, UK, 2000.View at: Google Scholar
R. M. Corless, “Defect-controlled numerical methods and shadowing for chaotic differential equations,” Physica D, vol. 60, no. 1–4, pp. 323–334, 1992.View at: Google Scholar
E. J. Doedel and J. P. Kernevez, “AUTO: software for continuation and bifurcation problems in ordinary differential equations,” Applied Mathematics Report, California Institute of Technology, 1986.View at: Google Scholar
Y. C. Fung, An Introduction to the Theory of Aeroelasticity, Dover, New York, NY, USA, 1993.
R. Seydel, From Equilibrium to Chaos, Elsevier, New York, NY, USA, 1988.