Derivation of ODEs and Bifurcation Analysis of a Two-DOF Airfoil Subjected to Unsteady Incompressible Flow

Abstract

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.

1. Introduction

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].

Figure 1 

Schematic of the airfoil.

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) [5]. 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. [6] 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

The equation of motion for the two-dimensional airfoil shown schematically in Figure 1 may be written in nondimensional form as [5],

𝜉^′′(𝜏)+𝑥_𝛼𝛼^′′(𝜏)+2𝜁_𝜉𝜔𝑈𝜉^′(𝜏)+(𝜔𝑈)²𝑥𝐹(𝜉)=𝑝(𝜏),_𝛼𝑟²_𝛼𝜉^′′(𝜏)+𝛼^′′(𝜏)+2𝜁_𝛼1𝑈𝛼^′1(𝜏)+𝑈²𝑀(𝛼)=𝑟(𝜏), (1) 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 𝐿𝑝(𝜏)=−𝑚𝑉²𝑀/𝑏,𝑟(𝜏)=−^𝛼𝑚𝑉²𝑟²_𝛼. (2) 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 [9], giving 𝐿(𝜏)=𝜋𝜌𝑏𝑉²[𝜉^′′(𝜏)−𝑎_ℎ𝛼^′′(𝜏)+𝛼^′𝑀(𝜏)+2(𝑋𝑇𝑀)],_𝛼(𝜏)=𝜋𝜌𝑏²𝑉²[𝑎_ℎ𝜉^′′1(𝜏)−(8+𝑎²_ℎ)𝛼^′′1(𝜏)−(2−𝑎_ℎ)𝛼^′1(𝜏)+2(2+𝑎_ℎ)(𝑋𝑇𝑀)], (3) where (𝑋𝑇𝑀)=𝜙(𝜏)[𝜉^′1(0)+(2−𝑎_ℎ)𝛼^′+(0)+𝛼(0)]_𝜏0𝜙(𝜏−𝜎)[𝜉^′′1(𝜎)+(2−𝑎_ℎ)𝛼^′′(𝜎)+𝛼^′(𝜎)]𝑑𝜎. (4) 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, 𝜙(𝜏)=𝑙−𝑎𝑒^−𝑏𝜏−𝑐𝑒^−𝑑𝜏. (5) where 𝑙=1,𝑎=0.165,𝑏=0.0455,𝑐=0.335 and 𝑑=0.3.

Previously [5], 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 [6] 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 [6] are 𝑦₁=∫^𝜏₀𝑒^{−𝑏(𝜏−𝜎)}𝛼(𝜎)𝑑𝜎,𝑦₂=∫^𝜏₀𝑒^{−𝑑(𝜏−𝜎)}𝛼(𝜎)𝑑𝜎,𝑦₃=∫^𝜏₀𝑒^{−𝑏(𝜏−𝜎)}𝜉(𝜎)𝑑𝜎, and 𝑦₄=∫^𝜏₀𝑒^{−𝑑(𝜏−𝜎)}𝜉(𝜎)𝑑𝜎. In the present analysis, the equations are converted into six autonomous ODEs with a better selection of two auxiliary variables as 𝑦₁=𝑎𝑒^−𝑏𝜏[𝜉^′1(0)+(2−𝑎_ℎ)𝛼^′+(0)+𝛼(0)]_𝜏0𝑎𝑒^{−𝑏(𝜏−𝜎)}[𝜉^′′1(𝜎)+(2−𝑎_ℎ)𝛼^′′(𝜎)+𝛼^′𝑦(𝜎)]𝑑𝜎,₂=𝑐𝑒^−𝑑𝜏[𝜉^′1(0)+(2−𝑎_ℎ)𝛼^′+(0)+𝛼(0)]_𝜏0𝑐𝑒^{−𝑑(𝜏−𝜎)}[𝜉^′′1(𝜎)+(2−𝑎_ℎ)𝛼^′′(𝜎)+𝛼^′(𝜎)]𝑑𝜎. (6) The following two auxiliary differential equations could then be constructed:

𝑦^′₁=−𝑏𝑎𝑒^−𝑏𝜏[𝜉^′1(0)+(2−𝑎_ℎ)𝛼^′(0)+𝛼(0)]−𝑏^𝜏₀𝑎𝑒^{−𝑏(𝜏−𝜎)}[𝜉^′′1(𝜎)+(2−𝑎_ℎ)𝛼^′′(𝜎)+𝛼^′(𝜎)]𝑑𝜎+𝑎[𝜉^′′1(𝜏)+(2−𝑎_ℎ)𝛼^′′(𝜏)+𝛼^′𝑦(𝜏)],_′2=−𝑑𝑐𝑒^−𝑑𝜏[𝜉^′1(0)+(2−𝑎_ℎ)𝛼^′(0)+𝛼(0)]−𝑑^𝜏₀𝑐𝑒^{−𝑑(𝜏−𝜎)}[𝜉^′′1(𝜎)+(2−𝑎_ℎ)𝛼^′′(𝜎)+𝛼^′(𝜎)]𝑑𝜎+𝑐[𝜉^′′1(𝜏)+(2−𝑎_ℎ)𝛼^′′(𝜏)+𝛼^′(𝜏)]. (7)(1)(2)(3)(4)(5)(6) By inspection, the differential equations (7) can be written in the following simple forms 𝑦^′₁=−𝑏𝑦₁+𝑎[𝜉^′′1(𝜏)+(2−𝑎_ℎ)𝛼^′′(𝜏)+𝛼^′𝑦(𝜏)],_′2=−𝑑𝑦₁+𝑐[𝜉^′′1(𝜏)+(2−𝑎_ℎ)𝛼^′′(𝜏)+𝛼^′(𝜏)]. (8) By substituting 𝑦₁ and 𝑦₂ from (6) into (3), the aerodynamics terms 𝑝(𝜏) and 𝑟(𝜏) become 1𝑝=𝜇{−𝜉^′′+𝑎_ℎ𝛼^′′−2𝑙𝜉^′1−(1+2𝑙(2−𝑎_ℎ))𝛼^′−2𝑙𝛼+2𝑦₁+2𝑦₂1},𝑟=𝜇𝑟²_𝛼{𝑎_ℎ𝜉^′′1−(8+𝑎²_ℎ)𝛼^′′1+2𝑙(2+𝑎_ℎ)𝜉^′1−(2−𝑎_ℎ1)(1−2𝑙(2+𝑎_ℎ))𝛼^′1+2𝑙(2+𝑎_ℎ1)𝛼−2(2+𝑎_ℎ)𝑦₁1−2(2+𝑎_ℎ)𝑦₂}. (9) After some algebra the equations of motion (1) are reformulated into the following second-order differential equations: 1(1+𝜇)𝜉^′′+(𝑥_𝛼−𝑎_ℎ𝜇)𝛼^′′+(2𝜁_𝜉𝜔𝑈+2𝑙𝜇)𝜉^′+(1+2𝑙(1/2−𝑎_ℎ)𝜇)𝛼^′+(2𝑙𝜇2)𝛼−(𝜇)𝑦₁2−(𝜇)𝑦₂+(𝜔²𝑈²(𝑥)𝐹(𝜉)=0,_𝛼𝜇−𝑎_ℎ𝑟²_𝛼𝜇)𝜉^′′+(1+1/8+𝑎_ℎ²𝑟²_𝛼𝜇)𝛼^′′−(𝑙(1+2𝑎_ℎ)𝑟²_𝛼𝜇)𝜉^′+(2𝜁_𝛼𝑈−(2𝑙𝑎_ℎ+𝑙−1)(1/2−𝑎_ℎ)𝑟²_𝛼𝜇)𝛼^′−(𝑙(1+2𝑎_ℎ)𝑟²_𝛼𝜇)𝛼+(1+2𝑎_ℎ𝑟²_𝛼𝜇)𝑦₁+(1+2𝑎_ℎ𝑟²_𝛼𝜇)𝑦₂1+(𝑈²)𝑀(𝛼)=0. (10) 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 [8] or integrated numerically using defect-controlled method. The defect-controlled method employed here is an 8th-order Runge-Kutta method [7], 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; 𝜔=0.2,𝜇=100,𝑎_ℎ=−0.5,𝑥_𝛼=0.25,𝑟_𝛼=0.5. 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.

Figure 2 

The real part of the eigenvalues as a function of 𝑈 for the linearized system.

The nonlinear set of ODEs are used for bifurcation analysis of the system taking into account a cubic nonlinearity of the form 𝑀(𝛼)=0.1𝛼+40𝛼³ 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 (𝜔𝑈)²𝜉(𝜏)+2𝑙𝜇(𝛼(𝜏)=0,−2𝑈²𝑙(𝑎_ℎ+0.5)𝜇𝑟²_𝛼)𝛼(𝜏)=𝑀(𝛼). (11) For 𝑎_ℎ=−0.5, the only real solution is 𝛼=𝜉=0. A bifurcation diagram showing both stable and unstable solutions is presented in Figure 3. As shown, the origin (𝛼=𝜉=0) is a stable fixed point until 𝑈/𝑈^∗=0.217, 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.

(a)

(b)

Figure 3 

(a) Bifurcation diagrams showing both stable and unstable solutions for pitch motion, (b) expanded view of partial region from (a). •, stable fixed point; ∘, unstable fixed point; —, stable limit cycle; - - - , unstable limit cycle.

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 [10]. 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 [10]. Although chaos was not obtained for this set of airfoil parameters, changing the parameters can lead to chaos through the torus bifurcations.

(a)

(b)

Figure 4 

The variation of Floquet multipliers around two points in Figure 3; (a) pitchfork bifurcation point 2, (b) bifurcation into torus point 7.

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 𝛼^′=0 for the particular initial conditions given by 𝛼(0)=7^∘ and 𝛼^′(0)=𝜉^′(0)=𝜉(0)=0; 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.

Figure 5 

Comparison of the bifurcation diagrams obtained via the defect-controlled method and AUTO. •, defect-controlled results; —, AUTO (stable solutions).

A bifurcation diagram was also constructed using AUTO for the same airfoil but with 𝜇=200 and 𝑀=0.01𝛼+50𝛼³. 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 𝑈/𝑈^∗=0.15 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.

(a)

(b)

Figure 6 

(a) Bifurcation diagrams showing both stable and unstable solutions for pitch motion, (b) expanded view of partial region from (a); 𝜇=200 and 𝑀(𝛼)=0.01𝛼+50𝛼3; •, stable fixed point; ∘, unstable fixed point; —, stable limit cycle; - - - , unstable limit cycle.

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.

Figure 7 

Comparison of the bifurcation diagrams obtained via the defect-controlled method and AUTO; •, defect-controlled results; —, AUTO (stable solutions).

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 𝑈/𝑈^∗=0.475 in the chaotic region, and 𝑈/𝑈^∗=0.5 in the periodic region of Figure 6. A zero value for the largest Lyapunov exponent at 𝑈/𝑈^∗=0.5 clearly indicates periodic motion; while, a small positive exponent at 𝑈/𝑈^∗=0.475 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, 0.455<𝑈/𝑈^∗<0.485 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.

(a)

(b)

Figure 8 

(a) Time evolution of the Largest Lyapunov exponents, (b) variation of the largest Lyapunov exponent with 𝑈/𝑈∗ for the same airfoil as Figure 6.

4. Conclusions

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.

Nomenclature

𝑎ℎ:	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, 𝑚/𝜋𝜌𝑏2
𝜉:	Nondimensional heave displacement, ℎ/𝑏
𝜏:	Nondimensional time, 𝑡𝑉/𝑏
𝜌:	Air density
𝜙:	Wagner's function
𝜔𝛼:	Uncoupled frequency in pitch, 𝐾𝛼/𝐼𝛼
𝜔𝜉:	Uncoupled frequency in heave, 𝐾ℎ/𝑚
𝜔:	Frequency ratio, 𝜔𝜉/𝜔𝛼

Acknowledgment

The authors would like to acknowledge the financial support of the Natural Sciences and Engineering Research Council of Canada.