DIPLOMARBEIT

Aeroelastic Analysis of a Glider Aircraft

Verfasst unter der Leitung von:Professor Ulf Ringertz (PhD)Department of Aeronautical and Vehicle EngineeringRoyal Institute of TechnologySE 100 44 Stockholm, Sweden

An der TU Wien beurteilt von:O.Univ.Prof. Dipl.-Ing. Dr.techn. Franz RammerstorferInstitut fur Leichtbau und Struktur BiomechanikTechnische Universitat Wien

1040 Wien, Osterreich

ausgefuhrt zum Zwecke der Erlangung des akademischen Grades einesDiplom-Ingenieurs und eigereicht an der Technischen Universitat Wien

Fakultat fur Maschinenwesen und Betriebswissenschaftenvon

Florian TOTHe0525414

Johann Teufelgasse 181230 Wien, Osterreich

Wien, am 14.6.2010

1

Die approbierte Originalversion dieser Diplom-/Masterarbeit ist an der Hauptbibliothek der Technischen Universität Wien aufgestellt (http://www.ub.tuwien.ac.at). The approved original version of this diploma or master thesis is available at the main library of the Vienna University of Technology (http://www.ub.tuwien.ac.at/englweb/).

Aeroelastic Analysis of a Glider Aircraft

Florian TOTHDepartment of Aeronautical and Vehicle Engineering

Royal Institute of TechnologySE 100 44 Stockholm, Sweden

June 3, 2010

Abstract

A simulation model for static aeroelastic phenomena is created, using a genericsailplane for the FAI1 18m-class as an example. The sailplane is modeled structurallyby a Nastran finite element model. Material testing of selected carbon-fiber-epoxy composites is conducted to verify the assumed material data. Structuraleigenmodes are used as a modal basis to couple aerodynamic and structural modelin a modal subspace approach. The aerodynamic forces are computed using anunsteady potential flow model. Equations of motion for a steady pull-up maneuverare formulated for the elastic aircraft. The obtained results of the static aeroelasticcalculation are compared with results for the rigid aircraft. Finally, the applicabilityof the used modeling tools for structural optimization with respect to aeroelasticeffects is discussed.

1 Introduction

Grass waving in the wind or leaves fluttering – aeroelasticity is engulfing us in our dailylives. If an elastic body is exposed to a fluid flow, it will deform under the imposed load.This in turn will change the aerodynamic forces. All physical phenomena arising in theinteraction between aerodynamic forces and elastic restoring forces are incorporated inthe field of aeroelasticity [2]. This term, however, is not completely descriptive, since fordynamic aeroelastical phenomena, inertial forces have an important impact too [3].

The field of aeroelasticity can be separated into two main areas. The first one involvesall problems, in which time is not an independent variable. Therefore it is termed staticaeroelasticity [2]. Steady aerodynamics can be used, and often it is the aim to investigatethe arising deformations in certain flight conditions, which affect aircraft performance [7].

1Federation Aeronautique Internationale

2

Control surface efficiency of elastic aircraft changes with flight speed. For example, de-flecting the aileron down will lead to increased camber of the airfoil, and therefore toan increased lift resulting in the desired rolling moment. On the other hand, deflect-ing the aileron downwards imposes a twisting moment on the wing, deflecting the wingnose downwards. Thus the effective angle of attack, and respectively, the effective rollingmoment is decreased [2]. As the aerodynamic twisting moment imposed by the ailerondeflection is growing with the airspeed, the desired effect of a control surface deflection isreduced. Finally, there will be a speed at which the rolling moment becomes zero and aaileron deflection has no effect at all. This is called the aileron reversal speed.

Another static aeroelastic phenomenon is divergence, which arises if a deformation incre-ment leads to a larger increment in aerodynamic forces than in elastic restoring forces[3]. Considering only linear theory, the aerodynamic forces are linearly dependent on thedynamic pressure and the wing’s deformation. Considering the frequently used two di-mensional example of the ’typical section’, it is shown that, if the aerodynamic center ofa wing’s cross section is located in front of it’s twist center, there will be a certain speedat which the aerodynamic forces and the elastic restoring forces grow equally fast [3, 2, 7].This speed is called the divergence speed. If it is exceeded, the aerodynamic forces cannotbe balanced by the elastic restoring forces any more, resulting in a stability problem.

If inertial forces have to be taken into account to describe certain phenomena, dynamicaeroelasticity has to be used. One of the most far reaching aeroelastic effects is perhapsflutter [3]. Modern high-speed aircraft are subject to many different kinds of flutter,involving separated flow, periodic breakaway of the flow and time lag effects between theaerodynamic forces and the motion of the wing [3]. Classical flutter is usually associatedwith unsteady potential flow and a coupling with structural eigenmodes, inducing self-amplified oscillations, which can lead to catastrophic structural failure within seconds[3, 2]. Therefore the stability of the dynamic system including unsteady aerodynamicloading has to be investigated to find the flutter speed.

It is to be observed, that exceeding the flutter speed or the divergence speed of an aircraftwill lead to fatal failure, as these two phenomena are stability problems. On the otherhand, control reversal won’t lead to structural failure, but possibly to serious handlingproblems for the pilot [7].

Although analytical solutions for certain aeroelastic problems exist, they are more suitableas guides in showing the relative influence of the various parameters [1]. Approximatemethods are frequently used to predict the aeroelastic behavior of given vehicles. Mostof them can be broken down into two steps. First the infinite-degree-of-freedom system,mathematically prescribing the problem, is approximated by a finite set of degrees offreedom, leading to a system of simultaneous equations [1]. Commonly the finite elementmethod is used to discretize the structure and boundary element methods are used toderive a discrete formulation of the aerodynamic equations. In a second step to find anapproximate solution the equations have to be solved simultaneously. The discretizedequations can be algebraic or differential with time as independent variable [1]. Staticaeroelastic problems often reduce to linear algebraic systems, solvable by the various

3

techniques of matrix inversions, or linear eigenvalue problems, where the roots of thecharacteristic polynomial can be found exactly. The set of differential equations describingdynamic aeroelastic problems can be Laplace-transformed, if it may be assumed thatexcitation at a certain frequency will lead to response at this frequency only. The resultingcomplex eigenvalue problem can then be used to study the stability of the system in thefrequency domain [1].

Modern high performance sailplanes are constructed using the principle of light weightdesign. Their structure, often consisting of composite materials, is highly flexible. In orderto achieve the maximum performance, the elastic deformations of the aircraft during flighthave to be taken into account, meaning that an static aeroelastic analysis is necessary.Furthermore all aircraft have to be well clear of stability problems like flutter or divergencein the specified flight envelope [13].

The purpose of the project is to investigate the suitability of the tools currently available atKTH’s Department of Aeronautical and Vehicle Engineering for the purpose of aeroelastictailoring of sailplanes. It is to be investigated to which extent structural optimization isfeasible. Therefore a simulation model of a gilder is developed. The design is based onan existing aircraft, which is currently available at the Department of Aeronautical andVehicle Engineering, and was analyzed in respect to geometry, airfoil shape, and material.A selection of the most frequently used carbon fiber composite materials is tested to obtainaccurate properties. Finally, a static aeroelastic trim analysis is carried out as a test casefor other aeroelastic investigations.

2 Sailplane and modeling tools

The investigated aircraft is a generic 18m glider. It is assumed to be manufactured formFRP2, a significant part of it being carbon fiber. A cargo bay located behind the mainwing should be usable to hold a computer to record measurement data during flight.Additionally is should be possible to install a power unit in the cargo bay. In the area ofthe intersection of main wing and fuselage, additional stiffness is provided by a interiorfiber reinforced composite structure. Furthermore additional bulkheads, and frames inlongitudinal direction are used to form closed cross sections in the area of the cargo bay.This is necessary to provide sufficient torsional stiffness.

To achieve a reasonable pressure distribution in the aerodynamic analysis, airfoil shapeshave to be carefully designed. This is necessary, because the pressure distribution is verysensitive to the local curvature. Based on measurements on an existing glider, suitableairfoils for main wing and stabilizers are chosen from a catalog [8].

The glider’s exterior geometry is modeled in sumo. This surface modeling tool, especiallydesigned for aircraft structures, was developed at KTH’s Department of Aeronautical andVehicle Engineering. The software should not be regarded as a full-fledged CAD–system3,

2fiber reinforced plastic3Computer Aided Design

4

but rather a easy-to-use sketchpad, highly specialized towards aircraft structures [14].The fuselage surface is modeled by defining cross sections. Their horizontal and verticaldimensions can be modified interactively using the fuselage’s side- and top view. Wingsurfaces are modeled by specifying airfoils at two or more desired locations. The glider’sfinished surface geometry can then be saved in the IGES–format4, to be further processed.

The structural modeling is done using NX6, a PLM-software5 by Siemens. It was mainlyused to add stiffeners and bulkheads to the exterior surfaces. Furthermore cutouts forthe canopy, gear and cargo bay, and a fairing between tailboom and vertical stabilizerwere added. On the finished geometry a 2D-mesh is generated, including the definition ofmaterial properties. Beam elements are added in certain areas, and the complete mesh isexported to a Nastran input-file6.

A normal modes solution is then used as input for the potential flow solver dwfs, developedby KTH’s Department of Aeronautical and Vehicle Engineering, to compute the aerody-namic forces. Using MATLAB, the linear equation system, modeling the deformableaircraft in a steady flight condition, is solved. The solution is then visualized using scope,a visualization tool to display the results of dwfs. Additionally it is capable of animatingtime-domain flight trajectories of elastic aircraft [14].

3 The finite element model

Different levels of detail can be applied when modeling an aircraft structure. The simplestpossible deformable model would feature beam elements for wings, fuselage and tailplane.This kind of model is a strongly reduced representation of the structure, and therefore alarge amount of heuristic approximations has to be made. Without experimental data,for example from ground vibration testing, obtaining accurate results is difficult [6]. Atearly design stages, without prototypes built, experimental validation is impossible.

Another possible modeling approach is a three dimensional shell model of the aircraft.Here it is possible to model the entire structure with just minor simplifications. Thedifficulty then arises in defining the material properties of the shell elements, and inmodeling the connections between different aircraft parts. However, using this method, itis comparatively easy to reach an accurate representation of the structural stiffness.

The finite element model is build in NX6, using Nastran as solution environment.Figure 1 shows the finite element mesh. The total model has roughly 786,000 degreesof freedom. It consists of about 133,000 shell elements, whereof more than 99.5% arequadrilateral. It also contains 1,200 beam elements. Around 450 connection-elementswere used to define constraints between the different parts of the model.

4Initial Graphics Exchange Specification5product lifecycle management6to improve readability Nastran is used to refer to NX Nastran by Siemens PLM Software.

5

Figure 1: Finite element mesh

3.1 Element types

The primary used element is CQUAD4, a quadrilateral element based on Mindlin-Reissner shell theory, able to represent in-plane bending and transverse shear behavior[9]. It is chosen over the higher order CQUAD8 element due to several reasons. In thefirst place, it performs better on doubly-curved surfaces [9], which appear throughout theaircraft structure. Secondly, a consistent coupling between second order shell elementsand beam elements used as stiffeners is not possible as Nastran just offers first orderbeam elements.

CTRIA3 elements, the triangular analog of CQUAD4, are used as rarely as possible, asthey may exhibit excessive stiffness [9]. Furthermore they show a rather poor accuracy[11]. In some cases, for example for mesh transitions, it is, however, advantageous to usetriangular elements in addition to quadrilateral elements.

3.2 Beam stiffeners

Local reinforcements, for example at the canopy cutouts edges, are conveniently modeledusing beam elements. Furthermore reinforcements by additional unidirectional fibers incertain areas are modeled by beam elements. This saves considerable work, arising by theneed of different shell properties in these areas. The beam element used is CBEAM, alinear beam element with two nodes, including extension, torsion, bending in two perpen-dicular planes and the associated shear [9]. It is chosen over the CBAR element, becauseit offers additional capabilities like distributed torsional mass moment of inertia.

6

3.3 Connections

The modeling of connections between the different parts of the aircraft is one of the moredifficult parts of the creation of the finite element model. This is due to the need toapply large simplifications to the real conditions. In Nastran this can be done by usingso called R–Type Elements, which impose fixed constraints between the points to whichthey are connected [9]. Thus, they are mathematically equivalent to one or more multi-point constraints, which are commonly used by various finite element codes to modelconnections.

In Nastran the assembly of the stiffness matrix differs from the classical approach,where the individual element stiffness matrices, commonly obtained by integrating theshape functions in the local element coordinate system, are transformed into a singleglobal coordinate system and assembled into the global stiffness matrix. In Nastranit is possible to specify individual displacement coordinate systems for arbitrary nodes.The stiffness matrix is then assembled similarly, but the element stiffness matrices aretransformed into the specified displacement coordinate systems. This approach facilitatesa variety of elegant modeling techniques [10], for example for the modeling of hinges withaxes in arbitrary direction.

The main concern when modeling the connections is to achieve a most accurate repre-sentation of the structure’s global behavior. It is important to consider, that the shellstructure has low bending stiffness. Therefore loads must be introduced at locations wheretwo shells intersect. To minimize the effects of point-loads, introduced by constraints cou-pling two single nodes, RBE3-spider elements can be used. These elements couple asingle dependent node to a set of independent nodes by calculating a weighted average ofthe corresponding degrees of freedom.

3.4 Material properties

The quality of any FEM model is directly dependent on the accuracy of the used materialproperties. Especially when dealing with composite materials, the acquisition of materialdata is difficult. The properties of the finished laminate depend on the used ply materi-als, resin type, material orientation, and the manufacturing process, including layup andcuring. The properties of most ply materials can be assumed to be orthotropic. Whenapplying classical lamination theory, the orthotropic shell properties can be computedaccording to the defined ply layup. In Nastran this is done by the entry PCOMP todefine the property of a shell element. The layup in terms of ply thickness, ply orientationin respect to the element reference orientation, and orthotropic material properties of thesingle plies have to be defined.

A lamina cannot withstand high stresses in any other direction than that of the fibers.Therefore the assumption of a plane stress state is not merely an idealization of reality, butinstead a practical and achievable objective of how to use laminates [4]. The stress-strainrelations for a orthotropic ply in plane stress state (σ3 = τ23 = τ31 = 0) are expressed by

7

the compliance matrix and can be written as ε1ε2γ12

=

1E1

−ν12E1

0

−ν21E2

1E2

0

0 0 1G12

σ1σ2τ12

. (1)

The engineering constants to be determined are the in-plane Young’s Moduli in 1- and2-direction E1 and E2, the in-plane shear modulus G12, and the in-plane Poisson’s ratioν12.

To experimentally determine the engineering constants, the unidirectional tension test isapplied. Here a thin flat strip of material, having a constant rectangular cross section, ismounted in the grips of a mechanical testing machine, and monotonically loaded in tensionwhile recording the force and the occurring strains [15]. The tests were conducted at theLightweight Structures Lab at KTH. The test specimen were prepared according to ASTMstandards7. Figure 2 shows obtained stress-strain plots for different carbon-fiber-epoxycomposites and the gage region of the corresponding test specimen.

Figure 2: stress-strain plots for unidirectional material loaded in fiber direction and wovenfabric loaded in warp-, weft- and 45-direction

Conducting unidirectional tension tests in the lamina’s 1-direction yields by definition

σ1 =P

A, E1 =

σ1ε1

and ν12 = −ε2ε1. (2)

7American Society for Testing and Materials

8

Here P is the applied load, and A the specimen’s cross section. Note, that ε1 and ε2 haveto be measured simultaneously. To determine E2 the same experiment is conducted in the2-direction. The obtained stiffness properties should then satisfy the reciprocal relation

ν12E1

=ν21E2

. (3)

To determine the Poisson’s ratio a third experiment, loading the specimen at 45 to the 1-direction, is necessary. Ex, the effective modulus in measurement direction is determinedby

σx =P

Aand Ex =

σxεx. (4)

Using the modulus transformation relation for lamina loaded in arbitrary direction yields

G12 =1

4Ex− 1

E1− 1

E2+ 2ν12

E1

. (5)

It is to be noted, that coupling between the normal stress σx and the shear strain γxyoccurs at the 45off-axis test. Therefore the ends of the lamina must be free to deform,or the specimen has to be sufficiently long, to ensure that no Saint-Venant effects arise inthe gage section [4].

4 Modal subspace formulation

Under the assumption of small, elastic deformations, linear equations of motion can beused. It is to be noted that effects like contact, arising for example in the connectionbetween wings and fuselage, cannot be treated properly by a linear model. These ef-fects, however, are local, and can be approximated through a linear model. The discreteequations of motion can be written as

Mu+Ku = fa(t,u, u, . . .), (6)

where M is the mass matrix, relating the inertial forces to the accelerations in the struc-tural degrees of freedom u. The stiffness matrix is denoted by K, and fa is the vector ofaerodynamic forces.

The general modal subspace method can be used to drastically reduce the number ofdegrees of freedom. A modal basis for the system (6) is found by solving the Eigenvalueproblem [

K − ω2M]z = 0. (7)

The eigenvectors zj corresponding to the eigenfrequencies ωj solving (7) are used toexpress the displacement vector as

u =n∑j=1

ηjzj, (8)

9

where ηj are the modal coordinates. Let n be the size of the original system. By justusing an appropriate subset of m nodal coordinates an adequate approximation of u isachieved, but the number of equations of motion is considerably reduced by n−m. Thisapproximation can be written in matrix form as

u ≈ Zη. (9)

The columns of the n × m matrix Z are the Eigenvectors, corresponding to the usedmodal coordinates contained in η. If the eigenvectors are mass-normalized, meaning thatzTjMzj = 1 and zTi Mzj = 0 for all i 6= j, the generalized mass matrix becomes theidentity matrix and the generalized stiffness matrix Ω is diagonal containing the squaresof the eigenfrequencies ωj. The motion of the structure is now restricted to the modalsub domain with m degrees of freedom. Inserting (9) into (6) and premultiplying withZT leads to the modal equations of motion

Iη + Ωη = f η(t,u, u, . . .). (10)

The load vector is projected into the modal subspace by f η = ZTfa. Note, that (10)expresses the balance between inertial and elastic restoring forces on the left hand side,and aerodynamic forces on the right. The left hand side of the equations of motion isdecoupled by the transformation into the modal subspace, however coupling in the modalequations of motion still arises by the aerodynamic forces on the right hand side of 10.

The selection of modes used to compose the modal subspace is crucial for the accuracyof the model. In dynamic analyses usually the first few lowest frequency eigenmodes areused, which in general leads to good results. In the case of static analyses, the usedmodes have to be able to represent the static deformations of the structure sufficientlywell. When investigating the global behavior of the structure, the first few eigenmodesare sufficient.

4.1 Transformation of modeshapes

The eigenvectors are calculated for the unconstrained system, including mechanisms forthe six independent control surfaces. This means that the stiffness matrix is singular. Intotal 12 of the system’s eigenfrequencies are zero. The 12 corresponding eigenvectors aredetermined by the Nastran solver and, in general, don’t span the space in a convenientway. A linear combination of the original nodal coordinates ηj can express the desired newmodeshapes, being rigid body displacement in coordinate directions, rotations around thecoordinate axes and relative rotations of the control surfaces around their hinge lines. The

10

linear combination is expressed by

Tx1 Tx2 . . . Tx12...

......

Rz1 Rz2 . . . Rz12

δ11 δ12 . . . δ112...

......

δ61 δ62 . . . δ612

0

0

1. . .

1

︸ ︷︷ ︸

A

η1η2...η12η13...ηm

︸ ︷︷ ︸η

=

tx...rzδ1...δ6η13...ηm

︸ ︷︷ ︸η′

. (11)

The coefficients of the transformation matrix A are found in the first 12 eigenmodes. Tx1is the 1st eigenmode’s x-translation component of a node with stiffness in all directionsattached, close to the center of gravity and not located on a control surface. The remainingcoefficients in the first row of A are corresponding to eigenmode 2 to 12. Likewise, for thecoefficients in row 2 to 6, the y- and z-displacements, and respectively, the rotations aroundx-, y- and z-axis of the same node have to be inserted. The coefficients in row 7 to 12,corresponding to relative control surface deflections, are obtained by taking characteristicdifferences between a displacement- or rotation component of the above used node, andthe corresponding component of a node located on the control surfaces. Modes 13 to n areelastic modes and do not have to be transformed. Nevertheless the same method couldbe applied to filter out control surface movements from these modes.

5 Aerodynamic model

A linear flow model is sufficient, as incremental motion from an equilibrium state is in-vestigated. This is consistent with the use of a linear finite element model. As sailplanesoperate at comparatively low speeds, the flow can be treated as incompressible. At rea-sonably high Reynold’s numbers viscous effects are confined to the boundary layer andwake [12], making a potential flow model accurate. The aerodynamic model is based onthe following assumptions:

• linear aerodynamics are used,

• the fluid is incompressible,

• inviscid as well as irrotational flow,

• the structure conducts simple harmonic oscillations with frequency ω in the gener-alized coordinates.

11

Therefore an unsteady potential flow model is useful. It is to be noted, that due to thelimitations of potential flow, an accurate prediction of the drag force is impossible. Asno performance calculations are conducted in the scope of this paper, this is no seriouslimitation.

To compute the generalized aerodynamic forces, the potential flow solver dwfs, developedat KTH’s Department of Aeronautical and Vehicle Engineering, is used. The aerodynamicmesh is created on the wetted surface of the aircraft, using the software sumo. Thetriangular mesh can be created automatically, using just a few input parameters. Figure 3shows the coarse aerodynamic mesh as well as the pressure distribution for the heavingmode. The central problem when using a modal subspace formulation, is the interpolation

Figure 3: Aerodynamic mesh and pressure distribution

of the aerodynamic mesh to the structural eigenmodes [6]. This is done automaticallyby the flow solver. The calculated vector of the generalized aerodynamic forces can beexpressed in the form

f q = f 0 +1

2ρU2Q0q. (12)

In addition to the generalized modal forces of the elastic modes f η = [Q1, . . . , Qm]T , itcontains the total aerodynamic forces and moments in the body reference frame in theform of

f q = [Fx, Fy, Fz,Mx,My,Mz, Q1, . . . , Qm]T . (13)

The forces are calculated by integrating the pressure distribution, and projecting it ontothe modeshapes [6], and respectively, the coordinate directions. It is to be noted, that thetotal aerodynamic moments Mx, My and Mz are calculated with respect to an arbitraryaerodynamic reference point. To represent the equivalent moments L, M and N in respectto the aircrafts center of mass, the columns in (13) have to be linearly combined, usingthe vector from the aircraft’s center of mass to the aerodynamic reference point ~r, in the

12

manner of LMN

=

Mx

My

Mz

+ ~r ×

FxFyFz

. (14)

In (12) ρ is the fluid density, and U the airspeed. The first term f 0 in (12) represents thestatic aerodynamic forces at a certain reference state q0. The static aerodynamic matrixQ0 contains the aerodynamic derivatives. It accounts for the change of the aerodynamicforces due to deviations form the reference state. The deviations are described by thestate vector q containing the modal coordinates ηi.

The flow solver dfwSolve works in the frequency domain to compute the complex valuedunsteady aerodynamic matrix Qk for a given reduced frequency k = ωb

U. The reduced

frequency, also known as Strouhal’s number, can be interpreted as a nondimensionalmeasure of the unsteadiness of the flow [7]. It is calculated using the semi-chord b and theairspeed U . The real part of the unsteady aerodynamic matrix, evaluated at low reducedfrequency, can be used as an approximation for the static aerodynamic matrix.

For the trim problem, investigated in the following, the state vector has to be augmentedto contain the changes in angle of attack and pitch rate. It then has the form

q = [∆α, q, η1, . . . , ηm]T . (15)

Due to the assumption of linear aerodynamics, the corresponding aerodynamic deriva-tives can be calculated by finite differences, using two independent static solutions. Thederivative vectors, forming the first two columns in the aerodynamic matrix, are therefore

1

2ρU2Cα =

fα1− f 0

α1 − α0

and1

2ρU2Cq =

f q1 − f 0

q1 − q0. (16)

The static aerodynamic forces fα1are evaluated at a state differing from the reference

state just in the angle of attack by ∆α = α1−α0. The same holds for f q1 in terms of thepitch rate.

6 The trim problem

Trimming an aircraft means to find a configuration in which steady-state flight is possible.In the scope of this paper, the 3-dimensional trim problem is restricted to the problem oflongitudinal control, just allowing a pitching motion. Furthermore wind is not includedin the model. The general equations of unsteady motion, stated in Etkin and Reid [5],are reduced to two dimensional space leading to

X +mg sin θ = m (u+ qw) ,Z −mg cos θ = m (w − qu) ,

M = Jy q,

θ = q,xE = u cos θ + w sin θ,zE = −u sin θ + w cos θ.

(17)

13

The total mass of the aircraft is denoted by m. The aerodynamic forces X and Z areorientated in the body frames x- and z-direction, and M is the aerodynamic momentaround the y-axis, with respect to the aircrafts center of mass, as shown in Figure 4.

zE

Z

XM

xE

~vθ

R z

x

g

1

Figure 4: Forces and coordinate systems

The relative velocity ~v = (u, v, w) of the body fixed frame against the inertial framedefines the angle of attack α by

tanα =w

v. (18)

It is to be noted, that the components of ~v are commonly defined in the moving bodyframe. In the current, two dimensional longitudinal case, the angular velocity q, prescrib-ing the pitch rate, is the derivative of the Euler angle θ, which defines the orientation ofthe moving body frame in respect to the fixed inertial (’earth’) frame.

The aerodynamic forces X and Z and the aerodynamic moment M are dependent on thetotal airspeed U = ‖~v‖, the change in the angle of attack ∆α and on the pitch rate q,as well as on the deformation of the structure. Using the above described aerodynamicmodel, taking the aircraft’s center of mass position into account, and omitting the columnscorresponding to lateral motion, the aerodynamic forces are formulated asZM

f η

= f +1

2ρU2

[cα cq Q

] ∆αqη

. (19)

Here f represents the static aerodynamic forces. The columns of the aerodynamic matrixcorresponding to changes in the angle of attack and the pitch rate are cα and cq. Thefollowing columns in the aerodynamic matrix, corresponding to structural modes, arecontained in Q.

In a steady-state flight condition u = w = q = 0 holds. Setting the load factor nz, definedas the ratio between normal force Z and the aircrafts weight mg, constant, and assumingθ = 0, the trim conditions are obtained from the equations of motion (17) giving

Z = nzmg, (20)

M = 0, (21)

q =(1− nz)g

u. (22)

14

Assuming that transient effects can be neglected, the acceleration term in the modalequations of motion (10) can be dropped. Inserting the aerodynamic forces (19) andadding the trim conditions (20) and (21) leads to the linear equation system0

0ΩA−1

− 12ρU2

[cα cq QA−1

]︸ ︷︷ ︸

S

∆αqη′

︸ ︷︷ ︸

x

= f +

−lfmg00

︸ ︷︷ ︸

a

. (23)

Note, that the transformation relation (11) has been used to enable the locking of modaldegrees of freedom related to control surface movements. The rigid body modes containedin η′ have to be constrained, fixing the aircraft in its body coordinate system. Additionally,all control surface deflections apart from the elevator deflection will be constrained.

The trim condition (22) requires the system vector x to be divided into a prescribed partxp and a free part xf . The system (23) then takes the form[

Spp SpfSfp Sff

] [xpxf

]=

[apaf

]. (24)

By adding the constrained modes as zero to the prescribed part, the constrains can beenforced. Additionally the antisymmetric modes can be added as zeros too, as they arenot relevant in the investigated case. The system to be solved is then

Sffxf = af − Sfpxp. (25)

7 Results

In addition to the aeroelastic equations of motion (23), the equations of motion for the two-degree-of-freedom system of the rigid aircraft were solved. These equations, containingthe change of angle of attack and the elevator deflection as variables, can be obtainedfrom (23) by constraining all elastic degrees of freedom.

7.1 Aeroelastic deformations during level flight

The case of level flight can be regarded as a special case of the previously described pull-upmaneuver. Setting the load factor nz = 1 means that lift and gravity forces are equal andthe pitch rate is zero.

Figure 5 shows the angle of attack α and the elevator deflection δ for various airspeedsfor the elastic, as well as for the rigid aircraft. The angle of attack necessary to maintainlevel flight is larger for the elastic aircraft. The effect becomes more pronounced withhigher airspeeds, because elastic deformations increase. This means that the elastic con-figuration’s performance is worse, as higher angles of attack increase lift-induced drag.

15

Figure 5: Angle of attack and elevator deflection during level flight

For the rigid aircraft, angle of attack and elevator deflection are both zero at an airspeedof 145km/h. Likely the aircraft, which served as a base for the model, was designed to flyat this speed with zero angle of attack of the fuselage and neutral elevator position. Ingeneral, a significant elevator-down command has to be applied to achieve level flight athigh speeds.

The wingtip deflection with respect to the fuselage in vertical direction increases almostlinearly with airspeed, as displayed together with other characteristic elastic deformationsof the glider during level flight in Figure 6. The twist angle of the wing around it’s

Figure 6: Characteristic deformations during level flight

spanwise axis is evaluated at 86% spanwise distance. It is positive for airspeeds greaterthan 104km/h, meaning that the leading edge is twisted upwards. For lower airspeedsthis effect is inversed. The bending of the tailboom is characterized by the angle between

16

the vertical stabilizer’s root and the wing root. It is lower than zero for all investigatedairspeeds, meaning that the end of the tailboom is bent downwards with respect to thewing root.

7.2 Aeroelastic behavior in a steady pull-up maneuver

A steady pull-up maneuver with load factor nz = 4 was simulated for a range of airspeeds.Figure 7 shows the deformed aircraft at two different airspeeds. Note that, for clarity,the deformations are scaled by factor 5. Wingtip deflections reach from 426mm, at anairspeed of 150km/h, to 526mm, at an airspeed of 280km/h (see Figure 9).

Figure 7: Deformations during a pull-up maneuver with load factor 4

The simulation shows that the maneuver is hardly possible for airspeeds lower than150km/h, as the necessary angle of attack and the elevator deflection have to be un-realistically large. This can be seen in Figure 8, showing the angle off attack and theelevator deflection during a steady pull-up maneuver with load factor 4.

Figure 8: Angle of attack and elevator deflection during a steady pull-up with load factor 4

17

Trying to fly such a maneuver would result in a stall, due to flow separation rapidly de-creasing lift. It has to be noted, that the simulation results for large angles of attack cannotbe regarded as valid, as linear theory is used. Furthermore effects like flow separation arenot included in the applied aerodynamic model. Observing the elastic configuration showsan interesting phenomenon. Contrary to the first expectation, the elevator is deflecteddownwards in pull-up maneuvers at airspeeds greater than 150km/h. This phenomenoncan easily be explained by a comparison with the elevator deflection at level flight (seeFigure 5), which is smaller (meaning deflected ’more’ up) at corresponding airspeeds.

The characteristic deformations of the glider during a pull-up maneuver with load factornz = 4 can be observed form Figure 9. When comparing the characteristic deformationsof the wing at different load factors (Figure 6 and 9) it is apparent, that higher wingloading leads to a twisting moment deflecting the leading edge downwards.

Figure 9: Characteristic deformations during a pull-up with load factor 4

7.3 Modeling effort and suitability for structural optimization

A large part of the overall modeling effort is consumed by the creation of the finite elementmodel. As the gliders exterior geometry is created in sumo and then imported into NX6,changes require considerable work. The interior geometry, like bulkheads and stiffeners,has to be recreated even if no changes in their position or size was intended. This alsomeans that the geometry has to be re-meshed.

Changes in the interior geometry usually just require re-meshing. The creation of a meshof good quality on the whole structure takes about 30 hours. To minimize the effort ofmeshing due to geometry changes, the structure is separated in sub-parts, which haveto be connected by constraints. The number of sub-parts shouldn’t be too high as theconstraints have to be redefined, if the mesh is changed in a sub-part.

As the geometry of the wing is likely to be changed a lot when conducting aeroelastic

18

tailoring, a C++ program is used to create a structured mesh of quadrilateral elementson the wing structure, using the surface data defined in sumo as an input. Control surfacepositions can be defined and stringers can be added in the desired locations. Furthermorethe program allows to specify different material properties in predefined areas. The createdmesh can then be imported into NX6 and has to be coupled with the fuselage.

If studies on the impact of material properties or material orientation are desired, this canbe accomplished by directly editing the Nastran input-file. A new modal solution hasto be computed after every change in the finite element model. The computational timeof about 10 minutes on a standard PC is negligible, compared to the modeling effort.

Changes in the exterior geometry also require a re-creation of the aerodynamic mesh. Thisis done automatically with SUMO. As the flow solver uses the structural modes as input,the aerodynamic matrix has to be recomputed if changes in the FE model are made. Thecomputational time for a solution on a coarse mesh is about 2 minutes on a standard PC.

Due to the application of a modal approach the system of aeroelastic equations to be solvedin MATLAB is very small, involving 10 to 20 degrees of freedom, depending on how manymodes are used. Computation time is of no concern in this area. The post-processing ofthe obtained results, and the time required to understand the complex relations betweencertain input parameters and results, should not be underestimated.

8 Conclusions

The static aeroelastic deformations calculated at the investigated flight conditions werereasonable. Wingtip deflection of 526mm was established at pull-up maneuver with loadfactor 4 and an airspeed of 280km/h. The elastic model was compared to a rigid model,which showed the same characteristic behavior in the investigated flight cases. The quan-titative differences between an elastic and a rigid model, however, were found to be sig-nificant, showing the importance of aeroelastic analysis.

The developed model is usable for structural optimization in respect to its aeroelastic be-havior, but structural changes in fuselage or tailplane require considerable effort. Changesin the wing geometry, material properties or orientation on the other hand, can be appliedrather quickly. It should, however, be considered to extend the functionality of the devel-oped C++ program to enable automatic coupling of the wing structure with the fuselagestructure.

References

[1] R. M. Bisplinghoff, H. Ashley. Principles of Aeroelasticity. John Wiley and Sons,Inc., 1962.

19

[2] H. W. Forsching. Grundlagen der Aeroelastik. Springer Verlag, Berlin ISBN 3-540-06540-7, 1974.

[3] R. L. Bisplinghoff, H. Ashley, and R. L. Halfman. Aeroelasticity. Dover Publivations,1955.

[4] Robert M. Jones. Mechanics of Composite Materials, Second Edition. Taylor &Francis, Inc., ISBN 1-56032-712-X, 1999.

[5] B. Etkin and L.D. Reid. Dynamics of Flight John Wiley and Sons, 1996.

[6] D. Eller. On an Efficient Method for Time-Domain Computational Aeroelasticity.Department of Aeronautical and Vehicle Engineering, KTH, TRITA/AVE 2006:1,December 2005.

[7] D. Borglund, D. Eller. Aeroelasticity of Slender Wing Structures in Low-Speed Air-flow. Division of Flight Dynamics, KTH, Lecture Notes, 2009.

[8] D. Althaus. Niedriggeschwindigkeitsprofile. Vieweg, ISBN 3-528-03820-9, 1996.

[9] Siemens PLM Software Inc. Nastran Help Library. Nastran Element Library Refer-ence, NX 6.0.4, September 2009.

[10] Siemens PLM Software Inc. Nastran User’s Guide. Nastran Element Library Refer-ence, NX 6.0.4, September 2009.

[11] Richard H. MacNeal. Perspective on finite elements for shell analysis. Finite Elementsin Analysis and Design, 30(3):175–186, August 1998.

[12] J. Katz and A. Plotkin. Low Speed Aerodynamics. Cambridge University Press,second edition, 2001.

[13] European Aviation Safety Agency. Certification Specifications for Sailplanes andPowered Sailplanes. EASA, CS–22/2, March 5th, 2009.

[14] http://www.larosterna.com/, May 31st, 2010.

[15] American Society for Testing and Materials. Standard Test Method for TensileProperties of Polymer Matrix Composite Materials. ASTM International, D 3039/D3039M–08, 2008

20