LOAD STATE OF AN AIRCRAFT WITH AN ELASTIC WING

BRNO UNIVERSITY OF TECHNOLOGY

VYSOKÉ UČENÍ TECHNICKÉ V BRNĚ

FACULTY OF MECHANICAL ENGINEERING

FAKULTA STROJNÍHO INŽENÝRSTVÍ

INSTITUTE OF AEROSPACE ENGINEERING

LETECKÝ ÚSTAV

LOAD STATE OF AN AIRCRAFT WITH AN ELASTIC WING

ZATÍŽENÍ LETOUNU PŘI VELKÝCH DEFORMACÍCH KŘÍDLA

DOCTORAL THESIS

DIZERTAČNÍ PRÁCE

AUTHOR

AUTOR PRÁCE

Ing. Pavel Schoř

SUPERVISOR

ŠKOLITEL

prof. Ing. Antonín Píštěk, CSc.

Abstrakt

Letadlo vybavené nosnými plochami je během letu vystaveno působení setr- vačných a aerodynamických sil. Aerodynamické síly působící na letadlo zav- isí na dynamickém tlaku, směru nabíhajícího proudu a tvaru letadla. Během návrhu letounu nebo kluzáku je vždy žádoucí, aby během letu bylo dosa- hováno optimálního napjatostního stavu. Dosavadní praxe používá metod které uvažují dokonale tuhou konstrukci letadla. Rozšířenou metodou je Prandtlova teorie nosné čáry. Tato teorie je ovšem použitelná pouze pro rovinná křídla. U většiny letadel lze použít předpokládat absolutně tuhé kon- stukce. Tento předpoklad ovšem nelze aplikovat na moderní kluzáky z kom- pozitní materiálů, nebo na letouny HALE (High-Altitude Long Endurance) s rozpětím křídla často překračujícím 30 m. Pokud je zatížení těchto letadel vypočítáno pomocí teorie nosné čáry, může docházet k velkým nepřesnostem.

V minulosti byly publikovány články týkající se deformace netuhého křídla.

Jako zdroj dat pro výpočet zatížení od aerodynamických sil bývá použita metoda vírové nebo dipolové mříže. Při použití této metody je křídlo mod- elováno pomocí plochých panelů umístěnými na střednicové ploše křídla.

Nevýhodou této metody jsou nepřesnosti při určování sklonu vztlakové čáry, klopivého momentu a úhlu nulového vztlaku. Tato metoda ovšem umožnuje řešit křídlo s obecnou geometrii a jeji implementace a použití je velmi jednoduché.

Další literatura uvádí řešení problému letadla s netuhým křídlem pomocí soudobých numerických metod výpočtu proudění (CFD) spojených s meto- dami strukturální analýzy, nejčastěji metodou konečných prvků. Při použití CFD je možné dosáhnou velmi přesných výsledků, avšak za cenu velmi velké náročnosti na výpočetní výkon. Během certikačního procesu kluzáku je běžné kontrolovat více než tři sta případů zatížení. Výpočet zatížení pomocí CFD je téměř nemožný, pokud není k dispozici výkoný počítačový cluster.

V této práci je navržena metoda výpočtu zatížení letadla s netuhým křídlem, založená na spojení panelové metody prvního řádu dle Katz and Plotkin, Low-Speed Aerodynamics, 2001 s metodou stukturální analýzy dle Píštěk et al., Pevnost a životnost letadel I, 1988 a Lebofsky,Numerically Generated Tangent Stiffness Matrices for Geometrically Non-Linear Struc- tures, 2013. Panelová metoda poskytuje přasná data pro výpočet zatížení křídla od vzdušných sil za předpokladu že lze dané proudění aproximovat po- mocí potenciálního proudění, Narozdíl metod založených na interakci s CFD metodami lze navrženou metodu používat i na bežném počítači.

Summary

During it’s flight, a winged aircraft is subjected to inertia and aerodynamic forces. What actually keeps the aircraft in the air is the aerodynamic force, which depends on dynamic pressure, relative wind vector and a shape of the aircraft. An ultimate goal of an airplane or sailplane designer is to design the internal structure of the aircraft so that an optimal stress state of the aircraft is achieved during its flight. Traditionally this is done by methods which assume rigid structure of the aircraft. Prandtl’s lifting line theory is widely adopted method to obtain the aerodynamic loads acting on a wing. This theory is however valid for planar wings only. The assumption of aircraft’s rigidity is valid for most aircraft, however for example modern composite gliders or High altitude, long endurance (HALE) airplanes with wingspan of 30 meters or more do not fulfill this assumption and the loads predicted by lifting-line model may become very inaccurate.

Some papers regarding deformations of an elastic wing have been pub- lished. In these articles a vortex/doublet lattice method (VLM/ DLM) is used as a source of aerodynamic loads. In the VLM/DLM methods the wing is modelled by flat panels with vertices located usually at local chamber lines.

A big disadvantage of these methods is inaccurate information about lift line slope, pitching moment and zero lift angle of attack. On the other hand, they can handle a non-planar geometry and they are easy to implement and use.

Other papers suggest a computational fluid dynamics (CFD) methods with conjunction with structural finite elements method to solve the problem of aircraft with elastic wing. The CFD methods can give very precise results, but they require much of computational power During a certification process of a sailplane, usually more than 300 load cases are evaluated. Therefore the evaluation of all load cases with CFD may become almost impossible, unless a big computer cluster is available.

This thesis proposes a methodology, how to calculate the load state of an aircraft with an elastic wing by first-order panel panel method according to Katz and Plotkin, Low-Speed Aerodynamics, 2001 coupled to a struc- tural analysis method according to Píštěk et al., Pevnost a životnost letadel I, 1988 and Lebofsky,Numerically Generated Tangent Stiffness Matrices for Geometrically Non-Linear Structures, 2013. For a flow, which is similar to a potential flow, panel methods can provide an accurate information about all characteristics required for a calculation of aerodynamic load. Unlike the CFD based methods, this method can be used on an average personal computer.

Klíčová slova

Letadlo, Letoun, Kluzák, Aerodynamika, Simulace letu, Výpočet zatížení, Elastické křídlo

Keywords

Aircraft, Aeroplane, Sailplane,Aerodynamics, Flight simulation, Load cal- culation, Elastic wing

SCHOŘ P.Zatížení letounu při velkých deformacích křídla. Brno: Vysoké učení technické v Brně, Fakulta strojního inženýrství, 2018. 65 s. Vedoucí dizertační práce prof. Ing. Antonín Pištěk, CSc..

I, PAVEL SCHOŘ, declare that this thesis titled, LOAD STATE OF AN AIRCRAFT WITH AN ELASTIC WING and the work presented in it are my own. I confirm that:

• This work was done mainly while in candidature for a research degree at Brno University of Technology.

• Where any part of this thesis has previously been submitted for a degree or any other qualification at this University or any other institution, this has been clearly stated.

• Where I have consulted the published work of others, this is always clearly attributed.

• Where I have quoted from the work of others, the source is always given. With the exception of such quotations, this thesis is entirely my own work.

• I have acknowledged all main sources of help.

Ing. Pavel Schoř

Contents

1 Introduction 4

1.1 Scope of this thesis . . . 4

1.1.1 Thesis structure . . . 4

1.2 Static aeroelasticity . . . 5

1.2.1 Motivation . . . 5

2 Governing equations 8 2.1 Domain decomposition . . . 8

2.2 Flow domain . . . 8

2.2.1 Navier-Stokes Equations . . . 8

2.2.2 Laplace’s Equation . . . 9

2.2.3 Boundary layer . . . 10

2.3 Structural domain . . . 10

2.4 Flight dynamics . . . 12

2.4.1 Coordinate systems . . . 12

2.4.2 Equations of motion . . . 12

3 Flow domain solution methods 14 3.1 Potential Flow . . . 14

3.1.1 Elements of potential flow . . . 14

3.1.2 Lifting line and vortex lattice method . . . 15

3.1.3 First order panel method . . . 16

3.1.4 Kutta condition . . . 18

3.1.5 Evaluation of panel forces . . . 20

3.2 Boundary layer . . . 21

4 Structural domain solution methods 24 4.1 Finite element method . . . 24

4.1.1 Beam element . . . 24

4.1.2 Coupling between axial and bending with torsional effects in a beam element . . . 25

4.2 Geometrical non linearity in FEM . . . 25

4.2.1 Newton-Rhapson method . . . 26

4.2.2 Numerically generated tangent stiffness matrix . . . . 27

5 Flight dynamics solution methods 28 6 Interaction between domains 29 6.1 Assumptions . . . 29

6.1.1 Beam model of a wing . . . 29

6.2 Interacting element . . . 30

6.2.1 Load . . . 30

6.2.2 Aerodynamic forces . . . 31

6.2.3 Inertia forces . . . 31

6.2.4 Deformation of the aerodynamic model . . . 31

6.3 Interaction modes . . . 32

6.3.1 Quasi static mode . . . 32

6.3.2 Unsteady mode . . . 32

7 Review of the recent state of the art of the static aeroelas- ticity 35 7.1 NASTRAN . . . 35

7.2 VLM based approach . . . 36

7.3 ASWING . . . 36

7.4 Navier-Stokes CFD based approach . . . 37

7.5 Immersed boundary method . . . 37

8 Numerical simulations 39 8.1 Validation of the FEM implementation . . . 39

8.1.1 Pure bending of an uniform beam . . . 39

8.2 Validation of the panel method . . . 39

8.2.1 NACA L5G10 . . . 40

8.3 Standard Cirrus sailplane . . . 43

8.3.1 Technical data . . . 43

8.3.2 Convergence study - airfoil FX-66-17AII-182 . . . 43

8.3.3 Standard Cirrus - integral aerodynamic data . . . 46

9 Elastic sailplane - Standard Cirrus 21m 48 9.1 Introduction and technical data . . . 48

9.2 Standard Cirrus 21m integral aerodynamic data . . . 49

9.3 Standard Cirrus 21m pull up maneuver at VA . . . 50

9.3.1 Load comparison . . . 51

9.4 Standard Cirrus 21m gust encounter at VA . . . 51

9.4.1 Load comparison . . . 52

9.5 Conclusion on Standard Cirrus 21m elastic wing . . . 53

10 Conclusion 57 10.1 Eligibility of the method . . . 57

10.1.1 Panel method . . . 57

10.1.2 Geometrically non-linear FEM . . . 57

10.2 An elastic sailplane . . . 58

10.2.1 Review of methods . . . 58

10.2.2 Summary . . . 58

10.3 Future work . . . 58

10.3.1 Improvements of the method . . . 58

10.3.2 Possible use cases . . . 59

List of Acronyms 64

List of Symbols 65

Chapter 1 Introduction

1.1 Scope of this thesis

A main goal of this thesis is to develop a methodology for computing aero- dynamic loads of aircraft with highly elastic wings. A key assumption is that the deformation progresses slowly in time. Therefore only static aeroelastic phenomena are studied. The method should be affordable for any aircraft design organization in matters of simplicity and computational costs. Al- though the name of the thesis is ”Load state of an aircraft with an elastic wing”, only modern sailplanes are evaluated in this thesis. A main reason for this is that modern sailplanes are kind aircraft, which mostly exhibit large deformations of their wings. Next, it is very easy to extend the presented method for analysis of airplane by simply adding the propulsion force into equations of motion. Also any powered airplane becomes a sailplane, once the engine stops. This method can be applied also to other kind of aircrafts, including rotorcrafts and flapping wings, but modification of the code would be necessary [4].

1.1.1 Thesis structure

The thesis is dived into five main chapters. In Chapter 1, the topic is briefly introduced by real life examples. In Chapter 2, a domain decomposition of the problem is presented. Also governing equations for each domain are introduced. Solution methods for each domain are presented in Chapter 3 to Chapter 5. Interaction between domains is explained in Chapter 6 A review of the most recent state of the art is presented in Chapter 7. Finally, numerical simulations are carried out in Chapter 8. Outcomes are discussed in Chapter 10

1.2 Static aeroelasticity

The aerodynamic forces acting on a surface of the sailplane’s wing may change its shape significantly when the stiffness of the wing is low. For instance, when witnessing a strength test of a glider wing as shown on Figure 1.1, one may naturally ask how the aerodynamic forces have changed due the deformation. According to airworthiness regulations, it is necessary to take this air load redistribution into account. The EASA CS22 regulation ex- plicitly states in paragraph CS22.301(C): ”If deflections under load would significantly change the distribution of external or internal loads, this redis- tribution must be taken to account.”[5]

Figure 1.1: A strength test of glider wing at IAE Brno University of Tech- nology

1.2.1 Motivation

A very strong motivation for detailed study of the static aeroelasticity was the event of destruction of NASA Helios HP03-2 in June 26, 2003. According to the mishap report, the Helios encountered an atmospheric turbulence, as a result of this, the wingtip deflection increased to 12 meters, as shown on Figure 1.2. After this, the Helios experienced several pitch oscillations, the last one resulted in the failure of the leading edge structure, the main spar remained undamaged as shown on Figure 1.3. In the mishap report, the NASA has described the methodology how the Helios was designed and how

the effect of the wing bending was simulated. First, the aerodynamic load of the undeformed wing was computed using VSAERO[4], then the deformed shape was computed using NASTRAN. Second, the aerodynamic load of the deformed wing was again computed using VSAERO. This load iterations were done manually. The NASA explicitly calls this approach as insufficient and sets a need for a new robust tool, which is able to compute the redistribution of the aerodynamic load on a flexible wing.

Figure 1.2: The Helios aircraft undergoing large deformation [6]

Figure 1.3: The Helios aircraft after structural failure during flight [6]

Chapter 2

Governing equations

2.1 Domain decomposition

In this thesis a problem of flight of an elastic aircraft is solved in three different domains. The forces acting on the surface of the aircraft are solved in Flow domain. The deformation of the aircraft is solved in Structural domain. When an equilibrium between the aerodynamic and internal forces in the structure is found, the resulting aerodynamic force is passed into Flight dynamics domain, where a flight condition for next time step is computed.

If a quasi static solution is desired, then only interaction between flow and structural domain is used. In this chapter only governing equations for each domain are presented, solution and interaction procedures are presented in Chapter 4.

2.2 Flow domain

2.2.1 Navier-Stokes Equations

A motion of an incompressible fluid can be described by Equation (2.1) and Equation (2.2) [1]. The first equation represent a conservation of mass and the second equation represent a conservation of momentum. These equations must be satisfied for each infinitesimally small volume of the fluid. These small volumes are called particles. The Equation (2.1) and Equation (2.2) are called Navier-Stokes equations for incompressible fluid (with constant viscosity coefficient) [1].

∇ ·u= 0 (2.1)

ρ ∂u

∂t +u· ∇u

=ρf− ∇p+µ∇²u (2.2) In Equation (2.1) and Equation (2.2):

u is velocity vector f is vector of body forces ρ is fluid density

µ is kinematic viscosity

2.2.2 Laplace’s Equation

Solving the Navier-Stokes equations provides an accurate model of the flow, but solution of these equations requires enormous effort. If the flow can be assumed as inviscid, incompressible and irrotational, then velocity potential Φcan be introduced and velocity at each point can be obtained as its gradi- ent. This flow can be described by Equation (2.4) and the complexity of the problem is reduced significantly. This flow is traditionally called a potential flow.The scalar function Φ is introduced in Equation (2.3). This function is called velocity potential and velocity−→u at each point can be obtained by its derivation:

u=∇Φ (2.3)

Substituting Equation (2.3) into Equation (2.1) gives Laplace’s Equa- tion (2.4), which also represents a continuity equation for the potential flow.

∇ ·u=∇ · ∇Φ = ∇²Φ = 0 (2.4) The momentum equation for the potential flow can be expressed as:

∇

E+p ρ +u²

2 + ∂Φ

∂t

= 0 (2.5)

In Equation (2.5), theE represents a potential of conservative body forces f:

f=−∇E (2.6)

The main advantage of the potential flow is a very low demand for com- putational power. Various methods can be used to obtain the velocity poten- tial, see Chapter 3, Section 3.1. Although the potential flow can not reflect all the physics of a real flow, the numerical solutions of Laplace’s equation

have been validated and used for decades as a primary tool for estimation of aerodynamic load [7] [8].

2.2.3 Boundary layer

In 1904 Ludwig Prandtl [9] suggested that a flow with small viscosity can be assumed as inviscid in region not very close to the surface. Close to the surface the viscous forces are dominant. The velocity changes from zero at the surface to the free stream velocity. Prandtl defined this region as Boundary layer. The thickness of the Boundary layer is defined as a distance measured from the surface where the magnitude of local velocity reaches99%

of magnitude of the free stream velocity.

Steady state two-dimensional boundary-layer flow is described by Equation (2.7) and Equation (2.8) These equations are a simplification of Equation (2.1) and Equation (2.2). A step by step derivation of Equa- tion (2.7) and Equation (2.8) can be found in literature [10], [11].

u∂u

∂x + ∂v

∂y = 0 (2.7)

u∂u

∂x +v∂v

∂y =−1 ρ

dp

dx +ν ∂²

∂y² (2.8)

2.3 Structural domain

As the aerodynamic load act on the surface of an aircraft, the aircraft struc- ture will deform. For this thesis, there is a need for a general method, which can predict a deformation of an arbitrary structure with arbitrary loads.

Such a method emerged in aerospace companies in cca. 1950’s. After several decades, it evolved to what is known as a Finite Element method (FEM).

Although the FEM is used for solving wide variety of problems such as fluid flow, heat conduction or magnetic field, it was originally developed for a structural analysis of swept wings. Therefore in this text the term FEM always refers to a solving of structural mechanics problems.

The virtual work principle The principle of virtual work states that

”The virtual work done by the body forces and surface tractions on the virtual displacement, is equal to the work done by stresses on the virtual strains”

δU =δW (2.9)

In Equation (2.9) and in Equation (2.10):

U represents the work done by stresses W represents the work done by body forces

The principle of stationary energy states that ”Of all continuous dis- placement satisfying the geometrical boundary conditions, those that satisfy the equilibrium conduction make the potential energy stationary.” For sta- ble structures the potential energy Π in this stationary condition is always minimum,

δΠ =δU −δW = 0 (2.10)

in order to satisfy Equation (2.10), the structure is discretized into large number of geometrically simpler sub-structures, called elements. Elements usually share their vertices with other elements. These vertices are called nodes. Using this discretization, the Equation (2.10) is for a steady case expressed as Equation (2.11) :

K·u=F (2.11)

Where: K is stiffness matrix of the structure u is vector of nodal dis- placements F is vector of applied nodal loads (forces and moments)

A more detailed explanation of the FEM can be found in literature, a minimal explanation is presented in Chapter 4, where a dealing with large deformations of beam elements is discussed.

2.4 Flight dynamics

2.4.1 Coordinate systems

Three different right-handed coordinate systems are used to solve the equa- tions of motion. The coordinate systems are designated as:

• E - Earth fixed, inertial

• Bf - body fixed, non-inertial

• S - body fixed, non-inertial

The coordinate systems are shown on Figure 2.1, however approach used is this thesis is different than used in literature [12], therefore the systems are briefly explained.

Earth fixed - E is the only inertial reference frame. A Flat Earth approx- imation is used - the origin is located at surface of the Earth. The Z-axis points outside the earth, X-axis and Y-axis point in arbitrary directions.

This is different to the traditional convention of NED (North East Down) reference frame [12].

Body fixed - Bf is used for solving of the Equations of motion - Equa- tion (2.12) and Equation (2.13). Its origin is always located at instant center of mass. Initially the axes are aligned with Earth fixed reference frameE. As the solution progresses in time, the Body fixed reference frame moves (and rotates) with the aircraft.

Body fixed Structural - S is used to define the geometry and position of Bfcoordinate system. It’s origin is defined by a geometry master model. A convention used here is: The X-axis points to tail, Y-axis points to the right wing, and Z-axis points upwards. As the structure deforms, origin of this axis system remains the same location. While a new position of Bf coordinate system is computed, it is always referenced to origin ofS.

2.4.2 Equations of motion

The equations of motion are solved in non-inertial reference frame Bf. The Newton’s second law can not be applied directly. A derivation of extensions of the right-hand side can be found in literature [12]. Then the Newton’s

second law in non-inertial reference frame is described by Equation (2.12) and Equation (2.13):

˙

v^bf = F^bf

m −ω˜^bfv^bf (2.12)

˙

ω^bf = J^bf−1

M−˜ω^bfJ^bfω^bf (2.13) Note that the force vectorF^bf in Equation (2.12) includes also the gravity force transformed into Body fixed reference frame.

Figure 2.1: Flight mechanics - coordinate systems

Chapter 3

Flow domain solution methods

A fluid motion is traditionally described by set of equations called Navier- Stokes equations [1]. The number of exact solution of Navier-Stokes equations is small. For practical problems, these equations are solved numerically [13].

3.1 Potential Flow

Traditionally Potential flow models are used to determine the aerodynamic load acting on a wing of an airplane. Prandtl’s Lifting line theory or Vortex lattice method are most commonly used methods. In this thesis a first order panel method [1] is used. When compared to traditional methods, the way how the load is determined is different. Both approaches are compared in following paragraphs.

3.1.1 Elements of potential flow

A uniform freestream and constant strength doublet-source panels are the only elements of potential flow used in this thesis to model the flow around an aircraft.

Uniform Freestream Velocity potential in uniform freestream at point P is expressed by Equation (3.1) [1]. Here theP₀ is arbitrary reference point

Φ(x, y, z) = Z P

P0

udx+vdy+wdz (3.1)

Constant strength source panel Velocity potential at pointP due con- stant strength source panel is expressed by Equation (3.2)[1].The cartesian

coordinates of point P, x,y,z are expressed in local coordinate system of the panel. The σ represent strength (or density) of the source.

Φ(x, y, z) = −σ 4π

Z

S

1

p(x−x₀)²+ (y−y0)²+z²dS (3.2) Constant strength doublet panel Velocity potential at point P due constant strength doublet panel is expressed by Equation (3.3)[1].The carte- sian coordinates of pointP, x,y,z are expressed in local coordinate system of the panel. The ν represents strength (or density) of the doublet.

Φ(x, y, z) = −ν 4π Z

S

z

[(x−x₀)²+ (y−y0)²+z²]³²dS (3.3) Constant strength doublet-source panel Velocity potential due con- stant strength doublet-source panel is expressed as a superposition of Equa- tion (3.3) and Equation (3.2).

3.1.2 Lifting line and vortex lattice method

The lifting line and vortex lattice methods are presented here only to demon- strate the fundamental difference between these two method and the panel method, which is in this thesis used to compute the aerodynamic loads.

Lift in VLM and lifting line theory The common denominator of Lift- ing line theory and VLM is the Kutta-Joukowski theorem, which is used to compute lift force :

L⁰ =ρvΓ⁰ (3.4)

In Equation (3.4): L⁰ is a vector of local lift force, perpendicular to vector of local velocity v. Γ is strength of a local bound vortex line segment, see below.

Lifting line theory Lifting line theory was first mathematical theory, which was capable to model a wing of finite span. It was derived by Prandtl in 1920’s. The model is shown on fig. 3.1. It consist of horse shoe vortices.

Each horse shoe vortex has three segments: Bound vortex and two wake segments. All segments have the same strength of circulation Γ. The bound vortex is attached to the wing, usually at 0.25c. Wake vortices are aligned to the free stream velocity. Wake vortices end in infinite distance behind the

wing. The boundary condition is a zero normal flow at 0.75c. The lift is determined by streght of bound vortex line segments Γ from eq. (3.4), and the local velocities vectors are given by all vortex lines segments.

Figure 3.1: Lifting line theory - vortex system

Vortex lattice method Vortex lattice method was developed by Falkner in 1943 [14]. In this method the wing is divided into segments, as shown on fig. 7.1. There is an unique horse shoe vortex attached to each wing segment.

As in Lifting line theory, the horseshoe vortex consists of bound and wake segment, where all segments of the vortex have the same circulation Γ. The control point is located at 0.75c, a zero normal flow boundary condition is prescribed at each control point. Multiple horseshoe vortices can be aranged in chordwise direction, providing both spanwise and chordwise distribution of the aerodynamic load. Unlike the lifting line theory, the vortex lattice method can be used for non-planar, swept wings and multiple lifting wings configurations. The biggest disadvantage of VLM is fact that the method requires large computational power and therefore it can by used only on computers.

3.1.3 First order panel method

An approach from VSAERO [4] is used to solve the Laplace’s equation - eq. (2.4). Assume a cross section of a wing as shown on fig. 3.3. The surface S∞ encloses the problem at infinity, surfaceS represents a body (wing) and surfaceW represents its wake. The surfaceS+W divides the domain into two regions: the external region with flow field of interest and velocity potential Φand the internal region with fictitious flow and velocity potential Φ_i. The surfaces are modeled by a doublet and source singularities.

Figure 3.2: Vortex lattice method - vortex system

Figure 3.3: Potential flow domain

Green’s theorem is applied to both outer and inner region and velocity potential Φ_P is obtained by combining both expressions. [4] The velocity potential Φ_P gives a velocity potential anywhere in the two regions. It is expressed in terms of surface integrals in terms of velocity potential and its normal derivative over the boundary surface as:

ΦP = 1 4π

Z Z

S+W+S∞

(Φ−Φi)−→n · ∇1 rdS

− 1 4π

Z Z

S+W+S∞

1 r

−

→n ·(∇Φ− ∇Φ_i)dS

(3.5)

Where r is the distance from the point P to the element dS on the sur- face, and −→n is unit normal vector of the element dS. [4] The first integral in eq. (3.5) represents disturbance potential from a surface distribution of doublets with density(Φ−Φ∞). The second integral represents contribution of sources with density −−→n ·(∇Φ− ∇Φ_i).[4]

Next, internal Dirichlet boundary condition is introduced. This boundary conditions sets the Internal flow equal to the Onset flow by:

Φ_i = Φ∞ (3.6)

Φ_P = Φ∞ (3.7)

In the panel method, the surface is discretized intonsurface panels andnw

panels in the wake region,as shown on Figure 3.4 with a constant distribution of singularities, which indicates the first order. When a linear or quadric distribution of panel’s sigularities is used, then the method is called second or higher order [1]. The Equation (3.5) is also discretized and the Dirichlet boundary condition is evaluated at centroid of each surface panel. These points are called collocation or control points. Result of this procedure is a linear equation [1] :

n

X

k=1

C_kµ_k+

nw

X

l=1

C_lµ_l+

n

X

k=1

B_kσ_k= 0 (3.8)

Where C_k, C_l are doublet influence coefficients of body and wake panels, B_k is source influence coefficient,µ_k is doublet strength of body panel, µ_l is doublet strength of wake panel,σ_k is source strength of body panel. Theσ_k is set as:

The Equation (3.8) is also used in well recoginzed and validated codes such as VSAERO [4] and PMARC [15]. The main advantage of the panels with constant strength sigularities is that lower quality surface meshes can be used, where the higher order methods require precise meshes without any irregularities [1].

σ_k =−→n_k·−→

Q∞ (3.9)

Where−n→_k is the panel normal unit vector and −→

Q∞ is the freestream velocity.

3.1.4 Kutta condition

The eq. (3.8) can be viewed in following manner: The term Pn

k=1B_kσ_k sat- isfies that the external flow does not penetrate the body surface - velocity

Figure 3.4: Panel method - body panels (red) and wake panels (grey) component normal to the panel is zero. The term Pn

k=1C_kµ_k represents a circulation, which is related to the lift of the body by Kutta-joukowski the- orem. However the value of the circulation is not known apriori. In order to specify the corculation, a Kutta condition is introduced. It states that:

”The flow leaves the sharp edge of an airfoil smoothly and the velocity there is finite” [1]. The assumption has been made that the pressure difference at the trailing edge is zero.

∆pT.E. = 0 (3.10)

If the circulation is modeled by a vortex distribution, then the eq. (3.10) is expressed as [1]:

γ_T.E. = 0 (3.11)

The circulation is then determined by setting doublet streght of wake panels as:

µ_W =µ_U−µ_L (3.12)

The Equation (3.11) and Equation (3.12) are widely used in traditional codes such as [1], [4], or [15]. However when dealing with Boundary layer interaction, it is more convenient to use formulation used in QUADPAN code [16], which uses modified eq. (3.12):

µ_W =µ_U −µ_L+ (−→

V_∞·nc_B)(−→r_{U L}·nc_B) (3.13) In the Equation (3.13), thenc_B represents an unit vector bisecting the trailing edge and −→r_{U L} represents vector which connects control points of upper and lower trailing edge panels.

The Equation (3.13) has been introduced because it allows a convenient interaction between potential flow and model of the boundary layer. See Section 3.2.

3.1.5 Evaluation of panel forces

The result of solution of Equation (3.8) is streght of source and doublet on each panel. To obtain aerodynamic forces acting on the body panel, the solution has to be post processed. Usual way is to compute panel’s local velocities by evaluating gradients of the doublet distribution. First the neigborhoud panels are determined, as show on Figure 3.5. Two surface patches are computed: l1 = P3−P1 and l2 = P4 −P2 Then the local perturbation velocities are computed:

q₁ = δµ

δl1 (3.14)

q2 = δµ

δl2 (3.15)

q₃ =σ (3.16)

The local velocity at k−th panel’s control point is then:

Q~_k =Q~_∞+~q_k (3.17)

In the Equation (3.17), the Q~_krepresentk−th panel’s local velocity,Q~∞

is freestream velocity transformed to k−th panel’s coordinate system and

~

q_k = [q₁, q₂, q₃]is local perturbation velocity computed by Equation (3.14) to Equation (3.16).

When the panel’s local velocity Q~_k is known, the panel forceF_k can be cumputed using the pressure coefficientCp:

Cp_k = 1− Q_k²

Q∞2 (3.18)

−→

F_k =Cp_k· 1

2·ρ·Q∞2·S_k· −→n_k (3.19)

Figure 3.5: Panel method - Evaluation of local panel velocity

In Equation (3.19): S_k is area of the panel and−n→_k is unit normal vector of the panel

An example of local forces acting on the surface of the wing is shown on Figure 3.6. On this figure, different color of each panel indicates different Cp value. Note that blue to green color indicates negative value ofCpand panel normals are oriented inwards the body. Yellow to red color indicates value of Cpfrom 0 to 1.

3.2 Boundary layer

An interaction between potential flow and model of a boundary layer is a common way how to extend introduce the viscous effects into the inviscid potential flow [17],[10]. Two most common methods are shown on Figure 3.7.

The top image represent the most intuitive approach, here the actual body surface is displaced by the displacement thickness of the boundary layerδ^∗.

The bottom of Figure 3.7 represents a more sophisticated approach. Here, the geometry of the body remains constant. The source term σ_k of doublet- source panels from Equation (3.9) is modified as in Equation (3.20), so that the dividing streamline is moved away from the body by theδ^∗.

σ_k =−→n_k·−→

Q_∞− ∂(U_eδ^∗)

∂s (3.20)

In Equation (3.20), the U_e represents velocity at the outer boundary of the boundary layer (Edge velocity). Thes represents distance along the stream- line.

Detailed description of the interaction scheme and the solution of bound- ary layer equations can be found in literature [10].

In this thesis, the analysis of the boundary layer is limited to quasi three dimensional analysis of the boundary layer. A following concept is used:

Figure 3.6: Panel method - wing geometry, forces

Figure 3.7: Panel method - Boundary layer interaction [17]

• An inviscid analysis is performed.

• The analyzed region is dived into span-wise sections. A stagnation point is found for each section.

• Each section is split by its stagnation point to top and bottom sub- surface.

• Boundary layer thickness and transpiration velocity is computed for each sub-surface.

• The potential flow is recomputed with updated transpiration velocities and boundary layer thicknesses.

• The procedure is repeated until integral aerodynamic characteristics stop to change between iterations.

Limitations of the boundary layer interaction used in this thesis are obvious: The solution of the boundary layer equations is valid only when there is not a pressure gradient in crosswise direction. Therefore this method can be only applied to a wing with a small dihedral and sweep angle. If there is a significant pressure gradient in crosswise direction, then full three dimensional solution of the boundary layer must be employed [11].

Chapter 4

Structural domain solution methods

4.1 Finite element method

A Lagrangian formulation of Finite element method for structural analysis was implemented. This method solves the Equation (2.11) as:

u =K⁻¹·F (4.1)

Whereuis an unknown vector of nodal displacements, later used to compute element nodal forces,K is stiffness matrix of the whole structure assembled from element stiffness matrices. The assembly process is well described in textbooks. Fis known vector of external load applied at nodes

4.1.1 Beam element

The only element implemented for structural analysis is a two node beam element. It is assumed that a beam element gives an accurate model of a section of a high aspect ratio wing - assuming one element per section. The definition of the beam element is shown on Figure 4.1.

Co-rotational framework The co-rotational framework is a useful tool for computation of element internal elastic forces. It tracks nodal position and rotation. Note that the nodal coordinate system is not aligned with the element coordinate system, as shown on Figure 4.1.

Stiffness matrix of a two node beam element used here is a standard formulation found in textbooks [18] [3]. The formulation in beam coordinate system is then:

Figure 4.1: Beam element [3]

K=









































EA

L 0 0 0 0 0 −^EA_L 0 0 0 0 0

0 ^12EIz_L3 0 0 0 ^6EIz_L2 0 −^12EIz_L3 0 0 0 ^6EIz_L2

0 0 ^12EIy_L3 0 −^6EIy_L2 0 0 0 −^12EIy_L3 0 −^6EIy_L2 0

0 0 0 ^GJ_L 0 0 0 0 0 −^GJ_L 0 0

0 0 −^6EIy_L2 0 2^2EIy_L 0 0 0 ^6EIy_L2 0 ^2EIy_L 0

0 ^6EIz_L2 0 0 0 2^2EIz_L 0 −^6EIz_L2 0 0 0 ^2EIz_L

−^EA_L 0 0 0 0 0 ^EA_L 0 0 0 0 0

0 −^12EIz_L3 0 0 0 −^6EIz_L2 0 ^12EIz_L3 0 0 0 −^6EIz_L2

0 0 −^12EIy_L3 0 ^6EIy_L2 0 0 0 ^12EIy_L3 0 ^6EIy_L2 0

0 0 0 −^GJ_L 0 0 0 0 0 ^GJ_L 0 0

0 0 −^6EIy_L2 0 ^2EIy_L 0 0 0 ^6EIy_L2 0 2^2EIy_L 0

0 ^6EIz_L2 0 0 0 ^2EIz_L 0 −^6EIz_L2 0 0 0 2^2EIz_L









































(4.2)

4.1.2 Coupling between axial and bending with tor- sional effects in a beam element

If the wing is made of composite material, then it may exhibit a coupling between axial tension and bending and torsion. This is due unbalanced composite layup [19]. These coupling effects can be easily introduced by a modification of the stiffness matrix. The resulting stiffness matrix of such an element can be found in reference [19]. These couplings were not used in this thesis - a stiffness matrix from Equation (4.2) was used.

4.2 Geometrical non linearity in FEM

As shown in Chapter 1 on Figure 1.1 and Figure 1.2, the structure of the wing may deform significantly. Using the Equation (4.1) would introduce large

errors, since it assumes small strains and deformations of the structure. For structures undergoing large deformations with small strains an equilibrium between external forcesP_extand internal elastic forcesF_int(u)is to be found, as described by Equation (4.3) [3].

P_ext =F_int(u) (4.3)

The equation is linearized by perturbing the external load and corre- sponding nodal displacements, then the Equation (4.3) becomes:

P_ext+dP=F_int(u₀+du) (4.4) Solution of Equation (4.4) is then [3] :

dP= ∂F_int

∂u du = [K_T(u)]du (4.5)

The termK_T(u)in Equation (4.5) is called Tangent stiffness matrix and is expressed as:

K_T(u) = ∂F_int

∂u =













∂Fint,1

∂u1

∂Fint,1

∂u2 . . . ^∂F_∂u^int,1

∂Fint,2 m

∂u1

∂Fint,2

∂u2 . . . ^∂F_∂u^int,2 ... ... ... ...m

∂Fint,m

∂u1

∂Fint,m

∂u2 . . . ^∂F_∂u^int,m

m













(4.6)

Note that while the stiffness matrix inEquation (4.1) is assembled for un- constrained degrees of freedom, the Tangent stiffness matrix in Equation (4.6) is assembled for all degrees of freedom in the structure.

4.2.1 Newton-Rhapson method

The Newton-Rhapson method is an incremental-iterative technique used to solve set of nonlinear equations, here the Equation (4.5). The solution process is illustrated on Figure 4.2. For the first iteration we compute the first displacement incrementδu¹₁:

δu¹₁ = [K_T(u₀)]⁻¹∆P = [K_T(u₀)]⁻¹P_ext,1 (4.7) Then the residualR¹₁ is computed. The residual is difference between the applied load and internal elastic forces.

R¹₁ =P_ext,1−F_int(u+δu¹₁) (4.8)

Then a new displacement increment is computed:

δu²₁ =

K_T(u₀+δu¹₁)−1

R¹₁ (4.9)

This process is repeated until the residual is close to zero. Then the load is incremented and whole process is repeated until the applied load is equal to the prescribed value.

Figure 4.2: Newton-Rhapson method [3]

4.2.2 Numerically generated tangent stiffness matrix

The hardest and most crucial part on geometrical non linearity is a forming of tangent stiffness matrix. It can be done analytically with much labor.

Other approach used here is to compute the derivatives from Equation (4.6) Numerically using complex step derivative method. This eliminates the an- alytical labor and increases accuracy, but also leads to longer computational times.

Chapter 5

Flight dynamics solution methods

A solution of aircraft’s equations of motion is straightforward. In this thesis, the forward Euler method is used to numerically integrate the Equation (2.12) and Equation (2.13). The integration process is described by Equation (5.1) and Equation (5.2)

v^bf_t+∆t = Z t+δt

t

˙ v^bf_t dt .

=v˙^bf_t ·∆t (5.1)

ω_t+∆t^bf = Z t+δt

t

˙ ω_t^bfdt .

=ω˙_t^bf ·∆t (5.2) As the new state vector v^bf_t+∆tω^bf_t+∆t is known, the new position of the aircraft in Earth fixed reference frame is computed.

The time step ∆t has to be selected small enough, to avoid unphysical results. Higher order integration schemes, such as Runge-Kutta method were not used due implementation difficulties, however a benefits might be gained from use of adaptive time step. The recent implementation uses only fixed size time step.

Chapter 6

Interaction between domains

6.1 Assumptions

6.1.1 Beam model of a wing

The first and most important assumption is that the wing is of a high aspect ratio and it can be modeled by using one beam (two node) element per section of the wing. Then the finite element model of the wing is assembled as shown on Figure 6.1 The internal elastic forces are evaluated at nodes with respect to element local coordinate system.

Figure 6.1: Finite element model of a wing and element coordinate systems

6.2 Interacting element

The interaction between flow and structural domain is explained on fig. 6.2.

Here a section of a wing is shown. The section consists of n aerodynamic panels, with only one element in span-wise direction. A beam element shares reference nodes of each tip cross section.

Figure 6.2: Panel method - FEM interaction

6.2.1 Load

The load is applied at the FEM nodes using Equation (6.1) and Equa- tion (6.2). The meaning and application region of f_i and r_i depends on what kind of forces is applied. Both aerodynamic and inertia forces can be applied simultaneously.

−

→f =

n

X

i=1

−

→f_i (6.1)

−

→m =

n

X

i=1

−

→f_i ×r_i (6.2)

6.2.2 Aerodynamic forces

The aerodynamic forces are always applied only at the node close to plane of symmetry. The meaning of f_i and r_i is:

fi is force on i-th panel in global coordinate system x,y,z, as shown on fig. 6.2 r_i is vector between reference node and control point of i-th panel.

6.2.3 Inertia forces

The application region of inertia forces depends on the actual design. Mul- tiple standalone masses can be used. The meaning of f_i and r_i is:

f_i is inertia force of i-th massr_i is vector between node and center of mass of i-th mass

6.2.4 Deformation of the aerodynamic model

The result of the FEM analysis is displacement vector u at each node. Co- ordinates of aerodynamic panels are updated at each load step using nodal translational and rotational displacements, according to eq. (6.3). Also the nodal coordinates are updated as explained in [3].

P T Snew= (Rx(∆φx)Ry(∆φy)Rz(∆φz))P T Sold+ ∆U (6.3) Where: ∆φ_x,∆φ_y,∆φ_z,are rotational increments in global coordinate system calculated from last step of Newton-Rhapson iteration, see bellow.

∆U is displacement increment andR_x(φ_x), R_y(φ_y), R_z(φ_z) are rotational matrices defined by eq. (6.4) to eq. (6.6) as:

R_x(∆φ_x) =





1 0 0

0 cosφ_x −sinφ_x 0 sinφx cosφx



 (6.4)

R_y(∆φ_y) =





cosφ_y 0 sinφ_y

0 1 0

−sinφ_y 0 cosφ_y



 (6.5)

R_z(∆φ_z) =





cosφz −sinφz 0 sinφ_z cosφ_z 0

0 0 1



 (6.6)

6.3 Interaction modes

6.3.1 Quasi static mode

In a quasi static mode, only the interaction between flow and structural domain is assumed. The mode is time independent. This mode is used to determine the angle of attack α, so that the wing produces enough lift that the load factor is equal to the load factor prescribed in the definition of load case.

For a rigid wing ,the task is reduced to determination of required lift coefficientC_Land interpolation of angle of attack αfrom a known lift curve.

For an elastic wing ,the task is more complicated. The aerodynamic forces acting on the wing depend on its shape. As the wing deforms, also the aerodynamic forces change. The lift is therefore a function of α, dynamic pressureq and stiffness of the wing.

The Newton-Rhapson method is used to compute the deformed shape.

Remind the Figure 4.2. Here the angle of attackαis used as the load control - increasing α increases the external load. The whole process is illustrated on Figure 6.3.

6.3.2 Unsteady mode

The unsteady mode is a time-dependent extension of the quasi static mode.

The meaning of ”increase load” from Figure 4.2 is however to advance to the next time step, respectively to next flight conditions computed from Equa- tion (5.1) and Equation (5.2). The whole process is illustrated on Figure 6.4.

Figure 6.3: Workflow diagram - quasi static mode

Figure 6.4: Workflow diagram - unsteady mode

Chapter 7

Review of the recent state of the art of the static

aeroelasticity

Some of recent research efforts are critically reviewed in this chapter. Please note that the review is made from a viewpoint of an aeronautical engineer, who is involved in a sailplane design process, and who is seeking for a right tool to predict the in-flight deformation and aerodynamic load redistribution of an sailplane with elastic wing.

7.1 NASTRAN

The NASTRAN computer code is well established standard in aero-space industry. It provides solutions to many structural analysis problems. There are several solution sequences related to both quasi-static and dynamic aero- dynamic problems.With solution sequence SOL 144 - Static Aeroelasticity [20] , the aerodynamic loads can be introduced into the FEM model. The SOL 144 provides four aero-elements for both subsonic and supersonic flow.

For subsonic flow, the doublet-lattice panels in a combination with slender bodies are available. Slender bodies are used to model the fuselage and doublet-lattice panels are used to represent lifting surfaces. As the doublet lattice method neglects all thickness effects, this solution is therefore better suitable for a high-speed jet plane, rather than a sailplane.

To couple the aerodynamic model with the structural model, an inter- polation is employed. In a brief explanation, one defines points, where the aerodynamic model interacts with the structural one and NASTRAN then redistributes the aerodynamic loads to these points using according to cho-

sen interpolation method. With the NASTRAN aero elements, also control surfaces can be modeled. With movable control surfaces, the whole aircraft can be trimmed.

7.2 VLM based approach

All methods based on Vortex lattice methods are very similar to the NAS- TRAN approach, a typical mesh is shown on Figure 7.1. Patil et al. [21]

showed the importance of nonlinear solution of the aerostatic problem. Yang et al.[22] developed a framework to trim an elastic aircraft similar to NAS- TRAN, but extended with non-linear deformations. They showed, that non- linear deformations affect the wing loads. However for all VLM based meth- ods there is always a question, if the pitching moment is introduced correctly while no actual wing surface (but only mean chamber surface) is modeled.

Figure 7.1: Vortex lattice method -mesh

7.3 ASWING

ASWING [23] is a computer code developed by Drela [23]. It is designed ”for the aerodynamic, structural, and control-response analysis of aircraft with flexible wings and fuselages of high to moderate aspect ratio.” It provides a complex framework, capable of solving both static and dynamic problems.

According to [24] and [23] the ASWING uses Lifting line theory, however the vortex lattice method is also referred as a source of aerodynamic data. Unlike Drela’s most popular code X-FOIL, the ASWING is not public available, therefore it is not possible to investigate the source of aerodynamic data.

The ASWING uses 2D sectional data to introduce drag and pitching moment of lifting surfaces,therefore is is theoretically possible to indirectly include

viscous effects. As the structural model, a non-linear beam is used. This makes possible to solve large deflection on sailplane’s wing.

7.4 Navier-Stokes CFD based approach

Many Computational Fluid Dynamics - Finite Element method (CFD-FEM) couplings are available today, mostly as a part of commercial packages like ABAQUS, or ANSYS workbench. For example Schutte et al. [25] used DLR- tau code coupled with finite element method and flight mechanics module to simulate unsteady aerodynamics of a free flying aeroelastic combat aircraft.

Solving the Navier-Stokes equations makes the CFD codes more general than simplified Euler or Potential flow codes. So one can precisely all aerodynamic properties of a sailplane, but even with the most recent hardware this proce- dure requires much of computational power and is therefore very expensive.

7.5 Immersed boundary method

The immersed boundary method [26] is an elegant solution of fluid-structure interaction problems. The flow domain consist of a fixed (regular) grid and moving boundary of the immersed object, as shown on Figure 7.2. The advan- tage is that the Navier Stokes equations - Equation (2.2) and Equation (2.2) can be solved directly on a regular grid, which reduces th computational time. Alhough the method has been used to solve real life problems [26], is is probably not suitable for solving flight of an elastic aircraft. As the Reynolds number increses, the mesh has to be denser and the computational cost increases rapidly.

Figure 7.2: Immersed boundary method - grid and boundary

Chapter 8

Numerical simulations

8.1 Validation of the FEM implementation

8.1.1 Pure bending of an uniform beam

A case of a pure bending of an uniform beam was selected for a validation of the FEM. The beam is loaded by a bending moment at it’s right tip.

The left tip is fully clamped. As suggested by other references [27], the beam should curl into a complete circle, when the value of bending moment reachesM = ^2πEI_L^zz A deformed shape of the beam is shown of Figure 8.1 for a sequence of four load steps, with maximum value of the bending moment M = ^2πEI_L^zz. As predicted, the beam curled into a complete circle.

8.2 Validation of the panel method

Two validation cases are presented in this section to verify if the 3D panel method is a suitable tool to derive the aerodynamic load and to estimate its errors, when compared to the experimental data. The boundary layer has not been employed for any case in this section.

NACA 0012 experimental data Chord-wise pressure distribution was compared to experimental data provided by Gregory [28] for Re = 3·10⁶. Since the data provided by Gregory are two dimensional, following wing geometry was used to avoid three dimensional effects:

• planform : rectangle

• span b= 100m

Figure 8.1: Pure bending of a beam

• chord c= 1m

Local pressure coefficients are plotted in plane of symmetry of the wing, results for angle of attack α = 15^◦ from panel method are compared to experimental data on Figure 8.2. Good agreement with the experiment is achieved, when comparing integrals Rc

0 cp(x) dx, the difference is 5.3%

8.2.1 NACA L5G10

The panel method was validated against experiments in wind tunnel. NACA wartime report L5G10 [29] was used as a source of experimental data. Ge- ometry of a wing used in the experiment is shown on Figure 8.3. In the panel method, the wing was span-wisely divided into 48 sections. Each section had 50 panels on the surface of the wing and 10 wake panels behind the wing Integral aerodynamic data Comparison of lift curve slope is shown on Figure 8.4

Span-wise lift distribution Comparison of local normal force coefficients C_N is shown on Figure 8.5.

Conclusion on NACA L5G10 Both comparisons indicates good agree- ment with the experimental data for M = 0.2. Based on the comparison of

Figure 8.2: Panel method validation - comparison of CP coefficients the numerical and experimental load, the 3D panel method can assumed as a suitable tool to provide aerodynamic loads.

Figure 8.3: Panel method validation - wing geometry [29]

Figure 8.4: Panel method validation - comparison of lift curves

Figure 8.5: Panel method validation - comparison of C_N coefficients

8.3 Standard Cirrus sailplane

A Standard Cirrus sailplane [30] is used for evaluation of the panel method on full aircraft configuration. This type of sailplane has been selected due availability of the geometry and flight test data [31].

8.3.1 Technical data

The geometry of the Standard Cirrus sailplane is shown on Figure 8.6 and briefly described in Section 8.3.1. The reference point for measuring the center of gravity position in Section 8.3.1 is located aft leading edge of wing’s root rib.

Figure 8.6: Standard Cirrus sailplane - 3D view

8.3.2 Convergence study - airfoil FX-66-17AII-182

A convergence study was performed to verify the sensitivity of the bound- ary layer solution to different chord-wise paneling. An airfoil FX-66-17AII-182 is analyzed in this study. This airfoil is used on wing tip section of the Stan- dard Cirrus. Three diffrefent geometries are compared to the experimental

Wing span 15m

Wing area 10m²

Root airfoil FX S02-196 Tip airfoil FX-66-17AII-182

Mean chord 0.65m

Dihedral 3^◦

Length 6.4m

Height 1.32m

Aspect ratio 22.5

Weight 330kg

CG position 250mm aft ref. pt.*

L

D max 35.2

Manoeuvring Speed 169kmh⁻¹ Never Exceed Speed 220kmh⁻¹

Table 8.1: Standard Cirrus technical data

data for Re = 4e6. On Figure 8.7 and Figure 8.8, the geometries are desig- nated as: ”Inviscid”, ”n190”, ’n60’. For all geometries a logarithmic spanwise paneling (denser mesh at leading edge) was used. Because of the 3D panel method, the airfoil was simulated using a wing of a very high aspect ratio and wing span. For this study, the wingspan was100m, while the chord was 1m. Only one section per span was used.

n60 A very rough mesh with 120 panels per wing section. It was proven that this rough mesh is not satisfactory for modeling boundary layer effects, due large errors in lift curve slope and over-predicted drag as shown on Figure 8.7 and Figure 8.8.

n190 A fine mesh with 380 panels per wing section. As shown on Fig- ure 8.7 and Figure 8.8, this fine mesh is suitable for boundary layer analysis, the predicted lift and drag match close to the experimental data.

Inviscid An optimized mesh with 220 panels per wing section. This mesh could be used for boundary layer analysis, but the drag would be over predicted. This mesh is later used for inviscid analyses. Note that the drag shown Figure 8.8 is a result of discretization errors rather than true viscous drag.

Figure 8.7: Airfoil FX-66-17AII-182 - lift coefficient

Figure 8.8: Airfoil FX-66-17AII-182 - drag polar

Conclusion on the convergence study The convergence study has shown that for the airfoil FX-66-17AII-182 operating atRe= 2e6and moderate an- gles of attack the boundary layer effect could be ignored with relatively small errors. With respect to that, the boundary layer was not used for later anal- yses to save the computational time. Usually to find a converged solution

(a) Comparison of speed polars (b) Lift curves

(c) Aerodynamic polars (d) Pitching moment Figure 8.9: Standard Cirrus - integral aerodynamic data

of boundary layer, it took five iterations, therefore five times increase in computational time.

8.3.3 Standard Cirrus - integral aerodynamic data

Next, the geometry of a whole Standard Cirrus was analyzed using and com- pared to the flight test data [31] as shown on Figure 8.9. Note that the deflection of the elevator was set to zero, therefore the pitch was untrimmed.

This is the most probable explanation for over predicted drag at high speeds shown on Figure 8.9d. There are two designation used on Figure 8.9, the

0BLCG250⁰ indicates a fine mesh (see Section 8.3.2 above) coupled with boundary layer analysis. The center of mass is located 250mm aft the da- tum reference point. The⁰CG250coarse⁰ indicates the optimized mesh from Section 8.3.2, center of mass is located250mmaft the datum reference point.

Conclusion The comparison of computed integral aerodynamic data to experiment indicates a reasonable agreement in regions, where the pitching moment is close to zero. It it reasonable to omit the boundary layer analysis, since the inviscid data is close to viscous data at angle of attack bellow7.5^◦.