Development of Longitudinal Flight Dynamics

AIAA AVIATION Forum June 25-29, 2018, Atlanta, Georgia 2018 Multidisciplinary Analysis and Optimization Conference

10.2514/6.2018-3425

Development of a Flight Dynamics Analysis Framework for Inclusion Into the MDAO process Rikin Gupta∗ , Nathan J. Love † , Rakesh K. Kapania‡ Virginia Polytechnic Institute and State University, Blacksburg VA, 24060

and David K. Schmidt§

Downloaded by VIRGINIA TECH on November 7, 2018 | http://arc.aiaa.org | DOI: 10.2514/6.2018-3425

D.K. Schmidt & Associates, Monument, CO 80132 A methodology for modeling the flight dynamics of flexible aircraft is presented. This methodology is ultimately to be incorporated into a multidisciplinary analysis and optimization (MDAO) framework. The modeling methodology is based on the mean-axis formulation of the equations of motion, and the embedded aeroelastic model is based on the use of non-dimensional stability derivatives, as used in modeling rigid vehicles. The methodology is presented and then applied to model a flexible flying-wing research drone for verification purposes. The results were shown to agree very well with previously published results, thus verifying that the methodology was implemented properly. It was then applied to model a larger flexible flying-wing vehicle, which has control surfaces along the entire trailing edge of the wing. The results indicated that this vehicle exhibited bodyfreedom flutter, occurring at 74.5 kts, with a flutter frequency of 4.01 Hz, and that the flutter mechanism involves the interactions between the rigid-body short-period (pitch) mode and the first symmetric aeroelastic mode. This first aeroelastic mode exhibits both bending and twist deformations, due to wing sweep. Below the flutter speed, the flight dynamics exhibit significant aeroelastic contributions to the pitch-rate responses of the center body from the various combinations of flap deflections, and that the eigenvalues of the dynamics of the flexible vehicle differ significantly from those of the rigid vehicle. Also, a reduced-order model of the dynamics that retains the effects of static-elastic deformations, shows that the aerodynamic coefficients adjusted for the static-elastic effects differ significantly from those obtained from a rigid model. A controllability/observability metric is offered for determining the relative importance of the system modes in the input-output dynamics of the vehicle. The resulting metrics indicate, for example, that the participation of the first aeroelastic mode in the center-body pitch-rate response is very significant, and that this contribution is maximized by using the symmetric outboard flap defections. In contrast, this contribution is minimized by using a pair of more-inboard flaps, which are located near the node line of the first symmetric vibration mode. It is concluded that incorporating this modeling methodology into an MDAO framework should be useful for capturing the effects of flight dynamics and control in optimizing vehicle designs.

Nomenclature m: Mass of the aircraft Iyy : Pitching moment of inertia about y-axis Ixx : Rolling moment of inertia about x-axis Izz : Yawing moment of inertia about z-axis V : Freestream velocity u: Surge velocity (rigid state) α: Angle of attack (rigid state) ∗ PhD

Candidate, Kevin T. Crofton Department of Aerospace Engineering, Student Member, AIAA Candidate, Kevin T. Crofton Department of Aerospace Engineering, Student Member, AIAA ‡ Norris and Wendy Mitchell Endowed Professor, Kevin T. Crofton Department of Aerospace Engineering, Associate Fellow, AIAA § Fellow, AIAA † PhD

1 of 27 American of Aeronautics and Astronautics Copyright © 2018 by the American Institute of Aeronautics and Astronautics, Inc. AllInstitute rights reserved.


θ: Pitch attitude (rigid state) q: Pitch rate (rigid state) δj : Deflection of jth control surface Xu , Zu , Mu : Rigid forces and moments due to surge velocity (u) Xα , Zα , Mα : Rigid forces and moments due to angle of attack (α) Xq , Zq , Mq : Rigid forces and moments due to pitch rate (q) Xδ , Zδ , Mδ : Rigid forces and moments due to deflection of control surface (δ) ηi : Modal displacement coordinate for the ith mode (elastic state) η˙ i : Modal displacement rate coordinate for the ith mode (elastic state) Migen : Modal generalized mass of ith mode ωi : Natural frequency of the ith mode ζi : Damping ratio of the ith mode Zηi : Elastic lift due to ith modal displacement Mηi : Elastic moment due to ith modal displacement Zη˙ i : Elastic lift due to ith modal displacement rate Mη˙ i : Elastic moment due to ith modal displacement rate Ξαi : Generalized force due to ith mode for the rigid body degrees of freedom (α) Ξqi : Generalized force due to ith mode for the rigid body degrees of freedom (q) Ξiδj : Generalized force due to ith mode for the rigid body degrees of freedom (δj ) Ξiηj : ith mode generalized force due to jth mode (modal displacement) Ξiη˙ j : ith mode generalized force due to jth mode (modal displacement rate) A: Plant matrix ARR : Rigid part of the plant matrix ARE : Rigid to elastic aerodynamically coupled part of the plant matrix AER : Elastic to rigid aerodynamically coupled part of the plant matrix AE : Elastic part of the plant matrix B: Input Matrix BR : Rigid part of the control matrix BE : Elastic part of the control matrix C: Output Matrix D: Feedthrough Matrix

I.

Introduction

This paper focuses on the development of a computational framework for the analysis of aircraft longitudinal flight dynamics based on the mean-axis formulation.1 Aircraft with high aspect ratio wings are aerodynamically more efficient, but there are trade-offs between increased aspect ratio (AR), stiffness, and weight. As the AR is increased and weight is decreased, the increased flexibility can induce performance issues such as high aeroelastic deflections and a decrease in the flutter speed. Therefore, it is imperative to develop a flight dynamics model which includes an interaction between rigid and elastic degrees of freedom. Normally a fully coupled six degree of freedom model is required for rigid body flight dynamics. To include the effects of elastic deformation, the physics based full order mathematical model uses a additional n Degrees of freedom (nDoF) to model the elastic effects. Previous work by Waszak and Schmidt1 studied and formulated the effects of wing flexibility on aircraft flight dynamics. They showed that the frequency separation between the rigid body modes and elastic modes decreased as the flexibility of the aircraft structure is increased. Buttrill and Zeiler2 developed a nonlinear flight mechanics model with linear nDOF aeroelastic dynamics using Lagrangian mechanics. It was noted that the flexiblity not only affects the dynamics but also changes the total aircraft mass distribution which alters the inertia tensor. Zerweckh3 developed a flexible aeroelastic model of an aircraft with a high aspect ratio wing. The analysis was compared with the flight test data which showed that the phugoid mode was slightly unstable because of the flexibility. Similarly, work done by Patil4 on investigating nonlinear aeroelasticity and flight dynamics showed that elastic modes with low frequencies affects both short-period and phugoid modes. Chang5 et al. performed flight dynamics analysis of a high aspect ratio flexible aircraft to study the effects of different design parameters like wing flexibility, and horizontal/vertical tail aerodynamics. Recently Schmidt6, 7 analyzed the longitudinal dynamics of the Body Freedom Flutter (BFF) vehicle8 and showed that the elastic modes with low natural frequencies affects both the short period and phugoid modes. It was also concluded that the short period 2 of 27 American Institute of Aeronautics and Astronautics


mode coupled with the 1st symmetric bending leading to BFF at high speeds. Aircraft flexibility can be leveraged to change the shape of the wing during flight. This is achieved by deflecting the control surfaces during flight to achieve the required wing shape. Therefore, it becomes important to consider the controllability and observability of the modes to ensure all the modes to be controlled and observable during flight. Suzuki9 et al. minimized wing structural weight by simultaneously designing the structure and control system. Control stability and maximum stresses were used as design constraints. Idan10 developed an aeroservoelastic design framework to study the interaction between the structure and control system and optimize closed loop control margins and flutter. Nam11 developed an aeroservoelastic framework for flutter suppression and gust load alleviation for an aircraft using two control surfaces with structural, aerodynamic, and control design variables to minimize the control performance index and root mean square (RMS) value of the gust response. Perez12 et al. used Multidisciplinary Design Optimization (MDO) methods to integrate flight dynamics and control surface sizing in the aircraft conceptual design process. Results showed better performance and flying qualities than the traditional approach of designing the control laws after the structure is fixed. Stanford13 demonstrated a technique to determine the optimal control surface layout. For a transport aircraft wing he used structural wing box sizing and actuation variables to study both the static and dynamic aeroelastic phenomenon. However, the controllability and observability aspects while performing the aeroserovoelastic based design optimization have not been studied previously, but will be addressed here. This study will focus on the development of a general flight dynamics modeling and analysis framework with integrated controllability and observability metrics that can be integrated into an MDAO framework. This modeling framework was then used to assess the dynamic characteristics of the mAEWing2 flying-wing aircraft14, 16 . The flight dynamics analysis requires information about the aircraft’s structure and aerodynamics. A FEM model was developed using MSC.NASTRAN 18 to determine the free vibration modes of the aircraft. The aerodynamics were determined using the Vortex Lattice Method code, Tornado.19 The process was verified by analyzing the precursor of the mAEWing2, the mAEWing120, 21 aircraft. The flight dynamics of mAEWing1 has previously been studied and results published,22 so it serves as a verification case of the developed modeling framework and software. Once the verification was performed, flight dynamics analysis was performed for the mAEWing2 aircraft following the approach presented in previous work.25

II.

Motivation

The ultimate goal of this work is to integrate the flight dynamics modeling and simulation into a aircraft multidisciplinary design analysis and optimization environment (MDAO). The MDAO framework has been previously developed and implemented to optimize the mAEWing2 aircraft.14 The developed MDAO framework consists of three major modules: modeling, analysis, and optimization. The modeling module consists of software to create the geometry and mesh of the aircraft as well as creation of the analysis models for structural analysis and aeroelasticity. The analysis module consists of flutter analysis, static aeroelastic trim analysis, and induced and profile drag computation. The analysis module also contains methods for post processing structural results such as calculation of material failure indices, buckling loads, etc. Lastly the optimization model contains the main optimizer algorithm, response computation, and design variable output to the modeling module. The optimization algorithm is a hybrid strategy using both gradient-based optimization (GBO) and particle-swarm optimization (PSO) algorithms.14 To generate the mAEWing2 aircraft model a base structural and outer-mold-line (OML) geometry was created. The aircraft OML was created using a 50% scaled geometry of the X56 aircraft.17 An optimization of the base aircraft model was performed with the objective of minimizing weight and induced drag while having constraints on material stresses, buckling, the stall lift limit, hinge moments, and flutter velocity. The structure was subjected to a 3.8g load at 0.85Vbf f . Design variables used in the optimization were spar and rib locations, the battery location, laminate fiber ply orientations, ply thicknesses, size of flaps, and flap rotations. The optimal structural weight was 35.54 lbs. This optimized structure was then used in a second optimization process to determine the flap deflection schedule that produced the minimum total drag on the aircraft.14 The MDAO framework previously developed has been expanded and improved. The new Bilevel Programming optimization process was used to optimize the mAEWing2 aircraft using a novel internal structural

3 of 27 American Institute of Aeronautics and Astronautics

layout called SpaRibs.23 The new BLP optimization showed an improved -17.7% reduction in weight from the baseline mAEWing2 aircraft.15 This MDAO framework will be further expanded through inclusion of the flight-dynamics modeling methodology presented here.

III.

Methodology


The Flight Dynamics Model is developed using the process described in the flowchart in Figure 1.24 Two types of information are needed to construct the model. Firstly, the structural free vibration mode shapes, frequencies, and generalized masses of the aircraft are required to capture elastic motion. These were determined using a FEM developed using MSC.NASTRAN consisting of shell and beam elements. Secondly, information regarding the aerodynamics of the aircraft is needed. The rigid body lift distribution and the longitudinal stability derivatives are determined using a vortex-lattice code.19 The vehicle geometry was determined from a CAD model of the aircraft. Once the geometric, structural, and aerodynamics information is gathered, the linear rigid body flight dynamics model can be set up using the state variable format. The free vibration mode shapes, frequencies, and generalized masses are used to add the elastic degrees of freedom to the rigid model to assemble the fully elastic model. The rigid lift distribution along with the free vibration mode shapes and generalized masses are used to determine the aeroelastic stability derivatives. The aeroelastic stability derivatives are used to capture the elastic effect on the aerodynamic forces, moments, and the generalized forces acting on the aircraft, and are the essence of the aeroelastic model. More details about calculating these stability derivatives can be found in previous publications.24 The generalized modal forces are determined using the principle of virtual work. The fully elastic model is then used to develop a static elastically adjusted model by residualizing the effects of the elastic-modal displacement rates and considering only the effects of static elastic deformation on the aircraft. The assembly of the the static-elastically adjusted model and the other two models is explained further in Section III.B.

Figure 1: Flight dynamics modeling (FDM) framework


A.

Aeroelastic Stability Derivatives

The aeroelastic model is constructed using aeroelastic stability derivatives, analogous to the dimensional stability derivatives, such as Mα , used in rigid-aircraft flight dynamics. The rigid-body, and aeroelastic stability derivatives model the forces due to rigid-body and elastic motion, and are expressed in terms of non-dimensional stability derivatives such as the contribution of the pitching moment due to change in angle of attack, CM α . The non-dimensional derivatives may be estimated numerically several ways, including strip theory, computational aerodynamic techniques, or flight and wind-tunnel tests. For this study, strip theory was utilized for estimating the aeroelastic derivatives, following the approach given in MFD24 Chapter 7, using the wing section lift distribution, clα (y) obtained from a VLM code. The rigid-body stability derivatives were estimated using the same vortex-lattice code. This modeling approach allows for stability derivatives for a given vehicle to be obtained from different sources, and to be updated if better numerical estimates become available. The longitudinal rigid and aeroelastic non-dimensional stability derivatives are defined in Table 1.


Table 1: Non-dimensional rigid and aeroelastic longitudinal stability derivatives Rigid

Aeroelastic

CM α =

∂CM ∂α

CXα =

∂CX ∂α

∂CZ c ∂q 2V

CZq =

CM δ =

∂CM ∂δ

CXδ =

∂CX ∂δ

∂CZ ∂α

CM ηi =

∂CM ∂ηi

CZηi =

∂CZ ∂ηi

CM q =

∂CM c ∂Cq 2V

CM η˙ i =

∂CM ∂ η˙ i

CZ η˙ i =

∂CZ ∂ η˙ i

CXq =

∂CX c ∂q 2V

CQαi =

∂CQi ∂α

CQqi =

∂CZ ∂δ

CQiδj =

∂CQi ∂δj

CQiηj =

CQiηj ˙ =

∂CQi ∂ η˙ j

CZα =

CZδ =

∂CQi ∂q ∂CQi ∂ηj

The dimensional stability derivatives defined in terms of the non-dimensional derivatives are shown in Table 2. Table 2: Dimensional rigid and aeroelastic longitudinal stability derivatives Rigid Mα =

1 ρV 2 Sc CM α 2Iyy

Aeroelastic 1 ρV 2 S CZα 2m

Mηi =

1 ρV 2 Sc CM ηi 2Iyy

Zηi =

1 ρV 2 S CZηi 2m

1 ρV 2 Sc CM q 2Iyy

Mηi ˙ =

1 ρV 2 Sc CM η˙ i 2Iyy

Zηi ˙ =

1 ρV 2 S CZ η˙ i 2m

Zα =

Xα =

1 ρV 2 S CXα 2m

Zq =

1 ρV 2 S CZq 2m

Xq =

1 ρV 2 S CXq 2m

Ξαi =

1 ρV 2 Sc CQαi 2Migen

Ξqi =

1 ρV 2 Sc CQqi 2Migen

1 ρV 2 Sc CM δ 2Iyy

Zδ =

1 ρV 2 S CZδ 2m

Ξiδj =

1 ρV 2 Sc CQiδj 2Migen

Ξiηj =

1 ρV 2 Sc CQiηj 2Migen

Ξiηj ˙ =

1 ρV 2 Sc CQiηj ˙ 2Migen

Mδ =

Xδ =

1 ρV 2 S CXδ 2m

Mq =


The generalized forces (Ξi. ) are determined using the principle of virtual work, where the virtual work is determined due to the forces and moments acting through a virtual elastic displacement, and can be expressed as, Z δW

b/2

=

(δL0 δw + δM0 δθ)dy

(1)

−b/2

δM0 δW

= δMac + eδL0 n n X X = −δL0 cos αv φi ηi + (δmac + eδL0 cos αv φ0i ηi ) i=1

(2) (3)

i=1


Here, δL0 and δM0 are the forces acting on the 2D section of the wing and αv is the total angle of attack of the vehicle. Using this expression the elastic effects on the generalized forces for the rigid body degrees of freedom, the control surfaces, modal displacement (ηj ), and the modal displacement rate (η˙ j ) are determined as: Ξqj

=

∂(δW ) ∂qj

(4)

More details on the derivation of these stability derivatives can be found in Chapter 7 of MFD24 . B.

Dynamics Models

The dynamic modeling of the aircraft consists of the fully elastic model and the static elastically adjusted model in addition to a rigid model. More details about the derivation can be found in the Chapter 8 of MFD.24

1.

Fully Elastic Model (nDOF Model)

The three equations governing perturbations in rigid body translational motion, are written as: fx m

u˙ =

(V0 r + R0 v) − (Q0 w + W0 q) − g cos Θ0 θ +

v˙

=

(P0 w + W0 p) − (R0 u + U0 r) + g(cos Θ0 cos Φ0 φ − sin Θ0 sin Φ0 θ) +

w˙

=

(Q0 u + U0 q) − (P0 v + V0 p) − g(cos Θ0 sin Φ0 φ + sin Θ0 cos Φ0 θ) +

(5) fy m

fz m

(6) (7)

where fx , fy , and fz are the linear perturbation aerodynamic forces (excluding propulsive forces) acting in the x, y, and z direction of the body-fixed axis. The equations of motion governing the perturbations of the rigid body rotational motion are written as: p˙ −

r˙ −

Ixz r˙ Ixx

=

q˙

=

Ixz p˙ Izz

=

1 Ixz (Q0 p + P0 q) + (Iyy − Izz )(R0 q + Q0 r) + l Ixx 1 ( Izz − Ixx )(R0 p + P0 r) + 2Ixz (R0 r − P0 p) + m Iyy 1 − Ixz (R0 q + Q0 r) + (Ixx − Iyy )(Q0 p + P0 q) + n Izz

(8) (9) (10)

where l, m, and n are the linear perturbation aerodynamic moments (excluding propulsive forces) acting about the x, y, and z axes. Including the aeroelastic model, the aerodynamics forces acting on the translational degrees of freedom are written in terms of the stability derivatives and the rigid-body motion and


elastic deformation, or fx m

= Xu u + Xα α + Xα˙ α˙ + Xq q + XδE δE +

fy m

= Yβ β + Yp p + Yr r + YδA δA + YδR δR +

fz m

= Zu u + Zα α + Zα˙ α˙ + Zq q + ZδE δE +

n X

(Xηi ηi + Xη˙i η˙i )

i=1 n X

(11)

(Yηi ηi + Yη˙i η˙i )

(12)

(Zηi ηi + Zη˙i η˙i )

(13)

i=1 n X i=1

Similarly, the aerodynamic moments acting on the aircraft rotational degrees of freedom can be written as: l


Ixx

=

Lβ β + Lr r + LδA δA + LδR δR +

n X

(lηi ηi + lη˙i η˙i )

(14)

i=1

m Iyy

=

n Izz

=

Mu u + Mα α + Mα˙ α˙ + Mq q + MδE δE + Nβ β + Np p + Nr r + NδA δA + NδR δR +

n X

(mηi ηi + mη˙i η˙i )

(15)

(nηi ηi + nη˙i η˙i )

(16)

i=1 n X i=1

The equations of motion governing the perturbations in the elastic degrees of freedom are written as, ηï + 2ζi ωi η˙i + ωi2 ηi =

Ξi Migen

(17)

These equations 5 - 17 may now be expressed in state-variable format, or x˙

= Ax + Bu

(18)

y

= Cx + Du

(19)

And noting that writing the state vector x in terms of rigid-body then elastic degrees of freedom, or states, the above model maybe partitioned as, " # AR ARE A= (20) AER AE " # BR B= (21) BE " # CR C= (22) CE h i D= D (23) Here (AR ) contains the rigid body dynamics, (ARE ) contains the rigid-to-elastic aerodynamic coupling, (AER ) contains elastic-to-rigid aerodynamic coupling, and (AE ) contains purely aeroelastic dynamics. For longitudinal dynamics, for example, including three elastic degrees of freedom, the state vector maybe written as: h i xT = u α θ q η1 η2 η3 η˙1 η˙2 η˙3 For these ten states the corresponding plant matrix (A) becomes:




Xu


   Zu  V     0  Mu A=  0    0   0   0    0 0

Xα

−g

Zα V

0

0 Mα 0 0 0 Ξ1α Ξ2α Ξ3α

0 0 0 0 0 0 0 0

Xq 1+

Zq V

1 Mq 0 0 0 Ξ1q Ξ2q Ξ3q

0

0

0

0

0

0

Z η1 V

Z η2 V

Zη3 V

Z η˙1 V

Z η˙2 V

Z η˙3 V

0 Mη1 0 0 0 Ξ1η1 − ω12 Ξ2η1 Ξ3η1

0 Mη2 0 0 0 Ξ1η2 Ξ2η2 − ω22 Ξ3η2

0 Mη3 0 0 0 Ξ1η3 Ξ2η3 Ξ3η3 − ω32

0 Mη˙1 1 0 0 Ξ1η˙1 − 2ζω1 Ξ2η˙1 Ξ3η˙1

0 Mη˙2 0 1 0 Ξ1η˙2 Ξ2η˙2 − 2ζω2 Ξ3η˙2

0 Mη˙3 0 0 1 Ξ1η˙3 Ξ2η˙3 Ξ3η˙3 − 2ζω3 (24)

Or, 

AR

   A=  0  AER

 ARη ARη˙ | {z }  ARE  0 I    Aη Aη˙  | {z }

(25)

AE

For the vehicle being considered, there are five control surfaces along each wing. Consider only symmetric deflections of these surface yields five control inputs to the longitudinal dynamics. So the B matrix is written as: 

Xδ1

   Zδ  1  V     0   Mδ 1 B=  0    0   0  Ξ  1δ1  Ξ2δ1 Ξ3δ1 2.

Xδ2

Xδ3

Xδ4

Zδ2 V

Zδ3 V

Zδ4 V

0 Mδ2 0 0 0 Ξ1δ2 Ξ2δ2 Ξ3δ2

0 Mδ3 0 0 0 Ξ1δ3 Ξ2δ3 Ξ3δ3

0 Mδ4 0 0 0 Ξ1δ4 Ξ2δ4 Ξ3δ4

Xδ5



  Zδ5   V     0   Mδ5   0    0   0   Ξ1δ5    Ξ2δ5  Ξ3δ5

Static-Elastically Adjusted Model

This model was developed from Chapter 7, Section 7.11 of MFD,24 and includes the effects of the deformed wing shape (i.e. the effects of static elastic deformation) on the rigid body degrees of freedom. Recall that, the flight dynamics equation as shown in Eq. 18-23 in state variable format can be written as, " # " #" # " # h i x˙R AR ARE xR BR = + (26) u x˙E AER AE xE BE " #" # h i h ih i CR xR + D u (27) y = C E xE Note that


                       

h x˙ TE = η˙ 1

η˙ 2

η˙ 3

η¨1

η¨2

i h η¨3 = η˙ i

ηï

h xTE = η1

η2

η3

η˙ 1

η˙ 2

i h η˙ 3 = ηi

η˙ i

i i

The AE (purely elastic) portion of the A matrix (in Eq. 24) is arranged such that it separates out the forces due to modal displacement (Aη ) and modal displacement rate (Aη˙ ) as: " # " # h i 0 0 I h i h ih i x˙E = (28) xR + xE + BE u AER Aη Aη˙ The static-elastic model is found by setting η˙i , and ηï to zero which yields x˙ E = 0 and xE = [η0 , 0], or 0


η0

= AER xR + Aη η0 + BE u =

−A−1 η (AER xR

+ BE u)

(29) (30)

which is solved for the static modal displacement η0 . The new aerodynamic forces and moments due to the static elastic effects can be determined by substituting the value of static modal deformation (η0 ) for the xE vector in the dynamics of the rigid body degrees of freedom (x˙ R ) in Eq. 26 as: x˙R = AR xR + ARE xE + BR u

(31)

x˙R = AR xR + [ARη ARη˙ ][η0 0]T + BR u

(32)

Also, recall from Eq. 25 So, the static-elastically adjusted state space model is written as: x˙R

= AR xR + ARη η0 + BR u

x˙R

=

(33)

AR xR + ARη (−A−1 η (AER xR + −1 (AR − ARη Aη AER )xR + (BR

x˙R

=

x˙R

= ASEA xR + BSEA u

BE u)) + BR u −

ARη A−1 η BE )u

(34) (35) (36)

Where, ASEA and BSEA are the static elastically adjusted plant and input matrices. To calculate the static-elastically adjusted C and D matrices, consider Eq. 27 written as, h ih i h ih i h ih i y = C R xR + C E xE + D u (37) Similar to Eq. 32, the output y can be rewritten as: y = CR xR + [Cη Cη˙ ][η0 0]T + Du

(38)

Substituting the value of η0 from Eq. 30 into Eq. 38 gives, y

=

CR xR + Cη (−A−1 η (AER xR + BE u)+ Du Cη A−1 η AER )xR

y

=

(CR −

y

=

CSEA xR + DSEA u

+ (D −

Cη A−1 η BE )u

(39) (40) (41)

Where, CSEA and DSEA are the static elastically adjusted output and feedthrough matrices. C.

Impulse Residues: Controllability and Observability Metrics

The impulse residues were developed from Chapter 10, Section 10.1.2 of MFD.24 The residue matrices (R) are collectively a measure of modal controllability and observability. The residues are important as they determine the magnitude of the contribution of each of the modes to the systems time response. The relationship between the residues and controllability and observability is as follows:, as if there is a pole-zero cancellation for a particular mode, in the transfer function for a given input and output, that mode is not controllable or not observable and the residue will be zero.


The matrix of the transer functions, TF(s), of the system defined by Eqs. 18 and 19, assuming D=0, can be represented as, y(s) T F (s)

[C(sI − A)−1 B]u(s) 1 = [CM [diag ]M −1 B] s − λi

=

(42) (43)

Here, M is the modal matrix of A with columns consisting of eigenvectors (νi ) and M −1 is the inverse of modal matrix with rows consisting of (µi ), or left eigenvectors. Expanding Eq. 43 for the (k, j) transfer function yields, T F (s)k,j


T F (s)k,j

= =

n X (ck νi )(µi bj ) i=1 n X i=1

(s − λi ) Rk,i,j (s − λi )

(44) (45)

which can be recognized as the partial-fraction expansion of T F (s)k,j Note that the impulse residue (Rk,i,j ) is, Rk,i,j

=

(ck νi )(µi bj ) = Ok,i Ci,j

(46)

That is, the impulse residue (R) is the product of controllability (Ci,j ) and observability (Ok,i ) metrics. The Ci,j is a measure of controllability of mode i from input j i.e. how the j th control input affects the ith modal response, and Ok,i is a measure of observability of mode i from response k i.e. how the ith modal response affects the k th physical response.

IV.

Longitudinal Flight Dynamics Analysis

Figure 2: Longitudinal flight dynamics analysis (LFDA) process The longitudinal flight dynamics analysis (LFDA) is now performed using the longitudinal flight dynamics modeling framework (see Figure 1). Figure 2 shows the LFDA flowchart. In the LFDA process, a response for the given control input, system poles, flutter analysis and controllability and observability metrics are calculated. The controllability and observability metrics are calculated using impulse residues as shown in the next section.


V. A.

Verification: mAEWing1

Aircraft Description & Model Parameters


In order to verify the FDM framework and flight dynamics analysis process, an analysis of the mAEWing1 aircraft was conducted and the results compared to earlier results. The mAEWing120 aircraft is a flying wing 10 ft span aircraft based closely on the BFF vehicle designed by Lockheed Martin.8 The CAD model of the aircraft is shown in Figure 3a. The propulsion battery and the flight computer are mounted inside the centerbody. The aircraft has a pusher propeller with the motor attached to the aft portion of the centerbody. The aircraft has three trailing edge flaps and one body flap on each side. The flaps are numbered as shown in Figure 3b. Basic information about the aircraft such as weight, wing span, and aspect ratio are given in Table 3.

(a)

(b)

Figure 3: a) mAEWing1 CAD geometry and b) mAEWing1 flaps numbering

Table 3: mAEWing1 aircraft parameters Weight 14.7 lb

Span 10 ft

Aspect Ratio 8.66

Area 11.55 ft2

Taper Ratio 0.285

Leading Edge Sweep (Λ) 22◦

The wing is constructed of a primary load bearing solid spar surrounded by foam that the provides aerodynamic shape to the wing. The foam is then covered in a thin layer of fiberglass to provide a smooth aerodynamic surface. The fiberglass skin is cut in chordwise segments so that it does not carry load. The solid spar was manufactured using composite laminated reinforcements surrounding the core. Table 4 contains information about the free vibration modes computed with a FEM model using MSC.NASTRAN .18 Three elastic modes are considered for the longitudinal flight dynamics modeling, the first, second, and third symmetric vibration modes. Table 4: mAEWing1 free-free symmetric vibration Modal Summary Mode 1 2 3

Type of Mode 1 Symm. Vibration 2nd Symm. Vibration 3rd Symm. Vibration st

Natural Frequencies 6.07 Hz 12.15 Hz 19.47 Hz

Damping Ratio (assumed) 2% 2% 2%

Generalized Masses 1.00 slug-ft/in 1.00 slug-ft/in 1.00 slug-ft/in

Table C (see Appendix) shows the longitudinal aeroelastic stability derivatives for the mAEWing1 aircraft. These were used to assemble the fully elastic model and the static-elastically adjusted model. They


compare favorably to the previous results.22 B.

LFDA & Flutter Results

The longitudinal flight dynamics analysis was performed for a steady level flight condition at a speed of 85.3 ft/s and flight altitude of 1000 ft. The system poles from a rigid model, and full and static-elastically adjusted models were determined and are tabulated in Table 5. The terms in the square brackets denote the natural frequency and damping ratio of the mode. These poles match well with previous results with less than 5% difference. Results computed using the current FDM framework and flight dynamics analysis process are labeled ”VT”, and the previously published results from Schmidt & Associates for comparison are labeled ”S & A”.22


Table 5: mAEWing1: System poles for rigid, fully elastic, and static-elastically adjusted model Mode 1 (VT) 1 (S & A) 2 (VT) 2 (S & A) 3 (VT) 3 (S & A) 4 (VT) 4 (S & A) 5 (VT) 5 (S & A)

Type of mode Short Period Short Period Phugoid Phugoid st 1 Aeroelastic 1st Aeroelastic 2nd Aeroelastic 2nd Aeroelastic 3rd Aeroelastic 3rd Aeroelastic

Rigid model [10.04,0.79] [10.0,0.80] [0.38,-0.019] [0.38,-0.020] [-,-] [-,-] [-,-] [-,-] [-,-] [-,-]

Fully elastic model [13.94,0.68] [14.0,0.70] [0.43,-0.017] [0.43,-0.020] [33.99,0.065] [33.4,0.07] [71.24,0.027] [71.1,0.030] [122.5,0.04] [122.0,0.04]

Static-elastically adjusted model [13.92,0.73] [14.2,0.73] [0.43,-0.018] [0.44,-0.020] [-,-] [-,-] [-,-] [-,-] [-,-] [-,-]

Figure 4: mAEWing1 centerbody pitch rate step response using L3-R3 (symmetric) -1◦ deg deflection

The pitch rate step responses of the centerbody for a -1◦ deg symmetric deflection of L3 - R3 flaps is shown in Figure 4. The response from all three models are also in very good agreement with the previous results (S&A).22 These pitch rate responses will be discussed in the following sections.


Results - Rigid Model The pitch rate step response (rigid model) of the centerbody for a -1◦ deg symmetric deflection of flaps L3 R3 is shown in Figure 4. For this model, two rigid modes are observed (see Table 5). The higher frequency short period mode involves oscillations of the pitch attitude, angle of attack, and a small contribution from the surge velocity. It has a natural frequency of 10.04 rads/sec with the damping ratio of 0.79. The eigenvectors of the low frequency phugoid mode reveal very low contribution from the angle of attack and pitch rate. The aircraft exhibits oscillation at almost constant angle of attack with large contribution from the surge velocity and moderate contribution from the pitch attitude. This slightly unstable phugoid mode has a natural frequency of 0.38 rads/sec and damping ratio of -0.019 which, as shown in Table 5.


Results - Fully Elastic Model The fully elastic model was developed using the aeroelastic stability derivatives tabulated in Table C in Appendix. The reulting system poles for the fully elastic model are tabulated in Table 5. Two rigid and three symmetric aeroelastic modes are observed. For the short period mode, the natural frequency has increased by almost 40% to 13.94 rads/sec with approximately a 10% decrease in the damping ratio to 0.7 compared to the rigid model. An inspection of the eigenvectors reveals the additional contribution of the three aeroelastic degrees of freedom (elastic pitch rate associated with three modes) namely θ˙E1 , θ˙E2 and θ˙E3 , to this mode. The second highest contribution came from the θ˙E1 degree of freedom, unlike the short period mode observed for the rigid model. So, this short-period-like mode can be characterized as an elastically coupled rigid-body mode or an elastic short period mode. The second rigid mode is the phugoid mode. It is a slightly unstable mode like the rigid model and has a natural frequency of 0.43 rads/sec with a damping ratio of -0.017. The remaining three modes are the first three symmetric aeroelastic modes. The first aeroelastic mode is slightly damped with a natural frequency of 34.0 rads/sec. The eigenvectors show that the highest contribution comes from the θ˙E1 degree of freedom and the second highest contribution comes from the rigid body pitch rate followed by the θ˙E2 degree of freedom. It also has small contributions from the surge velocity, angle of attack, and the pitch attitude. So, this mode can be characterized as a rigidelastically coupled mode. At higher speeds, the elastic short period mode couples with the first aeroelastic mode which leads to Body Freedom Flutter (BFF). The second aeroelastic mode is a slightly damped mode with a natural frequency of 71.2 rads/sec. The eigenvectors of this mode show that the highest contribution comes from the θ˙E2 degree of freedom and the second highest contribution comes from the θ˙E1 degree of freedom followed by the rigid body pitch rate. The coupling between the first and the second aeroelastic modes leads to Bending-Torsion Flutter at high speeds. The third aeroelastic mode is a slightly damped mode with a natural frequency of 122.4 rads/sec. This mode gets most of its contribution from the θ˙E3 degree of freedom and very little contribution from the second mode. Hence little modal coupling is observed regarding this mode. Results - Static-Elastically Adjusted Model This model was assembled by residualizing the effects of modal velocity and modal accelerations on the fully elastic model. This allows for the determination of the effects of static elastic deformation on the rigid body dynamics (i.e. the effects of deformed wing shape). The resulting static-elastically adjusted stability derivatives are tabulated in Table 6. In comparison to the rigid stability derivatives, there is a significant increase in the values of the CLα , CM α , CLq and CM q due to the wing static deformation. A similar trend is observed for the lift and pitching moment stability derivatives of the flaps. The pitch rate response of the centerbody using symmetric deflection of L3 - R3 flaps at -1◦ was shown in Figure 4. In comparison to the rigid model, the response for the static-elastically adjusted model has lower magnitude due to the reduced static-elastic pitch control power. The system poles for this model were also tabulated in Table 5. Two rigid modes with the effects of static elastic deformations are observed. The short period mode has a higher frequency of 13.9 rads/sec than that from the rigid model. This is a result of the flexibility of the wing as the aeroelastic increase the magnitude of CM α . This derivative acts like a torsional spring constant in the pitching motion of the vehicle. The phugoid mode is slightly unstable with a frequency of 0.43 rads/sec. The response of the static elastically adjusted model agrees well with the rigid-body or mean-axes response from the full order model. 13 of 27 American Institute of Aeronautics and Astronautics

Table 6: mAEWing1: Comparison of rigid and static elastically adjusted (SEA) longitudinal stability derivatives (/rad)

Derivatives

CL (Rigid)

CL (SEA)

CM (Rigid)

CM (SEA)

C.α *C.q C.δ1 C.δ2 C.δ3 C.δ4

4.59 4.13 0.77 0.61 0.50 0.43

5.91 5.76 0.58 0.50 0.59 0.86

-0.16 -1.85 0.021 -0.048 -0.20 -0.31

-0.39 -2.07 0.030 -0.061 -0.24 -0.35


* These coefficients have units of seconds (sec).

Results - Flutter Analysis Flutter analysis was also performed for the mAEWing1 aircraft. Figure 5 shows the velocity root locus plot, with the modes identified. The phugoid mode remains near the origin with very low damping. But with increasing velocity the short period moves further away from the imaginary axis while the branch associated with first aeroelastic mode crosses the imaginary axis, leading to Body Freedom Flutter (BFF). This occurs due to high interactions between the elastic short period mode and the first aeroelastic mode. The BFF occurs at a flight speed of 63 kts and a frequency of 4.8 Hz. The branch associated with the second aeroelastic mode leads to a classical bending - torsion flutter occurring at a flight speed of 90 kts at 9.75 Hz. The branch associated with the third aeroelastic mode moves away from the imaginary axis as the flight velocity increases.

Figure 5: mAEWing1 velocity root locus plot (30-90) kts Table 7 shows the comparison between the flutter results from VT NASTRAN, Schmidt & Associates (S & A) model, and the current longitudinal flight dynamics model. The results for the flutter analysis agree very well with the VT NASTRAN and (S & A). All these results show proper implementation of the modeling approach.


Table 7: mAEWing1: flutter analysis comparison Model

Flutter Speed (kts)

Flutter Frequency (Hz.)

VT NASTRAN (S & A) VT LFDA

63.6 63.0 63.0

5.10 4.80 4.85

VI.


A.

Flight Dynamics Analysis: mAEWing2

Aircraft Geometry and Configuration

Attention now turns to the modeling of the second vehicle, mAEWing2. The mAEWing214, 16 is a low speed swept back flying wing with a 14 ft. span and 37.5 lbs of weight. This vehicle was developed, in part, to support research on modeling and control of flexible aircraft. The vehicle has a total of twelve flaps for longitudinal and lateral/directional stability. The entire wing composes of four trailing edge flaps, one body flap, and one leading edge flap on each side. A CAD diagram of the aircraft planform geometry with flaps labeling can be seen in Figure 6. Basic properties for the aircraft are also tabulated in Table 8.

Figure 6: mAEWing2 flaps numbering

Table 8: mAEWing2 aircraft parameters Weight 37.48 lbs

Span 14.00 ft

Aspect Ratio 10.66

Area 18.38 ft2

Taper Ratio 0.26

Leading Edge Sweep (Λ) 22.00◦

The aircraft is electrically powered with a motor mounted on the rear of the centerbody. The avionics, flight computers, and batteries are contained in the centerbody as shown in Figure 7a. Unlike the mAEWing120 aircraft, mAEWing2 has a landing gear also which is modeled as a lumped mass. The mAEWing2 aircraft was designed at Virginia Tech using Multidisciplinary Design Analysis and Optimization (MDAO) tools to design the internal structural layout of the wing and the centerbody. The internal structure of the wing (Figure 7a) consists of straight spars and ribs, made of sandwich panels. The top and bottom skin of the wings and centerbody are constructed using symmetric carbon fiber fabric laminates. Both wings and winglets were designed to be removable to allow for interchangeability of the wings and centerbody.


The mAEWing214, 16 aircraft was fabricated at the University of Minnesota (UMN). This aircraft is a scaled version of the X56-A designed by Lockheed Matrin Skunk Works.17 Figure 7b below shows the internal structural layout of the mAEWing216 aircraft at UMN.


(a)

(b)

Figure 7: a) AEWing2 aircraft CAD model and b)mAEWing2 aircraft manufactured at UMN

B.

Structural Model

The internal structural layout of this aircraft was designed for a 3.8g manuever load factor using the MDAO14 framework. The FEM Model of mAEWing2 is shown in Figure 8a.14 The full details of the mAEWing2 modeling can be found in previous work.14 The free-free vibration analysis was performed using MSC. NASTRAN.18 Figure 8b shows the first six vibration modes from the FEM model.

(a)

(b)

Figure 8: a) mAEwing2 internal structural layout and b) mAEwing2 elastic mode shapes

Figure 9: mAEWing2 equivalent beam FEM model



To perform the dynamic analysis, an equivalent beam model was developed to determine the elastic axis and the mass normalized mode shapes along this elastic axis. Figure 9 shows the unique location of the elastic axis along the wing. Unlike the usual elastic axis that lies within the structure, the elastic axis of this aircraft follows a curved path and does not remain inside the wing near the outer span. This is an artifact of both designing the wing using composite materials and achieving a required flutter speed within the flight envelope. The mode shapes from the detailed FEM model were interpolated along the elastic axis using the finite element shape functions. Figure 10 shows the first six elastic modes along the elastic axis.

Figure 10: mAEwing2 elastic mode shapes along the elastic axis Mode 7, 9, and 10 are the first three symmetric vibration modes and Mode 8, 11, and 13 are the first three anti-symmetric vibration modes. For the longitudinal flight dynamics analysis, only symmetric modes are considered as asymmetric modes do not contribute towards the longitudinal stability of the aircraft. Modes 7, 9 and 10 were considered for the flight dynamics analysis. The natural frequencies are summarized in Table 9. Table 9: mAEWing2 free-free symmetric vibration modes Mode 1 2 3

C.

Type of Mode 1 Symm. Vibration 2nd Symm. Vibration 3rd Symm. Vibration st

Natural Frequencies (Hz.) 5.54 17.52 19.46

Aerodynamics Model

A model of the rigid-body aerodynamics was developed for the mAEWing2 aircraft using the Vortex Lattice Method (VLM) code.19 The planform of the aircraft and airfoil shapes used to generate the rigid-body aerodynammic model were determined from CAD geometry of the aircraft. The rigid body stability derivatives were also calculated using the VLM code and are shown in Table B (see Appendix). Figure 11a shows the aerodynamic mesh for the aircraft which consists of 552 panels, and the rigid 2-D section lift distribution (clα ) can be seen in Figure11b.


(a)

(b)


Figure 11: a) mAEWing2 aerodynamic mesh and b) mAEWing2 rigid lift distribution

D.

FDA & Flutter Results

A longitudinal flight dynamics analysis was performed for the rigid model, the fully elastic model, and the static-elastically adjusted model. To investigate all three models, an eigenanalysis was performed to calculate the system poles for each of the models. Next, the centerbody pitch rate response for the symmetric L4 R4 flap deflections was studied. To analyze the magnitude of each of the modes to the centerbody pitch rate response, the controllability and observability metrics were calculated for the four trailing edge flaps. Finally, a flutter analysis was performed to estimate the critical flutter speed. The results match very well with the previous results of (S&A).25 Results - Rigid Body Model The dynamic analysis for the rigid aircraft was performed for a steady level flight condition at a speed of 93.2 ft/s and flight altitude of 1000 ft. The centerbody pitch rate step response was determined for a -1◦ symmetric deflection of the L4-R4 flaps and is shown in Figure 12. This response represents a rapid, well damped pitch rate with little overshoot.

Figure 12: mAEWing2 rigid centerbody pitch rate step response The transfer function for the pitch attitude response is shown below, θcg (s) −65.67[7.28][0.09] = deg/deg −δ4 (s) [0.60, 11.49][0.085, 0.40] 18 of 27 American Institute of Aeronautics and Astronautics

(47)

The terms in the square brackets provide the location of the poles and zeros. Two terms in brackets list the damping ratio and natural frequency of two complex poles or zeros. The response is quite conventional for a rigid aircraft, reflecting both classical phugoid and short period modes. Results - Fully Elastic Model


To explore the effects of flexibility on the flight dynamics of mAEWing2, the fully elastic model was analyzed at a flight velocity of 92.3ft/s (below the flutter speed). This model was assembled using the aeroelastic stability derivatives model discussed in Section III.B.1. The aeroelastic stability derivatives are tabulated in Table D in the Appendix. An eigenanalysis was performed to determine the system poles for the fully elastic model, and they are tabulated in Table 10 with identified short period mode, phugoid mode, first AE mode, second AE mode, and the third AE modes. In comparison to the rigid model, both the short-period and phugoid mode from the fully elastic model have significantly less change in the frequency and damping. This is in contrast to what was observed for the mAEWing1 aircraft in Table 5 which had a large difference in these modes. Table 10: mAEWing2: System poles comparison from full-order and rigid model Mode

Type of Mode

Rigid Model [ω,ζ]

Full Order Model [ω,ζ]

1 2 3 4 5

Short Period Phugoid 1st AE 2nd AE 3rd AE

[11.49,0.60] [0.38,0.085] [-,-] [-,-] [-,-]

[12.34,0.69] [0.37,0.09] [28.81,0.079] [110.14,0.007] [121.78,0.045]

Figure 13: Polar representation of modes Figure 13 shows the polar diagram of eigenvectors for the phugoid and short period modes. A conventional phugoid mode was observed as the surge velocity dominates the modal response. From Table 10, this mode has a natural frequency of 0.37 rads/sec with a damping ratio of 0.09. The second mode for the fully elastic model in Table 10 is the short period mode. This has most of the contribution from rigid pitch rate and the second largest contribution from the elastic pitch rate (θ˙E1 ) associated with the first elastic degree of freedom. So, unlike a conventional short period mode, which is dominated by the rigid pitch rate response, this mode exhibits strong coupling with the first elastic degree of freedom and can be referred as elastic short period mode. The natural frequency for this mode is 12.34 rads/sec with damping ratio of 0.69. The third mode (in Table 10) is the first AE mode, which has a natural frequency of 28.8 rads/sec and a damping ratio of 0.08. This mode is dominated by the elastic pitch rate θ˙E1 degree of freedom. It is also coupled with the rigid body pitch rate, which has the second highest contribution. At high speeds, the elastic short period mode interacts with the first AE mode and leads to BFF. Finally, the third highest contribution comes from the θ˙E2 degree of freedom.


The fourth mode in Table 10 is the second AE mode which is a high frequency mode with a natural frequency of 110.14 rads/sec and a very low damping ratio of 0.007. The highest contribution comes from θ˙E2 degree of freedom, and the second largest contribution arises from the θ˙E1 degree of freedom, and the third highest from the rigid body pitch rate. Finally, the last mode in Table 10 is the third AE mode which is a slightly damped high frequency mode with a natural frequency of 121.78 rads/sec. The θ˙E3 degree of freedom dominates the response. The second largest contribution comes from the θ˙E1 degree of freedom followed by the θ˙E2 degree of freedom. A centerbody pitch rate step response for a -1◦ symmetric deflection of the L4 - R4 flaps for the rigid and fully elastic model are shown in Figure 14. The response for the fully elastic model was determined using Eq. 48, i=3 X (48) θEi (y = cg)ηi (t) qf lex (t) = qrig (t) + | {z } i=1 1 {z } | Downloaded by VIRGINIA TECH on November 7, 2018 | http://arc.aiaa.org | DOI: 10.2514/6.2018-3425

2

The first part of the above equation represents the contribution of the rigid pitch rate (qrig (t)) and the second part represents the effect of the local elastic deformation (θEi (y = cg)ηi (t)) at the cg location of the aircraft. The (qf lex (t)) represents the total flexible pitch rate response about the cg location of the mAEWing2 aircraft.

Figure 14: mAEWing2 flexible body pitch rate step response

Figure 15: Elastic response (θ˙E ) of centerbody from the three elastic degrees of freedom From Figure 14 the elastic contribution to this response is evident. To explore this elastic contribution 20 of 27 American Institute of Aeronautics and Astronautics

further, consider Figure 15. This figure shows the contribution of the three elastic degrees of freedom (θEi (y = cg)ηi (t)) to the pitch rate response in Figure 14. It can be seen that the first elastic degree of freedom contributes almost all of the elastic response. The contribution from the third elastic degree of freedom damps out eventually, however there is very little but constant contribution coming from the second elastic degree of freedom (pink color). Results - Impulse Residue: Controllability and Observability Metrics


To explore the extent of contribution of each of the modes in detail, the impulse residue metrics were calculated as discussed in Section IV. For the mAEWing2, the metrics were determined for the centerbody pitch rate response from symmetric deflections of the four trailing edge flaps (L2/R2-L5/R5) at a flight velocity of 92.3 ft/s. Table 11 shows the residue (controllability and observability) metrics for the four different inputs. Figure 16 shows the centerbody pitch rate step responses to a -1◦ symmetric flap deflections. The L2-R2 flap combination has the lowest response, and the L5-R5 flap combination has the highest pitch rate response. Table 11: mAEWing2 residue metrics for pitch rate response (rad/sec/rad

Modes

L2-R2

L3-R3

L4-R4

L5-R5

Phugoid Mode Short Period 1st AE Mode 2nd AE Mode 3rd AE Mode

0.16 9.05 22.15 0.702 3.43

0.58 30.51 14.31 3.98 10.61

0.88 47.89 49.64 0.41 10.73

1.03 56.88 85.47 1.85 1.34

Figure 16: mAEWing2 centerbody pitch rate step response using TE flaps From Table 11 and Figure 16, the following observations can be made: • For the L2-R2 flaps, the first AE mode with a residue of 22.15 rad/sec/rad dominates the pitch rate response, and is consistent with the step response shown in Figure 16. The short period modes has the second highest residue at 9.05 rad/sec/rad, and the third AE mode has the third highest residue at 3.43 rad/sec/rad. As expected (Table 11), the phugoid mode has almost negligible effect on the pitch rate response.


• For L3-R3 flap deflection response in Figure 16, the short period mode dominates the pitch rate response with a residue of 30.51 rad/sec/rad. From Table 11, the first and third AE modes have the second and third highest residues at 14.31 rad/sec/rad, and 10.61 rad/sec/rad. The second AE mode has the least contribution amongst all three AE modes. All these responses are consistent with the step response for L3-R3 in Figure 16 • Using the L4-R4 flaps excites both the short period mode and the first AE mode with almost equal residue magnitude at 47.89 rad/sec/rad and 49.64 rad/sec/rad. The third AE mode has the third highest contribution at 10.73 rad/sec/rad followed by very little contribution from the second AE mode. These results are consistent with the step response in Figure 16.


• The combination of the L5-R5 flaps excites the first AE mode the most at 85.47 rad/sec/rad followed by the short period mode at 56.88 rad/sec/rad, and very little and almost equal contributions from the second at 1.85 rad/sec/rad, and the third AE mode at 1.34 rad/sec/rad. The phugoid mode has very little excitation. From Figure 16, a non-minimum phase behavior is also evident as a -1◦ deflection should induce a positive pitch rate however, a negative pitch rate is seen earlier in the response. So, from the results in Table 11, it can be concluded that the first AE mode dominates the response for all the four symmetric flap deflection except for the L3-R3 combination, which are located near the node line of the first symmetric vibration mode. The short period mode has the second highest contribution except for the L3-R3 flaps, where the short period modes dominates. The third AE mode has the third highest contribution from all four flaps except from the L5-R5 combination. The second AE mode and phugoid mode have very little contribution for symmetric deflections of all the flaps. Results - Static-Elastically Adjusted Model The static elastically adjusted model was studied to investigate the effect of static deformation on the rigid body dynamics of mAEWing2. This model was assembled by residualizing the effects of modal velocity and modal acceleration as was discussed in Section III.B.2. The resulting rigid body stability derivatives were updated with static elastic corrections and are tabulated in Table 12. In comparison to the rigid body stability derivatives, there is a significant increase in the values of the CLα and CLq due to the static deformations of the aircraft. For the outermost trailing edge flaps (L4-R4 and L5-R5), the CLδ stability derivatives increased by more than two times. This large increase in lift stability derivatives highlights the effects of static elastic deformations. However, the CM α , CM q , and CM δi (for i = 1 to 5) derivatives for the aircraft and flaps are very similar to those of the rigid model. This is in contrast to what was observed for mAEWing1 aircraft in Table 6. Table 12: mAEWing2: Comparison of rigid and static elastically adjusted (SEA) longitudinal stability derivatives (/rad)

Derivatives

CL (Rigid)

CL (SEA)

CM (Rigid)

CM (SEA)

C.α *C.q C.δ1 C.δ2 C.δ3 C.δ4 C.δ5

4.61 4.19 0.64 0.48 0.41 0.37 0.29

7.19 6.52 0.377 0.203 0.345 0.681 0.973

-0.27 -2.02 0.0011 -0.042 -0.15 -0.23 -0.25

-0.25 -2.06 -0.0046 -0.045 -0.136 -0.207 -0.242


This difference is reflected in the centerbody pitch rate response of the static-elastically adjusted model for mAEWing2 (see Figure 17) in comparison to the mAEWing1 (see Figure 4). For mAEWing2 aircraft, the magnitude of the static elastic pitching moment coefficients are almost the same as those from the rigid



Figure 17: mAEWing2 centerbody pitch rate step response

model. Hence the static-elastic response is closer to that from the rigid model. The pitch attitude transfer function for the static-elastically adjusted model is shown on the left, and for comparison the transfer function for the rigid model is shown on the right in Table 13. Table 13: Transfer functions for static-elastically adjusted and rigid model

Static-elastically adjusted

Rigid

θcg (s) −64.5[10.78][0.092] = deg/deg −δ4 (s) [0.76, 11.90][0.09, 0.36]

θcg (s) −65.67[7.28][0.09] = deg/deg −δ4 (s) [0.60, 11.49][0.085, 0.40]

In comparison to the rigid model, the magnitude of the larger zero of the transfer function for the staticelastically adjusted has increased as it is directly proportional to the increase in value of CLα (see Table 12). The short period mode remains almost unchanged with a frequency of 11.90 rads/sec. The phugoid mode also remains the same with a natural frequency of 0.36 rads/sec with damping of 0.09.

Results - Flutter Analysis Flutter analysis was performed using the full-order model. Figure 18 shows the velocity root locus with the branches associated with the system modes identified. The phugoid mode remains close to the origin. With increasing velocity the short period moves further away from the imaginary axis. At higher flight speeds, this short period mode couples with the first symmetric aeroelastic mode leading to Body Freedom Flutter (BFF), which occurs at a flight speed of 74.5 kts with a frequency of 4.01 Hz. A similar trend was observed for the BFF7 aircraft and for the mAEWing1 aircraft, as discussed in the previous section.



Figure 18: mAEWing2 velocity root locus plot (30-80 kts)

Table 14 shows a comparison of the flutter results for the VT FEM model analyzed using NASTRAN and the current fully elastic model. The results for the flutter analysis computed using NASTRAN are in close agreement with the full-order model. Table 14: mAEWing2: flutter analysis comparison Model

Flutter Speed (kts)

Flutter Frequency (Hz.)

VT FEM Model (NASTRAN) VT Fully elastic model

78.5 74.5

3.62 4.01

VII.

Conclusion

A flight dynamics modeling and analysis methodology was developed for use in a MDAO framework for aircraft aeroservoelastic optimization. The flight dynamics methodology is based on the mean-axis formulation developed to include the elastic effects. Using this methodology three different dynamic models were created: 1) Rigid 2) Fully elastic 3) Static-elastically adjusted. The fully elastic model was developed using the aeroelastic stability derivatives. This model was residualized to develop the static-elastically adjusted model. The methodology was verified by applying it to the mAEWing1 aircraft and the results agreed well with the previously published results. Finally, a longitudinal flight dynamics analysis was performed for a second UAV aircraft mAEWing2. The rigid model of the mAEWing2 aircraft revealed conventional short period and phugoid modes. For the fully elastic model, an eigenanalysis revealed significant differences in the system modes compared to the rigid model. It was found that the short-period and the first AE mode exhibited coupling leading to BFF at a speed of 74.5 kts, with a flutter frequency of 4.01 Hz. Below the flutter speed, the centerbody elastic pitch rate response showed significant aeroelastic contributions to this response. The first elastic degree of freedom dominated the elastic pitch rate response. This effect of flexibility was also evident in the static-elastically adjusted model, which showed large differences in the static-elastically adjusted stability derivatives compared to the rigid stability derivatives. Finally, the controlability and observability metrics were analyzed


to determine the relative importance of the modes to the input-output dynamics of the aircraft. It was observed that for the centerbody pitch rate response, the first AE mode dominated the response using all but one flap pair. Hence, it is concluded that the fully elastic dynamics of the mAEWing2 differ significantly from the rigid dynamics and including them in the MDAO framework will be useful to capture the effects of flexibility early in the conceptual design stage of the aircraft.

Acknowledgments This research was funded by NASA as a part of the ”Performance Adaptive Aeroelastic Wing (PAAW)” project, NRA NNX14AL36A with Dr. Jeff Ouellette as the Technical Monitor. The authors would also like to thank Dr. Wei Zhao and Dr. Mohamed Jrad, both of Virginia Tech, for their help on this work.

References Downloaded by VIRGINIA TECH on November 7, 2018 | http://arc.aiaa.org | DOI: 10.2514/6.2018-3425

1 Waszak,

M. R., and Schmidt, D. K., ” Flight Dynamics of Aeroelastic Vehicles,” Journal of Aircraft, Vol. 25, No. 6, June 1988, pp. 563-571. 2 Buttrill, C. S., Zeiler, T. A., and Arbuckle, P. D., ”Nonlinear Simulation of a Flexible Aircraft in Maneuvering Flight,” AIAA Paper 87-2501-CP. AIAA Flight Simulation Technologies Conference (Monterey, CA), August 1987. 3 Zerweckh, S. H., vonFlotow, A. H., and Murray, J. E., ”Flight Testing of a Highly flexible Aircraft Case Study On The MIT Light Eagle,” Journal of Aircraft, Vol. 27, No. 4, 1990, pp. 342-349. 4 Patil, M. J., Hodges, D. H., and Cesnik, C. E. S., ”Nonlinear Aeroelasticity and Flight Dynamics of High-Altitude Long-Endurance Aircraft,” Proceedings of the 40th Structures, Structural Dynamics and Materials Conference, Saint Louis, MO, AIAA Paper 99-1470, Apr. 1999, pp. 2224-2232. 5 Chang, C. S., Hodges, D. H., and Patil, M. J., ”Flight Dynamics of Highly Flexible Aircraft,” Journal of Aircraft, Vol. 45, No. 2 (2008), pp. 538-545. 6 Schmidt, D. K., ”Stability Augmentation and Active Flutter Suppression of a Flexible Flying-Wing Drone,” AIAA Guidance, Navigation, and Control Conference, AIAA SciTech Forum, (AIAA 2016-2099) 7 Schmidt, D. K., ”MATLAB-Based Flight-Dynamics and Flutter Modeling of a Flexible Flying-Wing Research Drone,” Journal of Aircraft, Vol. 53, No. 4 (2016), pp. 1045-1055. 8 Burnett, E. L., Atkinson, C., Beranek, J., Sibbitt, B., Holm-Hansen, B., and Nicolai, L., ”NDOF Simulation Model for Flight Control Development with Flight Test Correlation,” Lockheed Martin Aeronautics Co., AIAA Paper No. 2010-7780, Toronto, August 2010. 9 Suzuki, S., and Yonezawa, S., ”Simultaneous structure/control design optimization of a wing structure with a gust load alleviation system,” Journal of Aircraft, Vol. 30, No. 2 (1993), pp. 268-274 10 Idan, M., Karpel, M., and Moulin, B., ”Aeroservoelastic Interaction Between Aircraft Structural and Control Design Schemes,” Journal of Guidance, Control, and Dynamics, Vol. 22, No. 4(1999), pp. 513-519 11 Nam, C., Chattopadhyay, A., and Kim, Y., ”Optimal Wing Planform Design for Aeroelastic Control,” AIAA Journal, Vol. 38, No. 8 (2000), pp. 1465-1470. 12 Perez, R. E., Hugh H. T. L., and Behdinan, K., ”Multidisciplinary Optimization Framework for Control-Configuration Integration in Aircraft Conceptual Design,”Journal of Aircraft, Vol. 43, No. 6 (2006), pp. 1937-1948. 13 Stanford, B., ”Optimal Control Surface Layout for an Aeroservoelastic Wingbox,” 58th AIAA/ASCE/AHS/ASC Structures, Structural Dynamics, and Materials Conference, AIAA SciTech Forum, (AIAA 2017-1813). 14 Zhao, W., Jrad, M., Gupta, R., and Kapania, R. K., ”Multidisciplinary Design Analysis and Optimization of Performance Adaptive Aeroelastic Wings,” AIAA SciTech Conference, Grapevine, TX, Jan. 9-13, 2017. 15 Zhao, W. and Kapania, R. K., ”BLP Optimization of Composite Flying-wings with SpaRibs and Multiple Control Surfaces”, 2018 AIAA/ASCE/AHS/ASC Structures, Structural Dynamics, and Materials Conference, AIAA SciTech Forum, (AIAA 2018-2150) 16 Regan,

C. D., ”mAEWing2: Conceptual Design and System Test,”AIAA Atmospheric Flight Mechanics Conference,


AIAA SciTech Forum, (AIAA 2017-1391). 17 http://www.lockheedmartin.com/us/aeronautics/skunkworks.html 18 MSC/NASTRAN

Quick Reference Guide: version 68. MacNeal-Schwendler Company, 2013.

19 Melin, T., A Vortex Lattice MATLAB Implementation for Linear Aerodynamic Wing Applications, Dept. of Aeronautics, Royal Institute of Technology, Sweden, Dec., 2000. 20 Regan, C. D., and Taylor, B., ”mAEWing1: Design, Build, Test,” Proceedings of the AIAA SciTech Conference, San Diego, CA, Jan. 4-8, 2016. 21 Zhao, W., Muthirevula, N., Kapania, R. K., Gupta, A., Regan, C. D., and Seiler, P. J., ”Finite Element Model Updating of A Small Flexible Composite UAV,” AIAA Atmospheric Flight Mechanics Conference, AIAA SciTech Forum, (AIAA 2017-1393).


22 Schmidt, D. K., Zhao, W., and Kapania, R.K., ”Flight-Dynamics and Flutter Modeling and Analyses of a Flexible Flying-Wing Drone - Invited,” AIAA Atmospheric Flight Mechanics Conference, AIAA SciTech Forum, (AIAA 2016-1748). 23 Locatelli, D., Mulani,B. S., and Kapania, R.K., ”Wing-Box Weight Optimization Using Curvilinear Spars and Ribs (SpaRibs)” Journal of Aircraft, Vol. 48, No. 5 (2011), pp. 1671-1684. 24 Schmidt,

D. K., Modern Flight Dynamics, McGraw-Hill, 2012.

25 Schmidt,

D. K., Dynamic Analysis mAEWing2 FEM 2.1V1.2, DKS PAAW Working Paper, 2018.

Appendix - Rigid and Aeroelastic Stability Derivatives for mAEWing1 and mAEWing2 Table B: mAEWing2 rigid body longitudinal stability derivatives (/rad)

Derivatives

CL Derivatives

CD Derivatives

CM Derivatives

C.α *C.q C.δ1 C.δ2 C.δ3 C.δ4 C.δ5

4.61 4.19 0.64 0.48 0.41 0.37 0.29

0.078 0.04 0.017 0.009 0.0001 0.002 0.003

-0.271 -2.02 0.0011 -0.042 -0.15 -0.23 -0.25




Table C: mAEWing1 aeroelastic stability derivatives (/rads)

CLη1 CLη2 CLη3 *CLη˙1 *CLη˙2 *CLη˙3 CM η 1 CM η 2 CM η 3 *CM η˙1 *CM η˙2 *CM η˙3

-1.543 -1.422 0.254 0.218 0.368 -0.676 0.015 1.157 -0.205 -0.474 -0.09 0.401

CQ1η1 CQ2η1 CQ3η1 CQ1η2 CQ2η2 CQ3η2 CQ1η3 CQ2η3 CQ3η3 *CQ1η˙1 *CQ2η˙1 *CQ3η˙1

0.142 0.157 -0.194Z 0.525 0.497 -0.529 -0.013 0.069 -0.061 -0.265 -0.105 0.116

*CQ1q *CQ2q *CQ3q CQ1α CQ2α CQ3α CQ1δ1 CQ2δ1 CQ3δ1 CQ1δ2 CQ2δ2 CQ3δ2

-0.415 -0.4277 0.315 -0.447 -0.465 0.586 0.100 0.003 -0.103 0.084 -0.060 0.002

*CQ1η˙2 *CQ2η˙2 *CQ3η˙2 *CQ1η˙3 *CQ2η˙3 *CQ3η˙3 CQ1δ3 CQ2δ3 CQ3δ3 CQ1δ4 CQ2δ4 CQ3δ4

-0.025 -0.096 0.103 0.147 0.232 -0.306 -0.008 -0.085 0.215 -0.185 -0.013 0.149

* These coefficients have units seconds.

Table D: mAEWing2 Aeroelastic Stability Derivatives (/rads)

CLη1 CLη2 CLη3 *CLη˙1 *CLη˙2 *CLη˙3 CM η 1 CM η 2 CM η 3 *CM η˙1 *CM η˙2 *CM η˙3 CQ1δ5

-1.43 0.15 0.32 -0.13 -0.15 -0.53 -0.034 0.12 -0.14 -0.22 0.12 0.16 -0.10

*CQ1q *CQ2q *CQ3q CQ1α CQ2α CQ3α CQ1δ1 CQ2δ1 CQ3δ1 CQ1δ2 CQ2δ2 CQ3δ2 CQ3δ5

-0.25 0.038 0.29 -0.35 0.11 0.47 0.031 -0.008 -0.005 0.036 -0.003 0.040 0.015

CQ1η1 CQ2η1 CQ3η1 CQ1η2 CQ2η2 CQ3η2 CQ1η3 CQ2η3 CQ3η3 *CQ1η˙1 *CQ2η˙1 *CQ3η˙1 CQ2δ5

0.10 0.009 -0.071 -0.032 -0.00035 0.009 -0.14 -0.007 0.024 -0.063 0.028 -0.019 0.012

* These coefficients have units of seconds


*CQ1η˙2 *CQ2η˙2 *CQ3η˙2 *CQ1η˙3 *CQ2η˙3 *CQ3η˙3 CQ1δ3 CQ2δ3 CQ3δ3 CQ1δ4 CQ2δ4 CQ3δ4 –

-0.0013 -0.009 -0.0348 -0.027 -0.047 -0.231 0.011 0.017 0.125 -0.044 0.0282 0.126 –