Numerical model for the simulation of fixed wings aeroelastic response

TECHNICAL PAPERS

Numerical model for the simulation of fixed wings aeroelastic response

G. R. Benini; E. M. Belo; F. D. Marques

Escola de Engenharia de São Carlos Universidade de São Paulo NPA - Núcleo de Pesquisas em Aeronáutica Lab. de Aeroelasticidade, Dinâmica de Vôo e Controle Av. Trabalhador Sancarlense, 400 13566-590 São Carlos, SP. Brazil

ABSTRACT

A numerical model for the simulation of fixed wings aeroelastic response is presented. The methodology used in the work is to treat the aerodynamics and the structural dynamics separately and then couple them in the equations of motion. The dynamic characterization of the wing structure is done by the finite element method and the equations of motion are written in modal coordinates. The unsteady aerodynamic loads are predicted using the vortex lattice method. The exchange of information between the aerodynamic and structural meshes is done by the surface splines interpolation scheme, and the equations of motion are solved iteratively in the time domain, employing a predictor-corrector method. Numerical simulations are performed for a prototype aircraft wing. The aeroelastic response is represented by time histories of the modal coordinates for different airspeeds, and the flutter occurrence is verified when the time histories diverge (i.e. the amplitudes keep growing). Fast Fourier Transforms of these time histories show the coupling of frequencies typical of the flutter phenomenon.

Keywords: Aeroelasticity, flutter, vortex lattice method

Introduction

Aeroelasticity deals with a class of fluid-structure interaction problems where airflows are involved. Aeroelastic problems in aircraft have been observed since the earliest days of flight. The most dangerous aeroelastic problem is the flutter phenomenon, which can be defined as a catastrophic dynamic instability. The improvements in aircraft performance and the utilization of lighter and consequently more flexible structures increase the susceptibility to aeroelastic problems such as flutter. In this context, the development of aeroelastic numerical models for utilization in the aircraft design phase and the conception of mechanisms to suppress the aeroelastic problems have become very important and have received special attention from the research community.

The numerical models for aeroelastic analysis can be divided in two vast categories, depending on if they treat the problem in the frequency or in the time domain. The solutions based on the frequency domain are the classical ones, but are valid only in the stability boundary, i.e., can be employed only for the prediction of critical conditions for flutter occurrence. In contrast, the solutions based on the time domain allow the determination of the structure aeroelastic response for any flight conditions, and have the additional advantages of allowing the inclusion of non-linear effects and the design of control systems for flutter suppression.

The aim of this work is to present a numerical model for aeroelastic analysis in time domain. The methodology is to treat the aerodynamics and the structural dynamics separately and then couple them in the equations of motion. The dynamic characterization of the structure is done by the finite element method (Zienkiewics, 1986) and the unsteady aerodynamic loads are predicted using the vortex lattice method (Katz and Plotkin, 1991). More details about the present model can be found in Benini (2002) and similar models were presented by Strganac and Mook (1990) and Preidikman and Mook (2000).

The vortex lattice method is based on potential aerodynamics and, as presented in this work, has limited applicability in predicting the actual loads on modern aircraft configurations operating at the transonic range. However the method has the advantage of a low computational effort when compared with modern CFD codes and can be utilized as an effective tool in the design of control systems for flutter suppression (Luton and Mook, 1993; Hall et al., 2000).

Nomenclature

a = aerodynamic influence coefficient

G,= trasformation coordinates matrices

I = identity matrix

K = stiffness matrix

K_m = modal stiffness matrix

L = aerodynamic loads vector

m = number of aerodynamic panels

M = mass matrix

M_m = modal mass matrix

n = normal vector

n = number of modes chosen for the simulation

N = number of structural degrees of freedom

p = static pressure

r = position vector in the Biot-Savart law

S = panel area

t = time

v_m = velocity of the wing motion

v_w = velocity induced by the wake

V

= induced velocity (or local velocity in the Bernoulli equation)

V¥ = airspeed

W = factor for local truncation error calculation

x, = physical displacements

= physical velocities

= physical accelerations

y = modal displacements and velocities

y₁ = modal displacements

y₂ = modal velocities

Greek Symbols

a = angle of attack

d_x= virtual displacements

D b = length of vortex ring in the spanwise direction

Dc = length of vortex ring in the chordwise direction

Dt = time interval

f = potential velocity

j = mode shape

F, = modal matrices

G = circulation

h = modal displacements

= modal velocities

= modal accelerations

r = air density

w = natural frequency

w

²

= diagonal matrix containing the squared natural frequencies

Subscripts

1,2 = refer to the ending points of a vortex segments

a = refer to the aerodynamic points

i,j = aerodynamic panel identifiers

K = refers to a vortex ring

l = refers to the lower side of the panel

L = refers to a control point

r = refers to a particular mode of vibration

s = refer to the structural points

TE = trailing edge

u = refers to the upper side of the panel

Superscripts

C = refers to the corrected solution

P = refers to the predicted solution

T = transpose of a matrix

Dynamic and Structural Model

The wing structural response is assumed to be linear and without internal damping. The equation of motion for the structure discretized in N degrees of freedom (DOF) is shown in Eq. (1), where M and K are N x N matrices, representing the mass and stiffness properties, and x(t), (t) and L(x,, t) are N x 1 vectors, representing the displacements, accelerations and external (aerodynamic) forces. It is important to note that the aerodynamic damping is provided implicitly by the aerodynamic forces.

The process of discretization is performed using the finite element method, and the eingenvalue problem (Eq. 2) provides the natural frequencies w_rand the mode shapes j_r.

The mode shapes can be arranged in a matrix according to Eq. (3). This matrix is named modal matrix and is used as a coordinate transformation matrix, according to Eq. (4), where h(t)represents the structural displacements in a modal domain and can be interpreted as a vector of coefficients which determines the influence of each mode shape in the physical structural response (Meirovitch, 1986).

Since the matrix F is constant with time, Eq. (5) can be written as

Substituting Eqs. (4) and (5) in Eq. (1) and pre-multiplying both sides by F^T (transpose of F) yields Eq. (6), where M_m = F^TMF and K_m = F^TKF are named modal mass and modal stiffness matrices, respectively.

Due to the orthogonality properties of the mode shapes, one can prove that the matrices M_m and K_m are diagonal matrices. In addition, it is possible to normalize the eingenvectors in a form that M_m = I, and then the division of both sides of Eq. (6) by the matrix M_myields Eq. (7), where w² is a diagonal matrix containing the squared natural frequencies.

In order to simplify the solution of Eq. (7), it is useful to consider only a few natural modes to describe the structural response. This is done truncating the summation in Eq. (4). In fact, only a few modes are necessary to obtain a solution with good precision (Meirovitch, 1986).

Finite Element Model

The wing employed in the aeroelastic simulations was designed by the EESC/USP Team for the 2000 SAE International AeroDesign^® East Competition, Florida, USA. The wing has a rectangular planform, with a 0.22 meter chord and a 3 meters span. The structural scheme consists of a Styrofoam^TM core, a Kevlar^TM shell in +/- 45º orientation and two unidirectional (UD) Carbon fiber spars, as shown in Fig. 1 (EESC/USP Team, 2000). The carbon fiber spars have constant width and vary from one layer at the wing tip up to three layers at the wing root. A finite element (FE) model of the semi-wing was generated using the software ANSYS^®. Two types of elements were employed: a quadrilateral laminated shell element for the composites, and a tetrahedral solid element for the Styrofoam^TM core. The FE mesh over the semi-wing is shown in Fig. 2, together with the schemes of the mode shapes selected to represent the dynamic response.

Aerodynamic Model: The Vortex Lattice Method

The vortex lattice method consists of distributing plane vortex singularities over the wing and over the wake. The plane vortex singularities satisfy the Laplace equation and when combined with the uniform stream can simulate incompressible and potential flows around the wing.

To implement the method, the wing is represented by a lifting surface without thickness and discretized in quadrilateral elements, called panels. A vortex ring is associated with each panel, being the leading segment of each vortex ring placed on the panel quarter chord line and its control point placed at the center of the three-quarter chord line. The wing discretization scheme is shown in Fig.3. To guarantee that the flow streamlines pass over the lifting surface, it is necessary to satisfy the boundary condition of zero normal velocity on the wing surface. This boundary condition is applied at the control points and results in the correct values for the vortex singularities (represented by the circulation G).

The boundary condition in each panel can be expressed according to Eq. (8), where the gradient of the potential velocity f corresponds to the perturbed velocities induced by the wing vortex singularities, v_m corresponds to the velocity of the wing motion (the freestream velocity relative to the wing plus the velocities of the wing structural deformations), v_w corresponds to the velocities induced by the wake, and n is the normal vector.

The velocity (V) induced by each straight vortex segment, extending from point 1 to point 2, at an arbitrary point P, is given using the Biot-Savart law (Eq. 9), where r₁ and r₂ are the vectors that define the position of point P in relation to the points 1 and 2.

It is important to note that the value of the circulation G is still not known in Eq. (9). So, only the values of the other terms will be calculated. This is done by assuming G=1. The velocity induced by each vortex ring at a point P is obtained adding the results obtained with Eq. (9) for the four corresponding vortex segments.

The velocity V_KL will be referred as the velocity induced by the vortex ring L on the control point K. Applying the zero normal velocity boundary condition at the control point K=1, Eq. (10) can be written, being the circulation in each vortex ring the unknowns and m the number of panels used in the wing aerodynamic discretization. Based on Eq. (10) the called influence coefficients (a_KL = V_KL. n_K ) are defined. Writing Eq. (10) as a function of the influence coefficients for each of the m control points and passing v_m and v_w to the right-hand side (RHS) of the equation, the linear system represented in Eq. (11) is obtained.

The evaluation of v_m consists of two steps: 1) the freestream velocity is obtained moving the wing in the aft direction, and 2) the velocities of the structural deformations are obtained solving the equation of motion (Eq. 7). The velocities induced by the wake (v_w vector) are obtained employing the Biot-Savart law (Eq. 9). It is important to consider that a portion of the wake is generated at each time interval, according to Fig.4. The circulation values of the last vortex rings generated are the same as those of the trailing edge vortex rings, to satisfy the three-dimensional Kutta condition. Thus, at each time interval new vortex rings are generated and the corresponding values of circulation are found. The value of circulation of each wake vortex ring remains the same during all the simulation time. In the present simulation, the wake rollup was not considered, so the wake is parallel to the freestream velocity plane.

The solution of the linear system of Eq. (11) provides the circulation values for the wing vortex rings, which will be employed for the aerodynamic loads calculation. The unsteady Bernoulli equation for each panel is written in Eq. (12), where p is the static pressure and the subscripts u and l refer to the upper and lower sides of the panel.

The last two terms in Eq. (12) refer to the unsteady case. The difference between them is given by Eq. (13), obtained from the definition of circulation (Katz and Plotkin, 1991).

If , Eq. (12) is analogous to the classical Bernoulli equation for the steady case, and the first two terms can be determined with the aid of the Kutta-Joukowski theorem (Eq. 14), where V_¥ is the freestream velocity, a is the local angle of attack, D_b is the length of the panel in the spanwise direction and S is the panel area.

Substituting Eqs. (13) and (14) in Eq. (12), the normal force in each panel can be computed and supplied as input to the equation of motion (Eq. 7). It is important to emphasize that the values of G in the above equations are given by G_i,j for the wing leading edge panels, and by (G_i,j - G_i-1,j) for the other panels (Katz and Plotkin, 1991).

Coupling Between Aerodynamic and Structural Meshes

Because the mode shapes and the aerodynamic forces are calculated using distinct meshes, it is necessary to convert the aerodynamic forces to the structural points and to supply the structural displacements to the aerodynamic points. The scheme chosen to exchange this information between the meshes is based on the fact that the aerodynamic and structural points are related through a coordinate transformation matrix, according to Eq. (15), where the subscript a refers to the aerodynamic points, the subscript s refers to the structural points and the matrix G is the transformation matrix. The same matrix can be employed to write the mode shapes in terms of the aerodynamic points, according to Eq. (16).

To write the aerodynamic forces in terms of the structural points, it is necessary to guarantee that the virtual works done by the forces are the same in both meshes, according to Eq. (17), where and are the virtual displacements. As these displacements are arbitraries, it is possible to write Eq. (18).

Substituting Eq.(18) into Eq. (7) and making use of Eq. (16) yields Eq. (19), which represents the conversion of forces between the two meshes.

It is important to note that the aerodynamic forces are applied at the control points, represented here by the vector x_a, and that the structural displacements must be supplied at the vortex rings corner points, grouped in the vector . So, to perform the conversion of displacements between the meshes, it is important to define a different transformation matrix, which will be called . Pre-multiplying both sides of Eq. (4) by the matrix , the modal displacements provided by the solution of Eq. (19) are converted to the physical displacements in the vortex rings corner points, according to Eq. (20).

The matrices F_a and _a can be obtained by the method of interpolation using surface splines (Harder and Desmarais, 1972). The method is based on the solution of the bending equation of a flat plate and have some limitations. Firstly, the aerodynamic and structural points should be projected in the wing plane, and secondly, only the out-of-plane displacements are interpolated. To implement the method, some structural points (65 in this work) were selected in the FE mesh to reproduce the structural modes. A schematic picture showing the superposition of these 65 structural points and a 52 panels aerodynamic mesh can be seen in Fig. 5.

Integration of the Equation of Motion

The equation of motion represented in Eq. (19) consists of a system of n second order ordinary differential equations (ODE), being n the number of modes chosen to represent the structural response. To facilitate the integration, this system will be re-written as a system of 2n first order ODE, utilizing Eqs. (21) and (22), where y₁ and y₂ are n x 1 vectors representing the modal displacements and modal velocities, respectively. The vectors y₁ and y₂ can be grouped together according to Eq. (23), where the vector y has dimension 2n x 1 and represents the modal displacements and velocities.

Differentiating Eqs. (21) and (22) and considering Eq. (19) yields Eqs. (24) and (25).

Eqs. (24) and (25) can be rearranged in a unique equation according to Eq. (26), where the vector is obtained in the same way than the vector y (Eq. 23) and the vector RHS represents the corresponding right hand side terms of Eqs. (24) and (25).

The 2n ODE of Eq. (26) will be solved employing a predictor-corrector method, implemented in PECLE form (Lambert, 1991). The letter P refers to the prediction of the vector y, the letter E to the evaluation of the vector RHS, the letter C to the correction of the vector y values, and the letter L to the evaluation of the local truncation error. The family of Adams-Bashforth methods is employed as the predictor equations, and the family of Adams-Moulton methods is employed as the corrector equations. The local truncation error is given by the Milne's estimate and modifies the corrected values, according to Eq. (27), where the superscripts P and C refer to the predicted and corrected values and the factor W depends on the order of the predictor-corrector equations. For the first time interval (i=1) only the E step is performed. For the second time interval (i=2) only the PE steps are performed. Then, for i=3 an up, all the PECLE steps are performed. The predictor-corrector equations for each time interval, together with the factor W value, are shown in Tab. 1.

Results and Discussion

Aeroelastic responses of the prototype aircraft wing model are shown in Figs. 6, 7 and 8, for different airspeeds. The aeroelastic responses are represented by time histories of the modal coordinates. The fast Fourier transforms (FFTs) of the aeroelastic responses show the oscillatory frequencies contained in each response. All of the FFT curves have been normalized to make their peak amplitudes equal to one.

The structural and aerodynamic meshes are the same as those represented in Fig. 5 and no camber was considered for the wing. The initial angle of attack is 5 degrees and the air density corresponds to the sea level condition. The modal displacements and velocities are set to zero as the initial conditions for the solution of Eq. (26) and the initial perturbation is due to a step input in the freestream velocity, in the aerodynamic model.

In Fig. 6 (V_¥ = 30 m/s) the aerodynamic damping is strongly evidenced for all modes and the responses reach the wing static equilibrium position. The FFTs show the peaks clearly separated, with exception of the 2^nd and 3^rd modes. In Fig. 7 (V_¥ = 70 m/s), all the modes are beginning to be excited by a common frequency (around 65 Hz), but the response is still stable. In Fig. 8 (V_¥ = 100 m/s), this common frequency dominates the motion of all modes and the response is clearly unstable. The general form of the curves seems to be in agreement with similar simulations (Preidikman and Mook, 2000), but because the structures are completely different it is impossible to have a quantitative comparison.

The wing tip displacements for the leading and trailing edges, and the wing tip angle of attack time histories are shown in Figs. 9 to 12 for the subcritical (V_¥ = 30 m/s) and supercritical (V_¥ = 100 m/s) conditions. These curves were obtained converting the modal displacements to physical displacements by using Eq. (4).

The influence of the time interval (Dt) in the numerical simulations is emphasized in Fig. 13, where the aeroelastic response of the 5^th mode for V_¥ = 10 m/s is shown for two different time intervals. It is observed that a numerical instability can occur depending on the value chosen for Dt. This indicates that one should be very careful when analysing the results, to avoid a numerical instability to be interpreted as a physical instability (flutter). In this work the simulations for each speed were performed for different time intervals. An initial value for Dt was assumed and then this value was decreased until the point where no difference was noticed between consecutive simulations. This procedure however can be very time consuming.

As the dependence of the time interval for the numerical stability is a characteristic of the predictor-corrector methods (Hughes, 1987), one of the options to improve the presented numerical model is the automation of the time interval choice, or the use of a different numerical method to integrate the equation of motion. It is important to point that predictor-corrector methods were also employed by Strganac and Mook (1990) and Preidikman and Mook (2000) in similar aeroelastic models, but no mention was made about problems with numerical instability.

Conclusions

The aeroelastic response of a prototype aircraft wing is obtained for different airspeeds. The results presented in this paper were obtained for three particular airspeeds, corresponding to a stable condition (30 m/s), an almost unstable condition (70 m/s) and a clearly unstable condition (100 m/s). It should be pointed that the exact flutter speed is anywhere between 70 and 100 m/s and that the only way to discover it is performing simulations for small increments in the airspeed. In this sense, if one is interested just in computing the flutter speed it is better to use a frequency domain model, which provides the flutter speed directly. However the time domain simulations can provide valuable additional information, such as the vibration levels suffered by the structure, and can be used as an effective tool to design control systems for flutter suppression or gust alleviation.

It was shown that the precision of the results depends strongly on the time interval choice, what can make the application of the numerical model very time consuming. Procedures to automate the time interval choice or to avoid this dependence would be of great value.

As a final remark, it is important to emphasize that the obtained responses are valid only for incompressible and potential flows.

Acknowledgements

The authors acknowledge the financial support of the São Paulo State Research Foundation, FAPESP (grant 00/00361-3).

References

Benini, G.R., 2002, "Numerical Model for the Simulation of the Aeroelastic Response of Fixed Wings" (In Portuguese), MSc. Dissertation, Escola de Engenharia de São Carlos, Universidade de São Paulo, 80p.

EESC/USP Team, 2000, "Technical Report - AeroDesign East", Deland.

Hall B.D., Mook D.T., Nayfeh A.H., Preidikman S., 2000, "A Novel Strategy for Suppressing the Flutter Oscillations of Aircraft Wings", 38^th AIAA Aerospace Sciences Meeting & Exhibit, Reno.

Harder R.L., Desmarais R.N., 1972, "Interpolation Using Surface Splines", Journal of Aircraft, Vol.9, No.2, pp.189-191.

Hughes, T.J.R., 1987, "The Finite Element Method: Linear Static and Dynamic Finite Element Analysis", Prentice-Hall.

Katz J., Plotkin A., 1991, "Low-Speed Aerodynamics: From Wing Theory to Panel Methods", McGraw-Hill.

Lambert J.D., 1991, "Numerical Methods for Ordinary Differential Systems: The Initial Value Problem", John Wiley & Sons.

Luton J.A., Mook D.T., 1993, " Numerical Simulations of Flutter and its Suppression by Active Control", AIAA Journal, Vol.31, No.12, pp.2312-2319.

Meirovitch L., 1986, "Elements of Vibration Analysis", McGraw-Hill, 2^nd ed.

Preidikman S., Mook D.T., 2000, "Time-Domain Simulations of Linear and Nonlinear Aeroelastic Behavior", Journal of Vibration and Control, Vol.6, pp.135-1175.

Strganac T.W., Mook D.T., 1990, "Numerical Model of Unsteady Subsonic Aeroelastic Behavior", AIAA Journal, Vol.28, No.5, pp.903-909.

Zienkiewics O.C., 1986, "The Finite Element Method", McGraw-Hill, 3^rd ed.

Paper accepted January, 2004. Technical Editor: Clóvis Raimundo Maliska

Part of this paper has been published at the XXIII Latin Ibero American Congress on Computational Methods for Engineering, 24-26 June 2002, Giulianova, TE, Italy, under the title: "Numerical Model for Aeroelastic Analysis in Time Domain".

Publication Dates

• Publication in this collection
  12 Aug 2004
• Date of issue
  June 2004

History

• Received
  Jan 2004