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

Authors: H. Alighanbari; S. M. Hashemi
Journal: International Journal of Aerospace Engineering
 (2009)

---

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


---

## Body

## 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 Figure1, 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 Figure1 may be written in nondimensional form as [5],(1)ξ′′(τ)+xαα′′(τ)+2ζξω¯Uξ′(τ)+(ω¯U)2F(ξ)=p(τ),xαrα2ξ′′(τ)+α′′(τ)+2ζα1Uα′(τ)+1U2M(α)=r(τ),
where ξ=h/b and ()′ denote differentiation with respect to nondimensional time τ; F(ξ) and M(α) 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 M(α)=α and F(ξ)=ξ); p and r are the nondimensional aerodynamic force and moment defined as
(2)p(τ)=-LmV2/b,r(τ)=-MαmV2rα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
(3)L(τ)=πρbV2[ξ′′(τ)-ahα′′(τ)+α′(τ)+2(XTM)],Mα(τ)=πρb2V2[ahξ′′(τ)-(18+ah2)α′′(τ)-(12-ah)α′(τ)+2(12+ah)(XTM)],
where
(4)(XTM)=ϕ(τ)[ξ′(0)+(12-ah)α′(0)+α(0)]+∫0τϕ(τ-σ)[ξ′′(σ)+(12-ah)α′′(σ)+α′(σ)]dσ.
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)ϕ(τ)=l-ae-bτ-ce-dτ.
where l=1,a=0.165,b=0.0455,c=0.335 and d=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 y1=∫0τe-b(τ-σ)α(σ)dσ,y2=∫0τe-d(τ-σ)α(σ)dσ,y3=∫0τe-b(τ-σ)ξ(σ)dσ, and y4=∫0τe-d(τ-σ)ξ(σ)dσ. In the present analysis, the equations are converted into six autonomous ODEs with a better selection of two auxiliary variables as
(6)y1=ae-bτ[ξ′(0)+(12-ah)α′(0)+α(0)]+∫0τae-b(τ-σ)[ξ′′(σ)+(12-ah)α′′(σ)+α′(σ)]dσ,y2=ce-dτ[ξ′(0)+(12-ah)α′(0)+α(0)]+∫0τce-d(τ-σ)[ξ′′(σ)+(12-ah)α′′(σ)+α′(σ)]dσ.
The following two auxiliary differential equations could then be constructed:(7)y1′=-bae-bτ[ξ′(0)+(12-ah)α′(0)+α(0)]-b∫0τae-b(τ-σ)[ξ′′(σ)+(12-ah)α′′(σ)+α′(σ)]dσ+a[ξ′′(τ)+(12-ah)α′′(τ)+α′(τ)],y2′=-dce-dτ[ξ′(0)+(12-ah)α′(0)+α(0)]-d∫0τce-d(τ-σ)[ξ′′(σ)+(12-ah)α′′(σ)+α′(σ)]dσ+c[ξ′′(τ)+(12-ah)α′′(τ)+α′(τ)].
By inspection, the differential equations (7) can be written in the following simple forms
(8)y1′=-by1+a[ξ′′(τ)+(12-ah)α′′(τ)+α′(τ)],y2′=-dy1+c[ξ′′(τ)+(12-ah)α′′(τ)+α′(τ)].
By substituting y1 and y2 from (6) into (3), the aerodynamics terms p(τ) and r(τ) become
(9)p=1μ{-ξ′′+ahα′′-2lξ′-(1+2l(12-ah))α′-2lα+2y1+2y2{12}},r=1μrα2{ahξ′′-(18+ah2)α′′+2l(12+ah)ξ′-(12-ah)(1-2l(12+ah))α′+2l(12+ah)α-2(12+ah)y1-2(12+ah)y2}.
After some algebra the equations of motion (1) are reformulated into the following second-order differential equations:
(10)(1+1μ)ξ′′+(xα-ahμ)α′′+(2ζξω¯U+2lμ)ξ′+(1+2l(1/2-ah)μ)α′+(2lμ)α-(2μ)y1-(2μ)y2+(ω¯2U2)F(ξ)=0,(xαμ-ahrα2μ)ξ′′+(1+1/8+ah2rα2μ)α′′-(l(1+2ah)rα2μ)ξ′+(2ζαU-(2lah+l-1)(1/2-ah)rα2μ)α′-(l(1+2ah)rα2μ)α+(1+2ahrα2μ)y1+(1+2ahrα2μ)y2+(1U2)M(α)=0.
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 velocityU*, is first obtained by setting F(ξ)=ξ, and M(α)=α 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,ah=-0.5,xα=0.25,rα=0.5. Figure 2 shows the variation with nondimensional velocity, U, 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 ofU 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 formM(α)=0.1α+40α3 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 y1 and y2 are set to zero as they are functions of speeds and accelerations.This yields(11)(ω¯U)2ξ(τ)+2lμα(τ)=0,(-2U2l(ah+0.5)μrα2)α(τ)=M(α).
For ah=-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 U/U*=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 U/U*, 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) 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.
(a)(b)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 Figure4. 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.The variation of Floquet multipliers around two points in Figure3; (a) pitchfork bifurcation point 2, (b) bifurcation into torus point 7.
(a)(b)Figure5 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 M=0.01α+50α3. 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 U/U*=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) Bifurcation diagrams showing both stable and unstable solutions for pitch motion, (b) expanded view of partial region from (a);μ=200 and M(α)=0.01α+50α3; •, stable fixed point; ∘, unstable fixed point;  —, stable limit cycle; - - -, unstable limit cycle.
(a)(b)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 Figure7. 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 Figure6. Figure 8(a) shows time evolutions of the largest Lyapunov exponent at U/U*=0.475 in the chaotic region, and U/U*=0.5 in the periodic region of Figure 6. A zero value for the largest Lyapunov exponent at U/U*=0.5 clearly indicates periodic motion; while, a small positive exponent at U/U*=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 U/U* is also shown in Figure 8(b). It is clear that for a narrow range of speed, 0.455<U/U*<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) Time evolution of the Largest Lyapunov exponents, (b) variation of the largest Lyapunov exponent withU/U* for the same airfoil as Figure 6.
(a)(b)


## 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 ofsix  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.



---

*Source: 248930-2009-08-16.xml*