Model Reduction of Nonlinear Aeroelastic Systems Experiencing Hopf Bifurcation

In this paper, we employ the normal form to derive a reduced-order model that reproduces nonlinear dynamical behavior of aeroelastic systems that undergo Hopf bifurcation. As an example, we consider a rigid two-dimensional airfoil that is supported by nonlinear springs in the pitch and plunge directions and subjected to nonlinear aerodynamic loads. We apply the center manifold theorem on the governing equations to derive its normal form that constitutes a simplified representation of the aeroelastic system near flutter onset (manifestation of Hopf bifurcation). Then, we use the normal form to identify a self-excited oscillator governed by a time-delay ordinary differential equation that approximates the dynamical behavior while reducing the dimension of the original system. Results obtained from this oscillator show a great capability to predict properly limit cycle oscillations that take place beyond and above flutter as compared with the original aeroelastic system.


Introduction
The need to model, capture and predict properly the nonlinear dynamics, such as limit cycle oscillations, bifurcations, chaos, ... (Dowell et al., 2003;Ghommem et al., 2010a;Ghommem et al., 2012;Abdelkefi et al., 2012a,b;Vasconcellos et al., 2012) associated with aeroelastic systems may lead to large and complex models (Beran and Silva, 2004;Beran et al., 2004).The complexity associated with such models and extensive computational resources and time required to numerically integrate them constitute an obstacle for running a large number of simulations to test different configurations and quantify the effect of the system's parameters as required for design purposes.To overcome this issue, one needs to develop robust reduced-order models that enable fast and simplified simulations by while neglecting irrelevant physics and response characteristics.
In this paper, we follow a different approach based on observing a particular behavior or phenomenon from data obtained from experiments or numerical simulations and then describing them with existing models such as self-excited oscillators.Towards this end, we employ perturbations techniques to derive the normal form of the Hopf bifurcation of these aeroelastic systems.This form constitutes an effective tool to capture the main physical behaviors of such systems near the Hopf bifurcation (flutter onset) (Nayfeh, 1973;Abdelkefi et al., 2012c,d,e).Several methods have been proposed for deriving this form and then analyzing the nature of Hopf bifurcation.For instance, Nayfeh et al. (2012) and Abdelkefi et al. (2013) applied the method of multiple scales to construct an approximation to the response of different wings near the Hopf bifurcation.Dimitriadis et al. (2004) used the center manifold theorem to a nonlinear aeroelastic system to reduce its dimensionality and predict its bifurcation and post-bifurcation behavior.Leng (1995) used the method of averaging to derive a reduced-order model for a nonlinear aeroelastic system.The model was used then to conduct a bifurcation analysis and showed a great capability to capture the nonlinear aspects associated with the aeroelastic system.
In our present work, we consider a rigid two-dimensional airfoil that is supported by nonlinear springs in the pitch and plunge directions and subjected to nonlinear aerodynamic loads and apply the center manifold theorem on the governing equations to derive its normal form.Then, we use this form to come up with a simplified representation based on time-delay differential equation to reproduce its nonlinear dynamical behavior.These types of equations have been used to model the behaviors of many physical systems in physiology (Heiden, 1979), biology (MacDonald, 1989), and population dynamics (Kuang, 1993).The emphasis of our study is to predict properly the type of dynamic instability (Hopf bifurcation) and reproduce limit cycle oscillations (LCOs) that take place beyond flutter onset.

Modeling of the aeroelastic system
The aeroelastic system under investigation is a two-dimensional rigid airfoil constrained to move with two degrees of freedom, namely the plunge (h) and pitch (α) motions, as shown in Fig. 1.The governing equations of motion of this system are (Strganac et al., 1999;Gilliat et al., 2003) 00 00 () where T m is the total mass of the wing with its support structure, W m is the wing mass alone, I  is the mass moment of inertia about the elastic axis, b is the semichord length, and x  is the nondimensional distance between the center of mass and the elastic axis.The viscous damping forces are described through the coefficients h c and c  for plunge and pitch, respectively.Table 1 presents the values of all of the parameters used in the following analysis.In addition, L and M are the aerodynamic lift and moment about the elastic axis.These aerodynamic loads are assumed to be given by the quasi-steady representation with a stall model (Strganac et al., 1999;Ghommem et al., 2010b) and written as: where U is the freestream velocity, l c  and m c  are the aerodynamic lift and moment coefficients, and s c is a nonlinear parameter associated with stall.The effective angle of attack due to the instantaneous motion of the airfoil is given by (Strganac et al., 1999;Ghommem et al., 2010b) where a is the nondimensional distance from the midchord to the elastic axis.

Fig 1. Schematic of a two-dimensional rigid airfoil aeroelastic system
Abdessattar Abdelkefi and Mehdi Ghommem / Journal of Modeling, Simulation, Identification, and Control (2013) 1: 57-77 60 The two spring forces for the plunge and pitch motions are represented by h k and k  , respectively.
The representative parameters of these stiffnesses are approximated in polynomial form by To express the equations of motion in state space form, we define the following state of variables: The equations of motion are then rewritten as Abdessattar Abdelkefi and Mehdi Ghommem / Journal of Modeling, Simulation, Identification, and Control (2013) 1: 57-77 In vector form, the equations of motions are expressed as where C(X,X,X) is a cubic vector function of the state variables and the matrix A(U) is expressed as Abdessattar Abdelkefi and Mehdi Ghommem / Journal of Modeling, Simulation, Identification, and Control (2013) 1: 57-77

62
The matrix A(U) has a set of four eigenvalues i  , for i=1,2,3,4.Based on these values, we can evaluate the stability of the trivial solution of equation ( 9).These eigenvalues are complex conjugates

Normal form of the Hopf bifurcation of the aeroelastic system
To derive the normal form of the Hopf bifurcation of the aeroelastic system near the onset of instability f U , we add a perturbation term, Uf U  , to the flutter speed to the appearance of the secular terms at the third order.Consequently, the matrix A(U) can be expressed as:  10) can be rewritten as: Then, we multiply equation ( 11) from the left by the inverse where  12) in a component form as follows: where the According to the center-manifold theorem, there exists a center manifold Moreover, the dynamics of the system are similar to those of this center manifold; that is, Because U  is small and 3 N is a cubic function of the components of Y, 1 H is zero to the third approximation.Thus, keeping only the resonance terms in equation ( 15) (Nayfeh and Blachandran, 1995;Abdelkefi et al., 2012c), we obtain the complex-valued normal form where e N depends on the cubic nonlinear spring coefficients Abdessattar Abdelkefi and Mehdi Ghommem / Journal of Modeling, Simulation, Identification, and Control (2013) 1: 57-77 We express 3 Y in the polar form as follows: where a is the amplitude of oscillations and  is its phase and substituting Equation (17) into Equation ( 16) and separating the real and imaginary parts, we obtain the following real-valued normal form of Hopf bifurcation: where 33 U K   and the subscripts r and i denote the real and imaginary parts, respectively.Equation ( 15) has three equilibrium solutions given by: Using the normal form, the amplitudes of the limit cycle oscillations for the pitch ( A  ) and the plunge ( h A ) motions are given by: where (.) r and (.) i denote the real part and imaginary part, respectively.

Time-delay representation
Different self-excited oscillators including the van der Pol, Rayleigh, and Duffing oscillators have been employed to model self-excitations in electrical and mechanical engineering applications (Linkens, 1974;Nayfeh et al., 2003;Akhtar et al., 2009b;Abdelkefi et al., 2012f).However, the use of these self-excited oscillators to predict the behavior of aeroelastic systems is only adequate near bifurcation for supercritical type of instability (Nayfeh et al., 2013).In a recent paper, Nayfeh et al. (2013) show that the van der Pol oscillator enables a good prediction of the plunge and pitch amplitudes only near the Hopf bifurcation.On the other hand, a discrepancy was obtained between the results of their developed model and the original results of the aeroelastic system beyond the flutter onset.Furthermore, Nayfeh et al. (2013) demonstrated that the use of a Duffing and van der Pol is unsuccessful to reproduce the sudden jump to large-amplitude LCO that takes place for the subcritical instability due to numerical issues.
It has been demonstrated in the literature that time-delay ordinary differential equations can be a very good candidate to model the behaviors of many physical systems in physiology (Heiden, 1979), biology (MacDonald, 1989), and population dynamics (Kuang, 1993).In a recent paper, Nayfeh and Nayfeh (2011) used similar model to analyze the nonlinear dynamics of chatter in a lathe cutting tool.In particular, the time-delay model showed a capability to capture different types of dynamic bifurcations depending on the system configuration.Consequently, to overcome the issues of predicting the supercritical branches at higher air speed and reproducing the sudden jump to large-amplitude that takes place for the subcritical Hopf bifurcation, we first introduce a time delay to the Duffing term and remove the van der Pol term as follows: To identify the parameters defining the time delay model as given by Equation ( 22), we derive its normal form.To this end, we follow Nayfeh and Balachandran (1995) and seek a third-order approximate solution of equation ( 22) in the form: where i i Tt  . We note that the solution does not depend on the slow scale 1 T since the secular terms first appear at where / ii DT    .Substituting Equations ( 23) and (24) into Equation ( 22) and equating coefficients of like power of , we obtain where The general solution of the first-order problem is given by Abdessattar Abdelkefi and Mehdi Ghommem / Journal of Modeling, Simulation, Identification, and Control (2013) 1: 57-77 66 where cc stands for the complex conjugate of the preceding term.Substituting Equation ( 27) into Equation ( 26), we obtain where NST stands for terms that do not produce secular terms.
Eliminating the secular terms in Equation ( 28) yields the following complex-valued normal form of the Hopf bifurcation: where Introducing the polar form 29) and separating real and imaginary parts, we obtain the real-valued normal form of the Hopf bifurcation where (.) r and (.) i stand for real and imaginary parts, respectively.
We first consider the supercritical cases.We identify the values for the parameters of the time delay model.Matching the normal form yields We plot in Figs.3(a), 3(b), and 3(c) the variations of the functions 1 () H  and 2 () H  with the time delay parameter , respectively.Both of 1 () H  and 2 () H  are periodic.Thus, solving for Equation (32) may lead to many values for .We present in Table 2 the values of the model parameters for different values of .In the rest of this paper, k is considered to be equal to zero.
Abdessattar Abdelkefi and Mehdi Ghommem / Journal of Modeling, Simulation, Identification, and Control (2013) 1: 57-77   We note that the time delay model given by Equation ( 22) fails to reproduce the subcritical behavior of the aeroelastic system.So, the second attempt was to introduce a time delay to the van der Pol term and remove the Duffing term as follows Employing the method of multiple scales and following similar approach as in the previous analysis, we obtain the following normal form Abdessattar Abdelkefi and Mehdi Ghommem / Journal of Modeling, Simulation, Identification, and Control (2013) 1: 57-77 Matching the normal forms yields We plot in Figs. 5 (a), 5(b), and 5(c) the variations of the functions 1 () G  and 2 () G  with the time delay parameter τ, respectively.Inspecting Figs.5(a) and (c), we note that there are no intersection between the two curves.This indicates that building a phenomenological model as given in Equation ( 33) to reproduce limit cycle oscillations of the aeroelastic system is not feasible for We present in Table 3 the values of the model parameters for different values of τ.
Table 3 Parameter values of the time delay model.m is an integer   The time delay models given by Equations ( 22) and ( 33) show capability to reproduce the supercritical behavior of the aeroelastic system only for specific values of 2 k  .In other words, none of them was able to cover all cases.Furthermore, they fail capture the subcritical instability due to numerical issues.In attempt to have a general reduced-order model, we propose to keep the van der pol term and introduce time delay in the Duffing term; that is, The normal form of the above equation is given by We consider the supercritical cases and present in Table 4 the parameter values obtained by matching the linear and nonlinear coefficients of the normal form.In Figs. 8, 9, and 10, we plot the variations of limit cycle oscillation amplitudes for the plunge and pitch motions with the air speed U for all cases.Unlike the previous time-delay models, the present model as given by Equation (37) is able to reproduce limit cycle oscillations of the aeroelastic system for all cases.In fact, it enables an efficient way to build a reduced-order model capable of capturing both amplitudes and frequency of the limit cycle oscillation.Abdessattar Abdelkefi and Mehdi Ghommem / Journal of Modeling, Simulation, Identification, and Control (2013) 1: 57-77 74 Next, we analyze a case where the aeroelastic system undergoes subcritical instability (i.e., 0 er N  ).We follow similar approach as above and identify the model parameters.Their values are given in Table 5.We plot in Fig. 11 the variations of limit cycle oscillation amplitudes of the pitch and plunge with the air speed obtained from the original aeroelastic system and time delay model for different values of τ.A sudden jump to large-amplitude LCO takes place at the flutter onset.Of interest is that the time-delay model is able to capture this jump.However, a mismatch between the two sets of data can be observed when we move away from the Hopf bifurcation.This is expected since the system identification is based on the normal form which predicts only the unstable branch of solutions for the subcritical instability (Nayfeh, 1973).

Conclusions
In this work, we derive a reduced-order model based on self-excited oscillators to approximate the response of a nonlinear aeroelastic system.We show that time-delay differential equations present the capability to capture perfectly both supercritical and subcritical behaviors near Hopf bifurcation and predict properly limit cycle oscillations of aeroelastic systems near bifurcation and for higher air speeds.Such reduced-order model can be used to quantify the impact of each of the system's parameters on the type of instability rather than going over high-fidelity simulation, perform a rapid and reasonably accurate exploration of a large design space, and implement control strategies to exploit desirable nonlinear dynamics.
Although the analysis here has been successfully implemented for a low-order system, it can be applied to higher-dimensional systems that experience Hopf bifurcation in a similar fashion.

Fig 2 .
Fig 2. Variations of the (a) real and (b) imaginary parts of the eigenvalues as a function of the air speed.
the matrix whose columns are the eigenvectors of the matrix () new vector Y such that X=PY, Equation ( a=0 is the trivial solution.The other solutions are nontrivial.The origin is asymptotically stable for 2

Fig 3 .
Fig 3. Variations of the functions 1 () H  (solid line) and 2 () H  (dashed line) with the time delay 
Fig7.Limit cycle oscillation amplitudes of plunge and pitch motions: normal form (solid blue line), numerical results of the aeroelastic system (dashed blue line), numerical integration of the time delay equation when τ=0.349s (thick dotted green line), τ=0.810s (thick dashed black line), τ=1.733s (thick dashed brown line), and τ=4.503s (thick dashed red line).

Table 1
Parameters of the considered configurations matrix whose elements are the eigenvalues   .Based on this formulation, we note that 21 YY  and 43 YY  .Consequently, we can rewrite Equation (

Table 2
Parameter values of the time delay model.m is an integer

Table 4
Parameter values of the reduced-order model (supercritical instability)

Table 5
Parameter values of the reduced-order model (subcritical instability)