Coupled Vortex-Lattice Flight Dynamic Model with Aeroelastic Finite-Element Model of Flexible Wing Transport Aircraft with Variable Camber Continuous Trailing Edge Flap for Drag Reduction

Coupled Vortex-Lattice Flight Dynamic Model with Aeroelastic

Finite-Element Model of Flexible Wing Transport Aircraft with

Variable Camber Continuous Trailing Edge Flap for Drag

Reduction

∗ Nhan Nguyen NASA Ames Research Center, Moffett Field, CA 94035 † Eric Ting Stinger Ghaffarian Technologies, Inc., Moffett Field, CA 94035 ‡ Daniel Nguyen University of California, Berkeley , CA 94720 § Tung Dao Stinger Ghaffarian Technologies, Inc., Moffett Field, CA 94035 ¶ Khanh Trinh Stinger Ghaffarian Technologies Inc., Moffett Field, CA 94035 This paper presents a coupled vortex-lattice flight dynamic model with an aeroelastic finite-element model to predict dynamic characteristics of a flexible wing transport aircraft. The aircraft model is based on NASA Generic Transport Model (GTM) with representative mass and stiffness properties to achieve a wing tip de- flection about twice that of a conventional transport aircraft (10% versus 5%). This flexible wing transport aircraft is referred to as an Elastically Shaped Aircraft Concept (ESAC) which is equipped with a Variable Camber Continuous Trailing Edge Flap (VCCTEF) system for active wing shaping control for drag reduction.

A vortex-lattice aerodynamic model of the ESAC is developed and is coupled with an aeroelastic finite-element model via an automated geometry modeler. This coupled model is used to compute static and dynamic aeroe- lastic solutions. The deflection information from the finite-element model and the vortex-lattice model is used to compute unsteady contributions to the aerodynamic force and moment coefficients. A coupled aeroelastic- longitudinal flight dynamic model is developed by coupling the finite-element model with the rigid-body flight dynamic model of the GTM.

I. Introduction The aircraft industry has been responding to the need for energy-efficient aircraft by redesigning airframes to be aerodynamically efficient, employing light-weight materials for aircraft structures, and incorporating more energy- efficient aircraft engines. Reducing airframe operational empty weight (OEW) using advanced composite materials is one of the major considerations for improving energy efficiency. Modern light-weight materials can provide less structural rigidity while maintaining sufficient load-carrying capacity. As structural flexibility increases, aeroelastic interactions with aerodynamic forces and moments can alter aircraft aerodynamics significantly, thereby potentially degrading aerodynamic efficiency.

Under the Fundamental Aeronautics Program at the NASA Aeronautics Research Mission Directorate, the Fixed Wing project is conducting multidisciplinary foundational research to investigate advanced concepts and technologies ∗ Research Scientist, AIAA Associate Fellow, Intelligent Systems Division, nhan.t.nguyen@nasa.gov † Engineer, Intelligent Systems Division, eric.b.ting@nasa.gov ‡ Student, AIAA Student Member, Mechanical Engineering Department, daniel.y.nguyen@berkeley.edu § Engineer, AIAA Student Member, Intelligent Systems Division, tung.x.dao@nasa.gov ¶ Engineer, Intelligent Systems Division, khanh.v.trinh@nasa.gov 1 of 38 American Institute of Aeronautics and Astronautics for future aircraft systems. A NASA study entitled “Elastically Shape Future Air Vehicle Concept” was conducted in 2010 to examine new concepts that can enable active control of wing aeroelasticity to achieve drag reduction. This study showed that highly flexible wing aerodynamic surfaces can be elastically shaped in-flight by active control of wing twist and vertical deflection in order to optimize the local angle of attack of wing sections to improve aerodynamic efficiency through drag reduction during cruise and enhanced lift performance during take-off and landing.

The study shows that active aeroelastic wing shaping control can have a potential drag reduction benefit. Conven- tional flap and slat devices inherently generate drag as they increase lift. The study shows that conventional flap and slat systems are not aerodynamically efficient for use in active aeroelastic wing shaping control for drag reduction. A new flap concept, referred to as Variable Camber Continuous Trailing Edge Flap (VCCTEF) system, was conceived 1, 2 by NASA to address this need. Initial results indicate that the VCCTEF system may offer a potential pay-off for drag reduction that will result in significant fuel savings. In order to realize the potential benefit of drag reduction by active aeroelastic wing shaping control, configuration changes in high-lift devices have to be a part of the wing shaping control strategy.

NASA and Boeing are currently conducting a joint study to develop the VCCTEF further under the research 3, 4 element Active Aeroelastic Shape Control (AASC) within the Fixed Wing project. This study built upon the devel- opment of the VCCTEF system for NASA Generic Transport Model (GTM) which is essentially based on the B757 airframe, employing light-weight shaped memory alloy (SMA) technology for actuation and three separate chord- wise segments shaped to provide a variable camber to the flap. This cambered flap has potential for drag reduction as compared to a conventional straight, plain flap. The flap is also made up of individual 2-foot spanwise sections which enable different flap setting at each flap spanwise position. This results in the ability to control the wing twist shape as a function of span, resulting in a change to the wing twist to establish the best lift-to-drag ratio (L/D) at any aircraft gross weight or mission segment. Current wing twist on commercial transports is permanently set for one cruise con- figuration, usually for a 50% loading or mid-point on the gross weight schedule. The VCCTEF offers different wing twist settings for each gross weight condition and also different settings for climb, cruise and descent, a major factor in obtaining best L/D conditions.

The second feature of VCCTEF is a continuous trailing edge flap. The individual 2-foot spanwise flap sections are connected with a flexible covering, so no breaks can occur in the flap platforms, thus reducing drag by eliminating these breaks in the flap continuity which otherwise would generate vorticity that results in a drag increase and also contributes to airframe noise. This continuous trailing edge flap design combined with the flap camber result in lower drag increase during flap deflections. In addition, it also offers a potential noise reduction benefit.

The VCCTEF serves multiple functions as: • A wing shaping control device to twist the flexible wing to obtain changes in lift-to-drag ratios that will reduce cruise drag throughout the flight envelope.

• A high-lift device for take-off, climb-out, let-down and final approach by using the full span cambered flap.

• A full span roll control effector in lieu of traditional ailerons using the aft section of the cambered flap.

• An aeroservoelastic (ASE) control device to compensate for reduced flutter margins of flexible wings.

This paper describes an aeroelastic formulation of a flexible wing aircraft based on a one-dimensional structural dynamic theory that models the wing structure as a beam in a coupled bending-torsion motion. The aeroelastic angle of attack is derived from kinematics of aircraft rigid-body velocities and wing structural deflection velocities. The resulting nonlinear aeroelastic equations of bending-torsion motion are coupled with the aircraft rigid-body flight dynamic equations of motion. The nonlinear aeroelastic formulation takes into account the engine thrust forces which are coupled with wing aeroelasticity as a force follower. The formulation therefore is an aero-propulsive-elasticity.

This inclusion of the propulsive effect may be important for aircraft with engine-mounted high flexible wing structures.

The aeroelastic analysis also takes into account wing pre-twist and dihedral angles which can cause a high degree of coupling between the wing aeroelastic deflections and the aircraft rigid-body motion.

In the present study, the aeroelasticity equations are transformed into a system of aeroelastic state-space equations using the finite-element method (FEM), a powerful numerical technique which converts a given system of partial dif- ferential equations into a truncated, discretized weak-form solution formulation which utilizes locally-defined basis functions (typically polynomials) to numerically approximate the solution to the governing partial differential equa- tions. In general, the standard finite-element method belongs to a class of numerical techniques known as weighted- residual methods, such as the Galerkin method, which seek to minimize the error between the “true” solution and the approximation space of basis functions. Mathematically, this property arises from the fact that the finite-element 2 of 38 American Institute of Aeronautics and Astronautics method solves a variational weak form of the initial boundary value problem using arbitrary test functions which be- long to the Hilbert space of functions which are square-integrable. Proper application of the finite-element method will result in the creation of a matrix system of equations which may be solved using standard numerical techniques.

A modal analysis is then performed to assess aeroelastic stability of the aircraft. Frequencies and damping ratios of the symmetric and anti-symmetric modes are computed.

A vortex-lattice model of the flexible wing GTM, herein referred to as Elastically Shaped Aircraft Concept (ESAC), is developed for coupling with the flight dynamic model and the aeroelastic finite-element model (FEM) to provide stability and control derivatives for the flight dynamic model and the aerodynamic loading information for the FEM. An automated geometry modeler developed in MATLAB provides an update of the aircraft wing deformed geometry for the vortex-lattice model.

II. Description of Elastically Shaped Aircraft Concept Fig. 1 - Generic Transport Model (GTM) and Remotely Piloted Vehicle at NASA Langley Fig. 2 - GTM Planform The elastically shaped aircraft concept is modeled as a notional single-aisle, mid-size, 200-passenger aircraft. The geometry of the ESAC is obtained by scaling up the geometry of NASA generic transport model (GTM) by a scale of 200:11. The GTM is a research platform that includes a wind tunnel model and a remotely piloted vehicle, as shown 3 of 38 American Institute of Aeronautics and Astronautics in Fig. 1. Figure 2 is an illustration of the GTM planform. The reason for selecting the GTM is that there already exists an extensive wind tunnel aerodynamic database that could be used for validation in the study. The benchmark configuration represents one of the most common types of transport aircraft in the commercial aviation sector that provides short-to-medium range passenger carrying capacities.

The aircraft has a take-off weight of 200,000 lbs for a typical operating load (gear up, flap up) that includes cargo, fuel, and passengers. Fuel weighs about 50,000 lbs for a range of about 3,000 nautical miles.

To compute the mass and inertia properties of the benchmark aircraft, a component-based approach is used. The aircraft is divided into the following components: fuselage, wings, horizontal tails, vertical tail, engines, operational empty weight (OEW) equipment, and typical load including passengers, cargo, and fuel. The fuselage, wings, horizon- tal tails, and vertical tail are modeled as shell structures with constant wall thicknesses. Based on publicly available data of component weight breakdown for various aircraft, an average wing mass relative to the total empty weight of the aircraft is taken to be 24.2% of the OEW.

To enable active wing shaping control, the wing structures of the ESAC are designed to increase wing flexibility.

The wing bending and torsional stiffness quantities are designed to achieve a wing deflection that is about double of that of a conventional aircraft wing. The VCCTEF is divided into 14 sections attached to the outer wing and 3 sections attached to the inner wing, as shown in Fig. 3. Each 24-inch section has three camber flap segments that can be individually commanded, as shown in Fig. 4. These camber flaps are joined to the next section by a flexible and supported material (shown in blue) installed with the same shape as the camber and thus providing continuous flaps throughout the wing span with no drag producing gaps. The flexible skin materials that cover the spanwise camber flap sections create constraints to the flap deflections. These constraints impose a certain relative flap deflection between any two adjacent spanwise flap sections.

Using the camber positioning, a full-span, low-drag, high-lift configuration can be activated that has no drag producing gaps and a low flap noise signature. This is shown in Fig. 5. To further augment lift, a slotted flap configuration is formed by an air passage between the wing and the inner flap that serves to improve airflow over the flap and keep the flow attached. This air passage appears only when the flaps are extended in the high lift configuration.

Because the wings are highly flexible, flutter margins can become a potential issue. Flight dynamics and control of a highly flexible wing aircraft must fully account for the effects of aeroelasticity.

Fig. 3 - GTM with VCCTEF Fig. 4 - Variable Camber Flap 4 of 38 American Institute of Aeronautics and Astronautics Fig. 5 - Cruise and High Lift VCCTEF Configurations III. Aeroelastic Analysis A. Reference Frames Fig. 6 - Aircraft Reference Frames Figure 6 illustrates three orthogonal views of a typical aircraft. Several reference frames are introduced to facilitate the rigid-body dynamic and structural dynamic analysis of the lifting surfaces. For example, the aircraft inertial reference frame A is defined by unit vectors a , a , and a fixed to the non-rotating earth. The aircraft body-fixed 1 2 3 reference frame B is defined by unit vectors b , b , and b aligned with the roll, pitch, and yaw axes, respectively. The 1 2 3 right wing elastic reference frame C is defined by unit vectors c , c , and c . The reference frames B and C are related 1 2 3 π by three successive rotations: 1) the first rotation about b by the sweep angle + Λ of the elastic axis that results in ′ ′ ′ ′ an intermediate reference frame B defined by unit vectors b , b , and b (not shown), 2) the second rotation about 1 2 3 ′′ ′ b by the dihedral angle Γ of the elastic axis that results in an intermediate reference frame C defined by unit vectors ′ ′ ′ ′ c , c , and c (not shown), and 3) the third rotation about c by an angle π that results in the reference frame C. This 1 2 3 1 5 of 38 American Institute of Aeronautics and Astronautics relationship is expressed as           b − sin Λ − cos Λ 0 cos Γ 0 sin Γ 1 0 0 c 1 1           =  b   cos Λ − sin Λ 0   0 1 0   0 − 1 0   c  2 2 b 0 0 1 − sin Γ 0 cos Γ 0 0 − 1 c 3 3     − sin Λ cos Γ cos Λ sin Λ sin Γ c     =  cos Λ cos Γ sin Λ − cos Λ sin Γ   c  (1) − sin Γ 0 − cos Γ c The left wing elastic reference frame D is defined by unit vectors d , d , and d . The reference frames B and D are 1 2 3 π related by three successive rotations: 1) the first rotation about − b by the elastic axis sweep angle + Λ that results ′′ ′′ ′′ ′′ in an intermediate reference frame B defined by unit vectors b , b , and b (not shown), 2) the second rotation about 1 2 3 ′′ ′ ′ b by the elastic axis dihedral angle Γ that results in an intermediate reference frame D defined by unit vectors d , 2 1 ′ ′ ′ d , and d (not shown), and 3) the third rotation about d by an angle π that results in the reference frame D. This 2 3 1 relationship is expressed as           b − sin Λ cos Λ 0 cos Γ 0 sin Γ 1 0 0 d 1 1             =         b − cos Λ − sin Λ 0 0 1 0 0 − 1 0 d 2 2 b 0 0 1 − sin Γ 0 cos Γ 0 0 − 1 d 3 3     − sin Λ cos Γ − cos Λ sin Λ sin Γ d     = (2)  − cos Λ cos Γ sin Λ cos Λ sin Γ   d  − sin Γ 0 − cos Γ d The aircraft velocity at the aircraft center of gravity (CG) in the aircraft body-fixed reference B is expressed in the reference frames C and D as     − sin Λ cos Γ cos Λ sin Λ sin Γ c [ ] 1     v = u v w  cos Λ cos Γ sin Λ − cos Λ sin Γ   c  − sin Γ 0 − cos Γ c = ( − u sin Λ cos Γ + v cos Λ cos Γ − w sin Γ ) c + ( u cos Λ + v sin Λ ) c 1 2 + ( u sin Λ sin Γ − v cos Λ sin Γ − w cos Γ ) c (3)     − sin Λ cos Γ − cos Λ sin Λ sin Γ d [ ] 1     v = u v w  − cos Λ cos Γ sin Λ cos Λ sin Γ   d  − sin Γ 0 − cos Γ d = ( − u sin Λ cos Γ − v cos Λ cos Γ − w sin Γ ) d + ( − u cos Λ + v sin Λ ) d 1 2 + ( u sin Λ sin Γ + v cos Λ sin Γ − w cos Γ ) d (4) where ( u , v , w ) are the aircraft velocity components in the forward, lateral, and downward directions defined by the unit vectors ( b , b , b ) , respectively.

1 2 3 Generally, the effect of the dihedral angle can be significant. In the analysis, the aeroelastic effects on the fuselage, horizontal stabilizers, and vertical stabilizer are not considered, but the analytical method can be formulated for ana- lyzing these lifting surfaces if necessary. In general, a whole aircraft analysis approach should be conducted to provide a comprehensive assessment of the effect of structural flexibility on aircraft performance and stability. However, the scope of this study pertains only to the wing structures.

B. Elastic Analysis In the subsequent analysis, the combined motion of the left wing is considered. The motion of the right wing is a mirror image of that of the left wing for symmetric flight. The wing has a varying pre-twist angle γ ( x ) commonly designed in many aircraft. Typically, the wing pre-twist angle varies from being nose-up at the wing root to nose-down at the wing tip. The nose-down pre-twist at the wing tip is designed to delay stall onsets. This is called a wash-out 6 of 38 American Institute of Aeronautics and Astronautics twist distribution. Under aerodynamic forces and moments, the aeroelastic deflections of a wing introduce stresses and strains into the wing structure. The internal structure of a wing typically comprises a complex arrangement of load carrying spars and wing boxes. Nonetheless, the elastic behavior of a wing can be captured by the use of equivalent stiffness properties. These properties can be derived from structural certification testing that yields information about wing deflections as a function of loading. For high aspect ratio wings, an equivalent beam approach can be used to analyze aeroelastic deflections with good accuracy. The equivalent beam approach is a typical formulation in many aeroelasticity studies. It is assumed that the effect of wing curvature is ignored and the straight beam theory is used to model the wing deflection. The axial or extensional deflection of a wing is generally very small and therefore can usually be neglected.

Consider an airfoil section on the left wing as shown in Fig. 7 undergoing bending and torsional deflections. Let ( x , y , z ) be the coordinates of point Q on a wing airfoil section. Then the undeformed local airfoil coordinates of point Q are [ ] [ ] [ ] y cos γ − sin γ η = (5) z sin γ cos γ ξ where η and ξ are local airfoil coordinates, and γ is the wing section pre-twist angle, positive nose-down.

Then differentiating with respect to x gives [ ] [ ] [ ] [ ] ′ ′ y − sin γ − cos γ η − z γ x = γ = (6) ′ z cos γ − sin γ ξ y γ x Fig. 7 - Left Wing Reference Frame of Wing in Combined Bending-Torsion Let Θ be a torsional twist angle about the x -axis, positive nose-down, and let W and V be flapwise and chordwise bending deflections of point Q, respectively. Then, the rotation angle due to the elastic deformation can be expressed as φ ( x , t ) = Θ d − W d + V d (7) 1 x 2 x 3 where the subscripts x and t denote the partial derivatives of Θ , W , and V .

Let ( x , y , z ) be the coordinates of point Q on the airfoil in the reference frame D and p be its position vector.

1 1 1 Then the coordinates ( x , y , z ) are computed using the small angle approximation as 1 1 1         x ( x , t ) x φ × ( y d + z d ) . d x − yV − zW 1 2 3 1 x x         = + = (8)  y ( x , t )   y + V   φ × ( y d + z d ) . d   y + V − z Θ  1 2 3 2 z ( x , t ) z + W φ × ( y d + z d ) . d z + W + y Θ 1 2 3 3 Differentiating x , y , and z with respect to x yields 1 1 1     ′ ′ x 1 − yV + z γ V − zW − y γ W 1 , x xx x xx x    ′ ′  = (9)  y   − z γ + V − z Θ − y γ Θ  1 , x x x ′ ′ z y γ + W + y Θ − z γ Θ 1 , x x x Neglecting the transverse shear effect, the longitudinal strain is computed as ds − ds s 1 1 , x ε = = − 1 (10) ds s x 7 of 38 American Institute of Aeronautics and Astronautics where √ √ ( ) ′ 2 2 2 2 2 s = 1 + y + z = 1 + ( y + z ) γ (11) x x x √ 2 2 2 s = x + y + z 1 , x 1 , x 1 , x 1 , x √ ( ) ( ) 2 2 ′ 2 ′ ′ 2 2 2 = s − 2 yV − 2 zW + 2 ( y + z ) γ Θ + ( x − 1 ) + y + z γ + z − y γ (12) xx xx x 1 , x 1 , x 1 , x x Ignoring the second-order terms and using the Taylor series expansion, s is approximated as 1 , x ( ) ′ 2 2 − yV − zW + y + z γ Θ xx xx x s ≈ s + (13) 1 , x x s x The longitudinal strain is then obtained as ( ) ′ 2 2 − yV − zW + y + z γ Θ xx xx x ε = s x [ ] [ ] [ ] ( ) ( ) ( ) 2 2 2 ( ) ( ) ( ) ( ) ′ ′ ′ ′ 2 2 2 2 2 2 2 2 ≈ − y 1 + y + z γ V − z 1 + y + z γ W + y + z γ 1 + y + z γ Θ (14) xx xx x ( ) ′ For a small wing twist angle γ , γ ≈ 0. Then ( ) ′ 2 2 ε = − yV − zW + y + z γ Θ (15) xx xx x The moments acting on the wing are then obtained as ( )       ( ) ′ 2 2 y + z γ + Θ M GJ Θ ∫∫ x x x       = + E ε dydz  M   0    y − z M 0 z − y   ( )   ′ ′ ′ Θ GJ + EB γ − EB γ − EB γ 1 2 3 x      ′  = (16)  W  xx  − EB γ EI − EI  2 yy yz ′ V xx − EB γ − EI EI 3 yz zz ′ where E is the Young’s modulus; G is the shear modulus; γ is the derivative of the wing pre-twist angle; I , I , and yy yz I are the section area moments of inertia about the flapwise axis; J is the torsional constant; and B , B , and B are zz 1 2 3 the bending-torsion coupling constants which are defined as     2 2 B ∫∫ y + z ( )     2 2 = y + z dydz (17)  B   z  B y The strain analysis shows that, for a pre-twisted wing, the bending deflections are coupled to the torsional deflection ′ via the slope of the wing pre-twist angle. This coupling can be significant if the wash-out slope γ is dominant for highly twisted wings.

C. Aeroelastic Angle of Attack The relative velocity of the air approaching a wing section includes the contribution from the wing elastic deflection that results in changes in the local angle of attack. Since aerodynamic forces and moments are dependent on the local angle of attack, the wing aeroelastic deflections will generate additional elastic forces and moments. The local angle of attack depends on the relative approaching air velocity as well as the rotation angle φ from Eq. (7). The relative air velocity in turn also depends on the deflection-induced velocity. The velocity at point Q due to the aircraft velocity and angular velocity in the reference frame D is computed as v = ¯ v + ω × r = ( u b + v b + w b ) + ( p b + q b + r b ) × ( − x b − y b − z b ) Q 1 2 3 1 2 3 a 1 a 2 a 3 = ( u + ry − qz ) b + ( v − rx + pz ) b + ( w + qx − py ) b = x d + y d + z d (18) a a 1 a a 2 a a 3 t 1 t 2 t 3 8 of 38 American Institute of Aeronautics and Astronautics where     x − ( u + ry − qz ) sin Λ cos Γ − ( v − rx + pz ) cos Λ cos Γ − ( w + qx − py ) sin Γ t a a a a a a     = (19)  y   − ( u + ry − qz ) cos Λ + ( v − rx + pz ) sin Λ  t a a a a z ( u + ry − qz ) sin Λ sin Γ + ( v − rx + pz ) cos Λ sin Γ − ( w + qx − py ) cos Γ t a a a a a a and ( p , q , r ) are aircraft angular velocity components in the roll, pitch, and yaw axes, and ( x , y , z ) is the coordinate a a a of point Q in the aircraft body-fixed reference frame B relative to the aircraft C.G. (center of gravity) such that x is a positive when point Q is aft of the aircraft CG, y is positive when point Q is toward the left wing from the aircraft a C.G., and z is positive when point Q is above the aircraft C.G.

a The local velocity at point Q due to aircraft rigid-body dynamics and aeroelastic deflections in the reference frame D is obtained as   x − ( z + W + y Θ ) W − ( y + V − z Θ ) V [ ] t xt xt   ˙ v = v + φ × p = v d + v d + v d = d d d  y + V − ( yV + zW ) V − ( z + W + y Θ ) Θ  (20) Q x 1 y 2 z 3 1 3 3 t t x x xt t z + W − ( yV + zW ) W + ( y + V − z Θ ) Θ t t x x xt t In order to compute the aeroelastic forces and moments, the velocity must be transformed from the reference frame D to the airfoil local coordinate reference frame defined by ( μ , η , ξ ) as shown in Fig. 2. Then the transformation can be performed using successive rotation matrix multiplication operations as           v 1 0 0 cos V sin V 0 cos W 0 sin W v μ x x x x x           =  v   0 cos ( Θ + γ ) sin ( Θ + γ )   − sin V cos V 0   0 1 0   v  η x x y v 0 − sin ( Θ + γ ) cos ( Θ + γ ) 0 0 1 − sin W 0 cos W v ξ x x z   cos V ( v cos W + v sin W ) + v sin V x x x z x y x   =  cos ( Θ + γ ) [ − sin V ( v cos W + v sin W ) + v cos V ] + sin ( Θ + γ ) ( − v sin W + v cos W )  x x x z x y x x x z x − sin ( Θ + γ ) [ − sin V ( v cos W + v sin W ) + v cos V ] + cos ( Θ + γ ) ( − v sin W + v cos W ) x x x z x y x x x z x   v + v V + v W x y x z x   ≈ (21)  − v [ V + W ( Θ + γ )] + v + v [( Θ + γ ) − V W ]  x x x y z x x v [ − W + V ( Θ + γ )] − v ( Θ + γ ) + v [ 1 + V W ( Θ + γ )] x x x y z x x for small deflections.

The local aeroelastic angle of attack on the airfoil section due to the velocity components v and v in the reference η ξ frame D, as shown in Fig. 8, is computed as ( ) v + w ¯ v + ∆ v + w v + w ¯ v + w ∆ v i i i i η ξ ξ ξ ξ ξ α = = = − (22) c v ¯ v + ∆ v ¯ v ¯ v η η η η η where w is the local downwash due to three-dimensional lift distribution over a finite wing and i ¯ v = ( u + ry − qz ) sin Λ sin Γ + ( v − rx + pz ) cos Λ sin Γ − ( w + qx − py ) cos Γ (23) a a a a a a ξ ∆ v = W − ( yV + zW ) W + ( y + V − z Θ ) Θ + v [ − W + V ( Θ + γ )] − v ( Θ + γ ) (24) ξ t x x xt t x x x y ¯ v = − u cos Λ (25) η ∆ v = − ( ry − qz ) cos Λ + ( v − rx + pz ) sin Λ + V − ( yV + zW ) V − ( z + W + y Θ ) Θ η a a a a t x x xt t − v [ V + W ( Θ + γ )] + v [( Θ + γ ) − V W ] (26) x x x z x x Assuming that the wing dihedral Γ is small and neglecting the local downwash w , the local aeroelastic angle of i 9 of 38 American Institute of Aeronautics and Astronautics attack is evaluated as ( u + ry − qz ) sin ΛΓ + ( v − rx + pz ) cos ΛΓ − ( w + qx − py ) a a a a a a α = − c u cos Λ + W − ( yV + zW ) W + ( y + V − z Θ ) Θ t x x xt t − − [ − W + V ( Θ + γ )] × x x u cos Λ [ − ( u + ry − qz ) sin Λ − ( v − rx + pz ) cos Λ − ( w + qx − py ) Γ − ( z + W + y Θ ) W − ( y + V − z Θ ) V ] a a a a a a xt xt × u cos Λ [ − ( u + ry − qz ) cos Λ + ( v − rx + pz ) sin Λ + V − ( yV + zW ) V − ( z + W + y Θ ) Θ ] ( Θ + γ ) a a a a t x x xt t + u cos Λ ( u + ry − qz ) sin ΛΓ + ( v − rx + pz ) cos ΛΓ − ( w + qx − py ) a a a a a a − × 2 2 u cos Λ { × − ( ry − qz ) cos Λ + ( v − rx + pz ) sin Λ + V − ( yV + zW ) V − ( z + W + y Θ ) Θ − [ V + W ( Θ + γ )] × a a a a t x x xt t x x × [ − ( u + ry − qz ) sin Λ − ( v − rx + pz ) cos Λ − ( w + qx − py ) Γ − ( z + W + y Θ ) W − ( y + V − z Θ ) V ] a a a a a a xt xt + [( Θ + γ ) − V W ] [( u + ry − qz ) sin ΛΓ + ( v − rx + pz ) cos ΛΓ − ( w + qx − py ) + W x x a a a a a a t } − ( yV + zW ) W + ( y + V − z Θ ) Θ ] (27) x x xt t Fig. 8 - Airfoil Local Coordinates The aeroelastic angle of attack is generally a nonlinear function of the structural deflections. A nonlinear aeroelas- tic theory based on this approach has been developed.

In this study, a linear aeroelasticity is developed by neglecting nonlinear deflection-dependent terms and letting u ≈ V , v ≈ V β , and w ≈ V α . Thus ∞ ∞ ∞ ( V + ry − qz ) sin ΛΓ + ( V β − rx + pz ) cos ΛΓ − ( V α + qx − py ) + W + y Θ ∞ a a ∞ a a ∞ a a t t α = − c V cos Λ ∞ [ − ( V + ry − qz ) sin Λ − ( V β − rx + pz ) cos Λ − ( V α + qx − py ) Γ ] W ∞ a a ∞ a a ∞ a a x + V cos Λ ∞ [ − ( V + ry − qz ) cos Λ + ( V β − rx + pz ) sin Λ ] ( Θ + γ ) ∞ a a ∞ a a + V cos Λ ∞ ( V + ry − qz ) sin ΛΓ + ( V β − rx + pz ) cos ΛΓ − ( V α + qx − py ) ∞ a a ∞ a a ∞ a a − × 2 2 V cos Λ ∞ { × − ( ry − qz ) cos Λ + ( V β − rx + pz ) sin Λ + V − z Θ a a ∞ a a t t + [( V + ry − qz ) sin Λ + ( V β − rx + pz ) cos Λ + ( V α + qx − py ) Γ ] V ∞ a a ∞ a a ∞ a a x } + [( V + ry − qz ) sin ΛΓ + ( V β − rx + pz ) cos ΛΓ − ( V α + qx − py )] ( Θ + γ ) (28) ∞ a a ∞ a a ∞ a a 10 of 38 American Institute of Aeronautics and Astronautics Assuming z ≈ 0, then evaluating the various partial derivatives of the local aeroelastic angle of attack yields α = − γ − tan ΛΓ ( 1 + tan ΛΓ γ ) (29) ∂ α 1 + 2 tan ΛΓ γ c = (30) ∂ α cos Λ ( ) ∂ α c = tan Λ γ − Γ sec Λ + 2 tan ΛΓ γ (31) ∂ β ( ) 2 2 z sec ΛΓ − tan Λ γ + 2 tan ΛΓ γ ∂ α y ( 1 + 2 tan ΛΓ γ ) a c a = − − (32) ∂ p V cos Λ V ∞ ∞ ( ) 2 2 z γ 1 + 2 tan ΛΓ ∂ α x ( 1 + 2 tan ΛΓ γ ) c a a = + (33) ∂ q V cos Λ V ∞ ∞ ( ) ( ) 2 2 2 2 x Γ sec Λ − tan Λ γ + 2 tan ΛΓ γ y γ 1 + 2 tan ΛΓ ∂ α c a a = − (34) ∂ r V V ∞ ∞ ∂ α γ c = − (35) 2 2 ∂ α cos Λ ∂ α c = − Γ ( tan Λ + Γ γ ) (36) ∂ β 2 2 ∂ α y γ z Γ ( tan Λ + Γ γ ) y z ( tan Λ + 2 Γ γ ) c a a a a = − − − (37) 2 2 2 2 2 ∂ p V cos Λ V V cos Λ ∞ ∞ ∞ 2 2 ∂ α x γ z tan ΛΓ ( 1 − tan ΛΓ γ ) x z ( 1 − 2 tan ΛΓ γ ) c a a a a = − + + (38) 2 2 2 2 2 ∂ q V cos Λ V V cos Λ ∞ ∞ ∞ ( ) 2 2 ∂ α x Γ ( tan Λ + Γ γ ) y tan ΛΓ ( 1 − tan ΛΓ γ ) x y Γ 1 − tan Λ − 2 tan ΛΓ γ a a c a a = − + − (39) 2 2 2 2 ∂ r V V V ∞ ∞ ∞ ∂ α tan Λ + 2 Γ γ c = (40) ∂ αβ cos Λ ∂ α 2 y γ z ( tan Λ + 2 Γ γ ) c a a = + (41) ∂ α p V cos Λ V cos Λ ∞ ∞ ∂ α 2 x γ z ( 1 − 2 tan ΛΓ γ ) c a a = − + (42) ∂ α q V cos Λ V cos Λ ∞ ∞ ∂ α x ( tan Λ + 2 Γ γ ) y ( 1 − 2 tan ΛΓ γ ) c a a = − − (43) ∂ α r V cos Λ V cos Λ ∞ ∞ ∂ α y ( tan Λ + 2 Γ γ ) 2 z Γ ( tan Λ + Γ γ ) c a a = − − (44) ∂ β p V cos Λ V ∞ ∞ ( ) z Γ 1 − tan Λ − 2 tan ΛΓ γ ∂ α x ( tan Λ + 2 Γ γ ) c a a = − (45) ∂ β q V cos Λ V ∞ ∞ ( ) y Γ 1 − tan Λ − 2 tan ΛΓ γ ∂ α 2 x Γ ( tan Λ + Γ γ ) c a a = + (46) ∂ β r V V ∞ ∞ ( ) 2 2 z Γ 1 − tan Λ − 2 tan ΛΓ γ ∂ α 2 x y γ x z ( tan Λ + 2 Γ γ ) y z ( 1 − 2 tan ΛΓ γ ) c a a a a a a a = − + + − (47) 2 2 2 2 2 ∂ pq V V cos Λ V cos Λ V cos Λ ∞ ∞ ∞ ∞ ( ) 2 2 y z Γ 1 − tan Λ − 2 tan ΛΓ γ ∂ α y ( 1 − 2 tan ΛΓ γ ) x y ( tan Λ + 2 Γ γ ) 2 x z Γ ( tan Λ + Γ γ ) a a c a a a a a = + + + (48) 2 2 2 2 ∂ pr V cos Λ V cos Λ V V ∞ ∞ ∞ ∞ ( ) ∂ α x ( tan Λ + 2 Γ γ ) x y ( 1 − 2 tan ΛΓ γ ) x z Γ 1 − tan Λ − 2 tan ΛΓ γ a a c a a a = − − + 2 2 2 ∂ qr V cos Λ V cos Λ V ∞ ∞ ∞ 2 y z tan ΛΓ ( 1 − tan ΛΓ γ ) a a − (49) V ∞ ∂ α ( V + ry − qz ) sin Λ + ( V β − rx + pz ) cos Λ + ( V α + qx − py ) Γ c ∞ a a ∞ a a ∞ a a = − × ∂ W V cos Λ x ∞ [ ] ( V + ry − qz ) sin ΛΓ γ + ( V β − rx + pz ) cos ΛΓ γ − ( V α + qx − py ) γ ∞ a a ∞ a a ∞ a a × 1 + (50) V cos Λ ∞ 11 of 38 American Institute of Aeronautics and Astronautics ∂ α 1 ( V + ry − qz ) sin ΛΓ γ + ( V β − rx + pz ) cos ΛΓ γ − ( V α + qx − py ) γ c ∞ a a ∞ a a ∞ a a = − − (51) 2 2 ∂ W V cos Λ V cos Λ t ∞ ∞ ∂ α ( V + ry − qz ) sin Λ + ( V β − rx + pz ) cos Λ + ( V α + qx − py ) Γ c ∞ a a ∞ a a ∞ a a = × ∂ V V cos Λ x ∞ [ ] ( V + ry − qz ) sin ΛΓ + ( V β − rx + pz ) cos ΛΓ − ( V α + qx − py ) ∞ a a ∞ a a ∞ a a × γ − (52) V cos Λ ∞ ∂ α γ ( V + ry − qz ) sin ΛΓ + ( V β − rx + pz ) cos ΛΓ − ( V α + qx − py ) c ∞ a a ∞ a a ∞ a a = − (53) 2 2 ∂ V V cos Λ V cos Λ t ∞ ∞ ∂ α ( V + ry − qz ) cos Λ − ( V β − rx + pz ) sin Λ c ∞ a a ∞ a a = − ∂ Θ V cos Λ ∞ [( V + ry − qz ) sin ΛΓ + ( V β − rx + pz ) cos ΛΓ − ( V α + qx − py )] ∞ a a ∞ a a ∞ a a − (54) 2 2 V cos Λ ∞ ∂ α y y γ [( V + ry − qz ) sin ΛΓ + ( V β − rx + pz ) cos ΛΓ − ( V α + qx − py )] c ∞ a a ∞ a a ∞ a a = − − (55) 2 2 ∂ Θ V cos Λ V cos Λ t ∞ ∞ These partial derivatives contribute to the aeroelastic angle of attack as follows: ∂ α ∂ α ∂ α ∂ α ∂ α ∂ α d α c c c c c c c α ( x , y , z ) = α + s + W + W + V + V + Θ + Θ (56) c 0 x t x t t ∂ s ∂ W ∂ W ∂ V ∂ V ∂ Θ d Θ x t x t t [ ] > 2 2 2 2 2 where s = is a α β p q r α β p q r αβ α p α q α r β p β q β r pq pr qr vector of the aircraft flight dynamic state variables.

The terms W , V , and Θ contribute the aerodynamic stiffness, while the terms W , V , and Θ contribute to the x x t t t aerodynamic damping.

The local angle of attack of an airfoil section at the aerodynamic center is evaluated at x = x , y = y , z = z , a ac a ac a ac y = − e , and z = 0, where x is the forward distance of the aircraft center of gravity from the aerodynamic center of a ac wing section and e is the forward distance of the aerodynamic center from the elastic axis. Then ∂ α ∂ α ∂ α ∂ α ∂ α ac ac ac ac ac α = α + s + U + U + Θ + Θ (57) ac 0 x t t ∂ s ∂ U ∂ U ∂ Θ ∂ Θ x t t There is another source of lift, called non-circulatory lift due to the reaction force of the air volume surrounding a wing section. The non-circulatory lift force is based on the aeroelastic angle of attack at the mid-chord location which is evaluated at x = x , y = y , z = z , y = e , and z = 0, where x is the forward distance of the aircraft center a mc a mc a mc m mc of gravity from the mid-chord location of a wing section and e is the forward distance of the aeroelastic center from m the mid-chord location. Then ∂ α ∂ α ∂ α ∂ α ∂ α mc mc mc mc mc α = α + s + U + U + Θ + Θ (58) mc 0 x t t ∂ s ∂ U ∂ U ∂ Θ ∂ Θ x t t D. Aeroelastic Equations of Coupled Bending-Torsion Motion The equilibrium conditions for bending and torsion are expressed as ∂ M x = − m (59) x ∂ x ∂ M ∂ m y y = f − (60) z ∂ x ∂ x ∂ M ∂ m z z = f − (61) y ∂ x ∂ x where m is the pitching moment per unit span about the elastic axis, f and f are the lift and drag forces per unit x z y span, respectively, and m and m are the bending moments per unit span about the flapwise and chordwise axes of the y z wing.

The wing section lift coefficient is given by c = c C ( k ) α + c δ (62) L L ac L ac α δ 12 of 38 American Institute of Aeronautics and Astronautics ω c where k = is the reduced frequency parameter, ω is the frequency of wing oscillations, c is the section chord, 2 V ∞ c is the section lift curve slope, c is a vector of the section lift derivatives due to the VCCTEF deflection δ = L L α δ [ ] > .

δ δ . . . δ 1 2 14 The function C ( k ) is the Theodorsen’s complex-valued function which is also expressed in terms of Bessel functions as C ( k ) = F ( k ) + iG ( k ) (63) where F ( k ) > 0 and G ( k ) < 0 are shown in Fig. 9.

When k = 0, the airfoil motion is steady and C ( k ) is real and unity.

Fig. 9 -Theodorsen’s Function The wing section lift coefficient due to harmonic motions is expressed as ˙ α c G ( k ) ac c = c α F ( k ) + c + c δ (64) L L ac L L ac α α δ 2 V k ∞ In addition, the apparent mass of the air contributes to the lift force acting at the mid-chord location as follows: ˙ π α c mc c = (65) L mc 2 V ∞ The total section lift coefficient is c = c + c (66) L L L ac mc The section pitching moment coefficient is evaluated as e e m c = c + c − c + c δ (67) m m L L m ac ac mc δ c c where c is the section pitching moment coefficient at the aerodynamic center and c is a vector of the section m m ac δ pitching moment derivatives at the elastic axis due to the VCCTEF deflection.

The section drag coefficient is expressed in a parabolic drag polar form as c L ac c = c + (68) D D π AR ε where c is the section parasitic drag coefficient, AR is the wing aspect ratio, and ε is the span efficiency factor.

D 13 of 38 American Institute of Aeronautics and Astronautics Expanding the expression and ignoring the nonlinear terms for the aeroelastic analysis, the section drag coefficient is obtained as ˙ α c G ( k ) ac c = c + c α F ( k ) + c + c δ (69) D D D ac D D 0 α α δ 2 V k ∞ where 2 c α L α c = (70) D α π AR ε 2 c c α L L 0 α δ c = (71) D δ π AR ε In addition, the propulsive effects of the aircraft engines must be accounted for in the analysis. Both the engine mass and thrust can contribute to the aeroelasticity. The propulsive force and moment vector are computed as       − sin Λ − cos Λ sin ΛΓ d ( − T sin Λ − m g Γ ) d 1 e 1 [ ]       f = δ ( x − x ) = δ ( x − x ) T 0 m g  − cos Λ sin Λ cos ΛΓ   d   − T cos Λ d  e e e 2 e 2 − Γ 0 − 1 d ( T sin ΛΓ − m g ) d 3 e 3 (72)   ( − Ty sin ΛΓ − T z cos Λ + m gy ) d e e e e 1   m = r × f = ( x d − y d − z d ) × f = δ ( x − x ) (73) e e e e 1 e 2 e 3 e e  ( − T x sin ΛΓ + m gx + T z sin Λ + m gz Γ ) d  e e e e e e 2 ( − T x cos Λ − Ty sin Λ − m gy Γ ) d e e e e 3 where T is the engine thrust, m is the engine mass, ( x , y , z ) is the coordinate of the engine thrust center such that y e e e e e is positive forward of the elastic axis and z is positive below the elastic axis, and δ ( x − x ) is the Dirac delta function e e such that ∫ δ ( x − x ) f ( x ) dx = f ( x ) (74) e e Transforming into the local coordinate reference frame and neglecting nonlinear contributions, the propulsive forces and moments are given by f = δ ( x − x ) [ − T sin Λ − m g Γ − T cos Λ V + ( T sin ΛΓ − m g ) W ] (75) x e e x e x e f = δ ( x − x ) [( T sin Λ + m g Γ ) V − T cos Λ + ( T sin ΛΓ − m g ) ( Θ + γ )] (76) y e e x e e f = δ ( x − x ) [( T sin Λ + m g Γ ) W + T cos Λ ( Θ + γ ) + T sin ΛΓ − m g ] (77) z e e x e e m = δ ( x − x ) [ − Ty sin ΛΓ − T z cos Λ + m gy + ( − T x sin ΛΓ + m gx + T z sin Λ + m gz Γ ) V x e e e e e e e e e e e x e − ( T x cos Λ + Ty sin Λ + m gy Γ ) W ] (78) e e e e x m = δ ( x − x ) [( Ty sin ΛΓ + T z cos Λ − m gy ) V − T x sin ΛΓ + m gx + T z sin Λ + m gz Γ y e e e e e x e e e e e e e − ( T x cos Λ + Ty sin Λ + m gy Γ ) ( Θ + γ )] (79) e e e e m = δ ( x − x ) [( Ty sin ΛΓ + T z cos Λ − m gy ) W − ( − T x sin ΛΓ + m gx + T z sin Λ + m gz Γ ) ( Θ + γ ) z e e e e e x e e e e e e e − T x cos Λ − Ty sin Λ − m gy Γ ] (80) e e e e The partial derivatives of the moment components are e ∂ m x = δ ( x − x ) [( − T x sin ΛΓ + m gx + T z sin Λ + m gz Γ ) V − ( T x cos Λ + Ty sin Λ + m gy Γ ) W ] (81) e e e e e e e xx e e e e xx ∂ x e [ ( )] ∂ m ′ y = δ ( x − x ) ( Ty sin ΛΓ + T z cos Λ − m gy ) V − ( T x cos Λ + Ty sin Λ + m gy Γ ) Θ + γ (82) e e e e e xx e e e e x ∂ x e [ ( )] ∂ m ′ z = δ ( x − x ) ( Ty sin ΛΓ + T z cos Λ − m gy ) W + ( T x sin ΛΓ − m gx − T z sin Λ − m gz Γ ) Θ + γ e e e e e xx e e e e e e x ∂ x (83) Using the sign convention as shown in Fig. 10 the lift and drag forces and pitching moment per unit span can be expressed as 14 of 38 American Institute of Aeronautics and Astronautics Fig. 10 - Airfoil Forces and Moment [ ] c G ( k ) π ce m 2 2 m = − cc + ec + ec α F ( k ) + ec ˙ α − ˙ α q cos Λ c + mge − mk Θ x m L L ac L ac mc ∞ cg tt ac 0 α α 2 V k 2 V ∞ ∞ [ ( ) ] 2 2 + me W − δ ( x − x ) m y + z Θ − m y W + m z V + δ ( x − x ) [ − Ty sin ΛΓ − T z cos Λ + m gy cg tt e e tt e e tt e e tt e e e e e e e + ( − T x sin ΛΓ + m gx + T z sin Λ + m gz Γ ) V − ( T x cos Λ + Ty sin Λ + m gy Γ ) W ] (84) e e e e e e x e e e e x [ ] c G ( k ) ˙ f = c + c + c α F ( k ) + c α q cos Λ c − mV + δ ( x − x ) ( − m V − m z Θ ) y D D D ac D ac ∞ tt e e tt e e tt α α 0 i 2 V k ∞ + δ ( x − x ) [( T sin Λ + m g Γ ) V − T cos Λ + ( T sin ΛΓ − m g ) ( Θ + γ )] (85) e e x e [ ] c G ( k ) π c f = c + c α F ( k ) + c ˙ α + ˙ α q cos Λ c − mg − mW + me Θ z L L ac L ac mc ∞ tt cg tt 0 α α 2 V k 2 V ∞ ∞ + δ ( x − x ) ( − m W + m y Θ ) + δ ( x − x ) [( T sin Λ + m g Γ ) W + T cos Λ ( Θ + γ ) + T sin ΛΓ − m g ] (86) e e tt e e tt e e x e where q is the dynamic pressure, m is the wing mass distribution, e is the eccentricity between the center of mass ∞ cg and the elastic axis (positive corresponding to the center of mass located forward of the elastic axis), k is the torsional radius of gyration, and the term cos Λ accounts for the wing sweep angle Λ as measured from the elastic axis.

The bending and torsion aeroelastic equations then become {[ ] } ( ) ∂ ′ ′ ′ GJ + EB γ Θ − EB γ W − EB γ V = 1 x 2 xx 3 xx ∂ x [ ] ( ) c G ( k ) π ce m 2 2 ˙ ˙ cc + ec α F ( k ) + ec α − α + ec + cc δ q cos Λ c − mge + mk Θ m L ac L ac mc L m ∞ cg tt ac α α δ δ 2 V k 2 V ∞ ∞ [ ( ) 2 2 − me W + δ ( x − x ) m y + z Θ − m y W + m z V + Ty sin ΛΓ + T z cos Λ − m gy cg tt e e tt e e tt e e tt e e e e e e ] − ( − T x sin ΛΓ + m gx + T z sin Λ + m gz Γ ) V + ( T x cos Λ + Ty sin Λ + m gy Γ ) W (87) e e e e e e x e e e e x ( ) ∂ ′ − EB γ Θ + EI W − EI V = 2 x yy xx yz xx ∂ x [ ] c G ( k ) π c ˙ ˙ c α F ( k ) + c α + α + c δ q cos Λ c − mg − mW + me Θ L ac L ac mc L ∞ tt cg tt α α δ 2 V k 2 V ∞ ∞ [ + δ ( x − x ) − m W + m y Θ + ( T sin Λ + m g Γ ) W + T cos Λ ( Θ + γ ) + T sin ΛΓ − m g e e tt e e tt e x e ( )] ′ − ( Ty sin ΛΓ + T z cos Λ − m gy ) V + ( T x cos Λ + Ty sin Λ + m gy Γ ) Θ + γ (88) e e e e xx e e e e x 15 of 38 American Institute of Aeronautics and Astronautics ( ) ∂ ′ − EB γ Θ − EI W + EI V = 3 x yz xx zz xx ∂ x [ ] c G ( k ) ˙ c + c α F ( k ) + c α + c δ q cos Λ c − mV D D ac D ac D ∞ tt α α 0 δ 2 V k ∞ [ + δ ( x − x ) − m V − m z Θ + ( T sin Λ + m g Γ ) V − T cos Λ + ( T sin ΛΓ − m g ) ( Θ + γ ) e e tt e e tt e x e ( )] ′ − ( Ty sin ΛΓ + T z cos Λ − m gy ) W − ( T x sin ΛΓ − m gx − T z sin Λ − m gz Γ ) Θ + γ (89) e e e e xx e e e e e e x Taking advantage of symmetry, the motion could be decomposed into symmetric and anti-symmetric motions subject to either symmetric-mode boundary conditions W ( 0 , t ) = V ( 0 , t ) = 0 or anti-symmetric mode boundary x x conditions Θ ( 0 , t ) = W ( 0 , t ) = V ( 0 , t ) = 0 at the left end. Half of the mass and mass inertia of the aircraft structure without the wings are added to the generalized mass of the system. Defining the vector quantities [ ] W U = (90) V [ ] c = (91) c D [ ] c L α c = (92) α c D α [ ] c L δ c = (93) δ c D δ [ ] π c 2 V ∞ c = (94) c [ ] g a = (95) [ ] e cg ε = (96) cg [ ] y e r = (97) e − z e [ ] T sin ΛΓ − m g e f = (98) − T cos Λ [ ] T cos Λ f = (99) Θ T sin ΛΓ − m g e [ ] T x cos Λ + Ty sin Λ + m gy Γ e e e e f = (100) Θ x − ( T x sin ΛΓ − m gx − T z sin Λ − m gz Γ ) e e e e e e [ ] > T x cos Λ + Ty sin Λ + m gy Γ e e e e f = (101) U x − ( − T x sin ΛΓ + m gx + T z sin Λ + m gz Γ ) e e e e e e [ ] B = (102) B B 2 3 [ ] I − I yy yz I = (103) − I I zy zz [ ] 0 1 J = (104) 1 0 16 of 38 American Institute of Aeronautics and Astronautics The aeroelastic partial differential equations are given as {[ ] } ( ) ∂ ′ ′ GJ + EB γ Θ − EB γ U = 1 x xx ∂ x [ ] ( ) c G ( k ) π ce m cc + ec + ec α F ( k ) + ec ˙ α − ˙ α + ec + cc δ q cos Λ c m L L ac L ac mc L m ∞ ac 0 α α δ δ 2 V k 2 V ∞ ∞ [ ( ) 2 2 2 > − mge + mk Θ − m ε U + δ ( x − x ) m y + z Θ − m r U cg tt cg tt e e tt e tt e e e + Ty sin ΛΓ + T z cos Λ − m gy + f U ] (105) e e e e U x x ( ) ∂ ′ > − EB γ Θ + EIU = x xx ∂ x [ ] c G ( k ) c + c α F ( k ) + c ˙ α + c ˙ α + c δ q cos Λ c − ma 0 α ac α ac c mc δ ∞ 2 V k ∞ − mU + m ε Θ + δ ( x − x ) [ − m U + m r Θ + ( T sin Λ + m g Γ ) U tt cg tt e e tt e e tt e x ( )] ′ + f − ( Ty sin ΛΓ + T z cos Λ − m gy ) JU + f ( Θ + γ ) + f Θ + γ (106) 0 e e e e xx Θ Θ x x IV. Finite Element Modeling A. Finite-Element Discretization The aeroelastic equations describe the wing bending and torsional deflections due to aerodynamic forces and moments.

Using the finite-element method, the structure can be discretized into n equally spaced one-dimensional elements.

Then the bending and torsional deflections can be approximated as n Θ ( x , t ) = Θ ( x , t ) (107) i

∑

i = 1 n U ( x , t ) = U ( x , t ) (108) i

∑

i = 1 where i is the i -th element, and n is the number of nodes.

For each element, the bending and torsional deflections are approximated as [ ] [ ] θ ( t ) i Θ ( x , t ) = = N ( x ) θ ( t ) (109) ψ ( x ) ψ ( x ) i θ i 1 2 θ ( t ) i   w ( t ) i  ′   w ( t )  i     v ( t ) i [ ]   ′   φ ( x ) φ ( x ) 0 0 φ ( x ) φ ( x ) 0 0 v ( t )   1 2 3 4 1 i U ( x , t ) = = N ( x ) u ( t ) (110)   i u i   0 0 φ ( x ) φ ( x ) 0 0 φ ( x ) φ ( x ) w ( t ) 1 2 3 4 2 i   ′   w ( t )  2  i    v ( t )  i ′ v ( t ) i where the subscripts 1 and 2 denote values at nodes 1 and 2, and ψ ( x ) , j = 1 , 2, and φ ( x ) , k = 1 , 2 , 3 , 4 are the linear j k and Hermite polynomial shape functions x ψ ( x ) = 1 − (111) l x ψ ( x ) = (112) l 17 of 38 American Institute of Aeronautics and Astronautics ( ) ( ) 2 3 x x φ ( x ) = 1 − 3 + 2 (113) l l [ ] ( ) ( ) 2 3 x x x φ ( x ) = l − 2 + (114) l l l ( ) ( ) 2 3 x x φ ( x ) = 3 − 2 (115) l l [ ] ( ) ( ) 2 3 x x φ ( x ) = l − + (116) l l L where x ∈ [ 0 , l ] is the local coordinate and l = is the element length.

n − 1 It can be shown that Hermite cubic shape functions result in exact nodal displacements and slopes, making them ideal candidates for problems involving beam bending.

The weak-form integral expressions of the dynamic aeroelastic equations are obtained by multiplying the bending > > and torsion aeroelastic equations by N ( x ) and N ( x ) , and then integrating over the wing span. This yields θ u l {[ ] } ∫ n ( ) d ′ ′ ′ ′′ > N GJ + EB γ N θ − EB γ N u dx = 1 i i

∑ θ θ u

dx i = 1 l [ ( ) ∫ n ∂ α ∂ α ′ ∂ α ∂ α ∂ α ac ac ac ac ac > ˙ N cc + ec α + s + N u + N ˙ u + N θ + N θ F ( k ) m L 0 i u i θ i θ i

∑ θ ac α u

∂ s ∂ U ∂ U ∂ Θ ∂ Θ x t t i = 1 ( ) ∂ α ∂ α ′ ∂ α ∂ α ∂ α c G ( k ) ac ac ac ac ac ˙ ¨ + ec ˙ s + N ˙ u + N ¨ u + N θ + N θ L i u i θ i θ i α u ∂ s ∂ U ∂ U ∂ Θ ∂ Θ 2 V k x t t ∞ ( ) ] ( ) π ce ∂ α ∂ α ′ ∂ α ∂ α ∂ α m mc mc mc mc mc ˙ ¨ − ˙ s + N ˙ u + N ¨ u + N θ + N θ + ec + cc δ q cos Λ cdx u i u i θ i θ i L m ∞ δ δ 2 V ∂ s ∂ U ∂ U ∂ Θ ∂ Θ ∞ x t t l ∫ n ( ) > 2 ¨ + N − mge + mk N θ − m ε N ¨ u dx cg i cg u i

∑ θ

θ i = 1 n [ ] ( ) ′ > 2 2 > ¨ + N m y + z N θ − m r N ¨ u + Ty sin ΛΓ + T z cos Λ − m gy + f N u (117) e i e u i e e e e U i

∑ e e θ e u

θ x x = x e i = 1 l ∫ n 2 ( ) d ′ ′ ′′ > > N − EB γ N θ + EIN u dx = i i

∑ u θ u

dx i = 1 l [ ( ) ∫ n ∂ α ∂ α ′ ∂ α ∂ α ∂ α ac ac ac ac ac > ˙ N c + c α + s + N u + N ˙ u + N θ + N θ F ( k ) 0 α 0 i u i θ i θ i

∑ u u

∂ s ∂ U ∂ U ∂ Θ ∂ Θ x t t i = 1 ( ) ∂ α ∂ α ′ ∂ α ∂ α ∂ α c G ( k ) ac ac ac ac ac ˙ ¨ + c ˙ s + N ˙ u + N ¨ u + N θ + N θ α u i u i θ i θ i ∂ s ∂ U ∂ U ∂ Θ ∂ Θ 2 V k x t t ∞ ( ) ] ∂ α ∂ α ′ ∂ α ∂ α ∂ α mc mc mc mc mc ˙ ¨ + c ˙ s + N ˙ u + N ¨ u + N θ + N θ + c δ q cos Λ cdx c i u i i i ∞ u θ θ δ ∂ s ∂ U ∂ U ∂ Θ ∂ Θ x t t l ∫ n n [ ( ) ′ > ¨ ¨ + N − ma − mN ¨ u + m ε N θ dx + − m N ¨ u + m r N θ + ( T sin Λ + m g Γ ) N u u i cg i e u i e e i e i

∑ u θ ∑ θ u

i = 1 i = 1 ( )] ′′ ′ ′ + f − ( Ty sin ΛΓ + T z cos Λ − m gy ) JN u + f ( N θ + γ ) + f N θ + γ (118) 0 e e e e i Θ θ i Θ i u x θ x = x e The expressions of the left hand sides can be integrated by parts upon enforcing the boundary conditions as l l {[ ] } {[ ] } ∫ ∫ ( ) ( ) 2 2 d ′ ′ ′ ′′ ′ ′ ′ ′ ′′ > > N GJ + EB γ N θ − EB γ N u dx = − N GJ + EB γ N θ − EB γ N u dx (119) 1 i i 1 i i θ θ u θ θ u dx 0 0 18 of 38 American Institute of Aeronautics and Astronautics l l ∫ ∫ ( ) ( ) d ′ ′ ′′ ′′ ′ ′ ′′ > > > > N − EB γ N θ + EIN u dx = N − EB γ N θ + EIN u dx (120) i i i i u θ u u θ u dx 0 0 Upon substitution, the final form of the aeroelastic equations are produced as l {[ ] } ∫ n ( ) ′ ′ ′ ′ ′′ > N GJ + EB γ N θ − EB γ N u dx 1 i i

∑ θ θ u

i = 1 l [ ( ) ∫ n ∂ α ∂ α ′ ∂ α ∂ α ∂ α ac ac ac ac ac > ˙ + N cc + ec α + s + N u + N ˙ u + N θ + N θ F ( k ) m L 0 i u i θ i θ i

∑ θ ac α u

∂ s ∂ U ∂ U ∂ Θ ∂ Θ x t t i = 1 ( ) ∂ α ∂ α ′ ∂ α ∂ α ∂ α c G ( k ) ac ac ac ac ac ˙ ¨ + ec ˙ s + N ˙ u + N ¨ u + N θ + N θ L i u i θ i θ i α u ∂ s ∂ U ∂ U ∂ Θ ∂ Θ 2 V k x t t ∞ ( ) ] ( ) π ce ∂ α ∂ α ′ ∂ α ∂ α ∂ α m mc mc mc mc mc ˙ ¨ − ˙ s + N ˙ u + N ¨ u + N θ + N θ + ec + cc δ q cos Λ cdx i u i i i L m ∞ u θ θ δ δ 2 V ∂ s ∂ U ∂ U ∂ Θ ∂ Θ ∞ x t t l ∫ n ( ) > 2 ¨ + N − mge + mk N θ − m ε N ¨ u dx cg i cg u i

∑ θ

θ i = 1 n [ ] ( ) ′ > 2 2 > ¨ + N m y + z N θ − m r N ¨ u + Ty sin ΛΓ + T z cos Λ − m gy + f N u = 0 (121) e i e u i e e e e U i

∑ e e θ e u

θ x x = x e i = 1 l ∫ n ( ) ′′ ′ ′ ′′ > > N − EB γ N θ + EIN u dx i i

∑ u θ u

i = 1 l [ ( ) ∫ n ∂ α ∂ α ′ ∂ α ∂ α ∂ α ac ac ac ac ac > ˙ − N c + c α + s + N u + N ˙ u + N θ + N θ F ( k ) 0 α 0 i u i θ i θ i

∑ u u

∂ s ∂ U ∂ U ∂ Θ ∂ Θ x t t i = 1 ( ) ∂ α ∂ α ′ ∂ α ∂ α ∂ α c G ( k ) ac ac ac ac ac ˙ ¨ + c ˙ s + N ˙ u + N ¨ u + N θ + N θ α u i u i θ i θ i ∂ s ∂ U ∂ U ∂ Θ ∂ Θ 2 V k x t t ∞ ( ) ] ∂ α ∂ α ′ ∂ α ∂ α ∂ α mc mc mc mc mc ˙ ¨ + c ˙ s + N ˙ u + N ¨ u + N θ + N θ + c δ q cos Λ cdx c i u i i i ∞ u θ θ δ ∂ s ∂ U ∂ U ∂ Θ ∂ Θ x t t l ∫ n n [ ( ) ′ > > ¨ ¨ − N − ma − mN ¨ u + m ε N θ dx − N − m N ¨ u + m r N θ + ( T sin Λ + m g Γ ) N u u i cg i e u i e e i e i

∑ u θ ∑ u θ u

i = 1 i = 1 ( )] ′′ ′ ′ + f − ( Ty sin ΛΓ + T z cos Λ − m gy ) JN u + f ( N θ + γ ) + f N θ + γ = 0 (122) 0 e e e e i Θ θ i Θ i u x θ x = x e The mass matrix, damping matrix, stiffness matrix, and force vector corresponding to each finite element i are then established as [ ] [ ] l ( ) ∫ > 2 > > 2 2 > > N k N − N ε N N y + z N − N r N θ cg u θ u e e e θ θ θ θ M = m dx + m i e > > > > − N ε N N N − N r N N N cg θ u e θ u u u u u x = x 0 e  ( ) ( )  l G ( k ) G ( k ) ∂ α c π ce ∂ α ∂ α c π ce ∂ α ∫ > ac m mc > ac m mc N ec − N N ec − N L θ L u θ α θ α ∂ Θ 2 V k 2 V ∂ Θ ∂ U 2 V k 2 V ∂ U t ∞ ∞ t t ∞ ∞ t 2  ( ) ( )  + q cos Λ cdx (123) ∞ G ( k ) G ( k ) ∂ α c ∂ α ∂ α c ∂ α > ac mc > ac mc − N c + c N − N c + c N α c θ α c u u u ∂ Θ 2 V k ∂ Θ ∂ U 2 V k ∂ U 0 t ∞ t t ∞ t 19 of 38 American Institute of Aeronautics and Astronautics  ( )  l ∂ α ∂ α G ( k ) ∂ α ∫ > ac ac c > ac N ec F ( k ) + N N ec F ( k ) N L θ L u α α θ ∂ Θ ∂ Θ 2 V k θ ∂ U t ∞ t 2  ( )  C = q cos Λ cdx i ∞ ∂ α ∂ α G ( k ) ∂ α > ac ac c > ac − N c F ( k ) + N − N c F ( k ) N α θ α u u u ∂ Θ ∂ Θ 2 V k ∂ U t ∞ t [ ] [ ] l l ′ ′ ∫ G ( k ) ∫ ∂ α c π ce ∂ α π ce ∂ α > ac > m mc > m mc 0 N ec N − N N − N N L θ θ α u 2 θ θ u 2 ∂ U 2 V k 2 V ∂ Θ 2 V ∂ U x ∞ ∞ ∞ x + q cos Λ cdx + q cos Λ cdx ′ ∞ ′ ∞ G ( k ) ∂ α c ∂ α ∂ α > ac > mc > mc 0 − N c N − N c N − N c N α c θ c u u u u u ∂ U 2 V k ∂ Θ ∂ U x ∞ x 0 0 (124) [ ]   ( ) l ′ ′ ′ ′ ′ ′′ ∫ > > N GJ + EB γ N − N EB γ N u θ θ θ   K = dx i ′′ ′ ′ ′′ ′′ > > > − N EB γ N N EIN u u u θ [ ] ′ > 0 N f N U x u θ ( ) [ ] + ′ ′ ′′ > > − N f N + f N − N ( T sin Λ + m g Γ ) N − ( Ty sin ΛΓ + T z cos Λ − m gy ) JN Θ θ Θ e e e e e u x u u u θ x = x e [ ] l ′ ∫ > ∂ α > ∂ α ac ac N ec F ( k ) N N ec F ( k ) N L θ L u θ α θ α 2 ∂ Θ ∂ U x + q cos Λ cdx (125) ′ ∞ ∂ α ∂ α > ac > ac − N c F ( k ) N − N c F ( k ) N α θ α u u u ∂ Θ ∂ U x [ ] [ ] l ∫ > > − N ( Ty sin ΛΓ + T z cos Λ − m gy ) e e e e N ge cg θ ( ) θ F = m dx + i ′ > > N f + f γ + f γ − N a 0 u Θ Θ u x 0 x = x e [ ] l ∫ > − N ( cc + ec α ) m L 0 ac α 2 θ + q cos Λ cdx (126) ∞ > N ( c + c α ) 0 α 0 u [ ] l ∫ > ∂ α ac ∂ F − N ec F ( k ) i L 2 θ α ∂ s = q cos Λ cdx (127) ∞ ∂ α > ac ∂ s N c F ( k ) L u α ∂ s  ( )  l G ( k ) > ∂ α c π ce ∂ α ∫ ac m mc − N ec − ∂ F L θ α 2 V k 2 V i ∂ s ∂ s ∞ ∞ 2  ( )  = q cos Λ cdx (128) ∞ G ( k ) > ∂ α c ∂ α ac mc ∂ ˙ s N c + c α c u 2 V k ∂ s ∞ ∂ s [ ] l ( ) ∫ > ∂ F − N ec + cc i L m θ δ δ = q cos Λ cdx (129) ∞ > ∂ δ N c u δ The discrete global approximation system is formed by enforcing equilibrium conditions at element interfaces then summing each element matrix during the assembly process, resulting in n M = M (130) i

∑

i = 1 n C = C (131) i

∑

i = 1 n K = K (132)

∑ i

i = 1 ( ) n ∂ F ∂ F ∂ F i i i F = F + s + ˙ s + δ (133)

∑ i

∂ s ∂ ˙ s ∂ δ i = 1 The resultant matrix equation is of the form ∂ F ( k ) ∂ F ( k ) ∂ F M ( k ) ¨ x + C ( k ) ˙ x + K ( k ) x = F + s + ˙ s + δ (134) e e e ∂ s ∂ ˙ s ∂ δ [ ] > ′ ′ ′ ′ where x = is the nodal displacement vector, e θ w w v v · · · w θ w w v v 1 1 1 n n + 1 n + 1 n + 1 1 1 n + 1 n + 1 M is the mass matrix, C is the damping matrix, K is the stiffness, and F is the force vector, all as functions of the reduced frequency parameter k .

20 of 38 American Institute of Aeronautics and Astronautics B. Structural Damping It is important to note that the aerodynamic damping matrix can be augmented by a structural damping matrix. The structural damping matrix can be obtained by transforming the generalized coordinates into modal coordinates via the eigenvalue analysis as follows: Consider the zero-speed structural dynamic equations ( ) ∂ F ∂ F ∂ F − 1 − 1 − 1 ¨ x + M C ˙ x + M K x = M F + s + ˙ s + δ (135) e s e s e s s s ∂ s ∂ ˙ s ∂ δ where M is the structural mass matrix, C is the structural damping matrix, K is the structural stiffness matrix at zero s s s − 1 speed ( V = 0), and F is the force vector. Assuming that the eigenvalues of the matrix M K are positive real and ∞ s s − 1 distinct, the matrix M K may be simplified using the similarity transformation s s − 1 2 − 1 M K = X Ω X (136) s s s s s where X is the eigenvector matrix and Ω = diag ( ω , ω , . . . , ω ) is a diagonal matrix whose entries are the frequencies s s 1 2 n of the structural dynamic modes.

The structural dynamic matrix equation may therefore be expressed as ( ) ∂ F ∂ F ∂ F − 1 2 − 1 − 1 ¨ x + M C ˙ x + X Ω X x = M F + s + ˙ s + δ e s e s e s s s s ∂ s ∂ ˙ s ∂ δ − 1 − 1 Multiplying through by X and letting ϕ = X x be the modal coordinates gives the transformed structural e dynamics equation ( ) ∂ F ∂ F ∂ F − 1 − 1 2 − 1 − 1 ¨ ˙ ϕ + X M C X ϕ + Ω ϕ = X M F + s + ˙ s + δ (137) s s s s s s s ∂ s ∂ ˙ s ∂ δ which can also be expressed in modal coordinates as ∂ f ∂ f ∂ f i i i ¨ ϕ + 2 ζ ω ˙ ϕ + ω ϕ = f + s + ˙ s + δ (138) i i i i i i i ∂ s ∂ ˙ s ∂ δ where ζ is the viscous damping ratio of the i -th mode and is typically a parameter obtained by ground vibration testing i and similar methods.

If ζ = diag ( ζ , ζ , . . . , ζ ) gives the diagonal viscous damping ratio matrix, then the structural damping matrix is 1 2 n computed as − 1 C = 2 M X ζ Ω X (139) s s s s s The total damping matrix is then given by the linear superposition of the structural and aerodynamic damping matrices C = C + C (140) s a where C is the aerodynamic damping matrix computed from the previous section.

a The system of equations is then translated into a state space form [ ] [ ] [ ] [ ] ˙ x 0 I x e e ( ) = + (141) ∂ F ∂ F ∂ F − 1 − 1 − 1 M F + s + ˙ s + δ ¨ x − M K − M C ˙ x e e ∂ s ∂ ˙ s ∂ δ with I representing the identity matrix, whereupon the eigenvalues of the matrix equation yield the vibrational fre- quencies and mode shapes of the aeroelastic system. The flutter boundary is defined to be the airspeed at which the real parts of the eigenvalues of the systems become zero.

C. Time Integration Methods The state space equations may be used to perform time integration of the nodal displacement values given an initial deflection profile. For instance, the first-order, explicit (forward) Euler scheme of integration at time step i + 1 is given by [ ] [ ] ([ ] [ ] [ ]) ˙ x ˙ x 0 I ˙ x e e e ( ) = + ∆ t + (142) − 1 ∂ F ∂ F ∂ F − 1 − 1 ¨ x ¨ x − M K − M C ¨ x M F + s + ˙ s + δ e e e ∂ s ∂ ˙ s ∂ δ i + 1 i i 21 of 38 American Institute of Aeronautics and Astronautics where ∣ ∆ t ≤ ([ ])∣ (143) ∣ ∣ ∣ 0 I ∣ max ∣ eig ∣ − 1 − 1 ∣ ∣ − M K − M C must be satisfied in order to maintain numerical stability of the solution.

This integration scheme can become prohibitively expensive to enforce in problems with large frequency magni- tudes. One approach is to apply modal truncation method so that the time increment may be artificially increased by truncating the number of modes obtained from the analysis. This is due to the fact that it is generally difficult to excite 1 2 higher-frequency modes because the energy in a structure is proportional to Kx , which is in turn proportional to the squared angular frequency ω . However, there are alternative time integration methods which can be employed n to relax such restrictions or increase solution accuracy. One such method is the implicit backward Euler scheme of integration which is similarly expressed as [ ] [ ] ([ ] [ ] [ ]) ˙ x ˙ x 0 I ˙ x e e e ( ) = + ∆ t + (144) ∂ F ∂ F ∂ F − 1 − 1 − 1 M F + s + ˙ s + δ ¨ x ¨ x − M K − M C ¨ x e e e ∂ s ∂ ˙ s ∂ δ i + 1 i i + 1 from which it follows that [ ] ( [ ]) ([ ] [ ]) − 1 ˙ x 0 I ˙ x e e ( ) = I − ∆ t + ∆ t (145) − 1 ∂ F ∂ F ∂ F − 1 − 1 M F + s + ˙ s + δ ¨ x − M K − M C ¨ x e e ∂ s ∂ ˙ s ∂ δ i + 1 i This method is fairly stable and permits larger time steps to be used. It is important to note, however, that the explicit scheme is time-accurate, as compared to the implicit scheme which effectively damps out higher-frequency oscillations. Due to the first-order nature of both the Euler methods, higher-accuracy methods may be desired. One popular solution method, Newmark integration, bypasses the state space form and instead algebraically updates the nodal degrees of freedom of the aeroelastic finite element equations: M ¨ x + C ˙ x + Kx = F (146) Given specified initial values of ¨ x , ˙ x , and x (assumed to be zero in this study), one can freely select a time step ∆ t ( ) 1 1 1 and the parameters δ ≥ and β ≥ δ + . Then, the integration coefficients are defined as 2 4 2 a = (147) β ∆ t δ a = (148) β ∆ t a = (149) β ∆ t a = − 1 (150) 2 β δ a = − 1 (151) β ( ) ∆ t δ a = − 2 (152) 2 β a = ∆ t ( 1 − δ ) (153) a = δ ∆ t (154) Next, the effective stiffness matrix is formed as ˜ K = K + a M + a C (155) 0 1 22 of 38 American Institute of Aeronautics and Astronautics Then, for each time step i , the integration is performed as follows: ˜ F = F + M ( a x + a ˙ x + a ¨ x ) + C ( a x + a ˙ x + a ¨ x ) (156) i i 0 i 2 i 3 i 1 i 4 i 5 i − 1 ˜ ˜ x = K F (157) i + 1 i ¨ x = a ( x − x ) − a ˙ x − a ¨ x (158) i + 1 0 i + 1 i 2 i 3 i ˙ x = ˙ x + a ¨ x + a ¨ x (159) i + 1 i 6 i 7 i + 1 1 1 It can be shown that for δ = and β = , the Newmark time integration method is second-order accurate and 2 4 unconditionally stable. However, due to the numerous integration parameters and multi-step nature of the method, the Newmark scheme can be computationally expensive.

V. Flight Dynamic Coupling Consider the rigid aircraft flight dynamics in the stability axes described by m ( ˙ u + qw − rv ) = ( C sin α − C cos α ) q S + T − m g sin θ (160) a L D ∞ a m ( ˙ v + ru − pw ) = C q S + mg cos θ sin φ (161) Y ∞ m ( ˙ w + pv − qu ) = ( − C cos α − C sin α ) q S + mg cos θ cos φ (162) L D ∞ ( ) 2 2 ¯ ¯ ¯ ¯ ¯ ¯ ¯ ¯ I ˙ p − I ˙ q − I ˙ r − I pq + I pr + ( I − I ) qr + I r − q = C q Sb (163) xx xy xz xz xy zz yy yz l ∞ ( ) 2 2 ¯ ¯ ¯ ¯ ¯ ¯ ¯ ¯ − I ˙ p + I ˙ q − I ˙ r + I pq − I qr + ( I − I ) pr + I p − r = C q S ¯ c + T ¯ z (164) xy yy yz yz xy xx zz xz m ∞ e ( ) 2 2 ¯ ¯ ¯ ¯ ¯ ¯ ¯ ¯ − I ˙ p − I ˙ q + I ˙ r − I pr + I qr + ( I − I ) pq + I q − p = C q Sb (165) xz yz zz yz xz yy xx xy n ∞ ˙ φ = p + q sin φ tan θ + r cos φ tan θ (166) ˙ θ = q cos φ − r sin φ (167) ˙ ψ = q sin φ sec θ + r cos φ sec θ (168) where m is the aircraft mass, S is the aircraft reference wing area, ¯ I , ¯ I , ¯ I , ¯ I , ¯ I , ¯ I are the aircraft principal a xx yy zz xy xz yz moments of inertia, ¯ c is the mean aerodynamic chord, b is the wing span, ¯ z is the offset of the thrust line below the e aircraft CG, and ( φ , θ , ψ ) are the aircraft Euler angles.

Note that the aerodynamic coefficients C , C , C , C , C , and C are influenced by the aeroelastic deflections of L D y l m n the aircraft wings. So, the equations of motion of rigid aircraft are coupled with the aeroelastic equations.

The wing contribution to the aircraft unsteady lift coefficient is evaluated as ( ) ∫ L ˙ ˙ 1 α c G ( k ) π α c ac mc ∆ C ( k ) = c α F ( k ) + c + + c δ cos Λ cdx (169) L L ac L L α α δ S 2 V k 2 V − L ∞ ∞ This can be also written as ( ) n ˙ ¨ ∆ C = C + C s + C ˙ s + C u + C ˙ u + C ¨ u + C θ + C θ + C θ + C δ + C δ (170) L L L L L i L i L i L i L i L i L L e

0 s ˙ s ∑ u ˙ u ¨ u ˙ ¨

θ δ δ θ θ e i = 1 where ∫ L C = c α cos Λ cdx (171) L L 0 α S − L ∫ L 1 ∂ α ac C = c F ( k ) cos Λ cdx (172) L L s α S ∂ s − L ( ) ∫ L 1 c G ( k ) ∂ α π c ∂ α ac mc C = c + cos Λ cdx (173) L L ˙ s α S 2 V k ∂ s 2 V ∂ s − L ∞ ∞ ∫ l 2 ∂ α ′ ac C = c N F ( k ) cos Λ cdx (174) L L u α u S ∂ U x 23 of 38 American Institute of Aeronautics and Astronautics ( ) ∫ L 2 ∂ α c G ( k ) ∂ α ′ π c ∂ α ′ ac ac mc C = c N F ( k ) + c N + N cos Λ cdx (175) L L u L ˙ u α α u u S ∂ U 2 V k ∂ U 2 V ∂ U − L t ∞ x ∞ x ( ) ∫ l 2 c G ( k ) ∂ α π c ∂ α ac mc C = c N + N cos Λ cdx (176) L L u u ¨ u α S 2 V k ∂ U 2 V ∂ U ∞ t ∞ t ∫ l 2 ∂ α ac C = c N F ( k ) cos Λ cdx (177) L L θ α θ S ∂ Θ ( ) ∫ l 2 ∂ α c G ( k ) ∂ α π c ∂ α ac ac mc C = c N + c N + N cos Λ cdx (178) L L θ L θ θ ˙ α α θ S ∂ Θ 2 V k ∂ Θ 2 V ∂ Θ t ∞ ∞ ( ) ∫ l 2 c G ( k ) ∂ α π c ∂ α ac mc C = c N + N cos Λ cdx (179) L L θ θ ¨ α θ S 2 V k ∂ Θ 2 V ∂ Θ ∞ t ∞ t ∫ L C = c cos Λ cdx (180) L L δ δ S − L Assuming a parabolic drag polar, the aircraft unsteady drag coefficient is evaluated from the unsteady lift coeffi- cient as C ( k ) L 2 2 C ( k ) = C + + C δ + C δ + C δ + C δ (181) D D D e D D r D 0 2 e 2 r δ δ e δ r δ π AR ε e r where C , C , C , and C are the drag derivatives due to the elevator and rudder deflections.

D D D D δ 2 δ 2 e r δ δ e r Neglecting the drag contribution, the wing contribution to the aircraft unsteady pitching moment coefficient is evaluated as [ ] ∫ L ( ) 1 ˙ α c G ( k ) π ˙ α c ac mc ∆ C ( k ) = cc − x c α F ( k ) − x c − x + c − x c δ cos Λ cdx (182) m m ac L ac ax L mc m ac L ac α α δ δ S ¯ c 2 V k 2 V − L ∞ ∞ This is expressed as ( ) n ˙ ¨ ∆ C = C + C s + C ˙ s + C u + C ˙ u + C ¨ u + C θ + C θ + C θ + C δ (183) m m m m m i m i m i m i m i m i m

0 s ˙ s ∑ u ˙ u ¨ u θ ˙ ¨

δ θ θ i = 1 Due to symmetry, the partial derivative contributions to the aircraft unsteady lift and pitching moment do not involve aircraft lateral-directional states β , p , and r since these terms cancel out in the integration over both the left and right wings.

The aircraft unsteady rolling moment coefficient is evaluated as ( ) ∫ L 1 ˙ α c G ( k ) π ˙ α c ac mc C ( k ) = y c α F ( k ) + y c + y + y c δ cos Λ cdx (184) ac L ac ac L mc ac L l α α δ Sb 2 V k 2 V − L ∞ ∞ The aircraft unsteady yawing moment coefficient is evaluated as [ ( )] ∫ L 1 c L 2 C ( k ) = − y c + cos Λ cdx + C δ (185) n ac D n r 0 δ r Sb π AR ε − L Due to symmetry, the partial derivative contributions to the aircraft unsteady rolling and yawing moments do not involve aircraft longitudinal states α and q since these terms cancel out in the integration over both the left and right wings.

VI. Vortex-Lattice Aerodynamic Model Coupling Vorview is a computational aerodynamic tool that is used for the development of the aeroelastic computational capability. Vorview provides a rapid method for estimating aerodynamic force and moment coefficients as well as aerodynamic stability and control derivatives for a given aircraft configuration. It is based on the vortex-lattice lifting line aerodynamic theory. The vehicle configuration is constructed within Vorview by a series of panels that are formed by spanwise and chordwise locations of bound vortices. Vorview computes the vehicle aerodynamics in both the longitudinal and lateral directions independently. The longitudinal and lateral aerodynamics are then combined to produce overall aerodynamic characteristics of the vehicle at any arbitrary angle of attack and angle of sideslip. Due to the inviscid nature of any vortex-lattice method, the drag prediction by Vorview is most reliable 24 of 38 American Institute of Aeronautics and Astronautics for induced drag prediction. For viscous drag due to boundary layer separation or wave drag due to shock-induced boundary layer separation, the prediction may be less reliable. Vorview can provide a rapid estimation of aerodynamic derivatives including dynamic derivatives due to angular rates. Owing to the computationally efficient vortex-lattice method, aerodynamic derivatives can be estimated in Vorview fairly quickly. A flight dynamic model for a given vehicle configuration can be easily developed with Vorview that supplies the model with all necessary aerodynamic 7 16 information for the vehicle. Vorview has been validated by both wind tunnel data as well as the NASA Cart3D tool, which is a high-fidelity inviscid (Euler) CFD analysis code targeted at analyzing aircraft performance in conceptual and preliminary aerodynamic design. In general, both Vorview and Cart3D seem to have similar predictive capabilities when compressibility is not a factor.

Figure 11 illustrates an aerodynamic model of the GTM in Vorview.

Fig. 11 - Vorview Aircraft Model A. Automated Vehicle Geometry Modeling Tool To enable a coupled aeroelastic solution, the aircraft deformed geometry must be generated at each iteration. An automated vehicle geometry modeling tool has been developed in MATLAB to update the aircraft deformed geometry.

The vehicle geometry modeler directly outputs a geometry input file that can be read by Vorview during a solution cycle.

Fig. 12 - GTM Coordinate Systems With reference to Fig. 12, the coordinate reference frame ( x , y , z ) defines the Body Station (BS), the Body Butt B B B Line (BBL), and the Body Water Line (BWL) of the aircraft, respectively. The coordinate reference frame ( x , y , z ) V V V 25 of 38 American Institute of Aeronautics and Astronautics is the translated coordinate system attached to the nose of the aircraft. The stability reference frame ( x , y , z ) is attached to the CG such that x = ¯ x − x , y = y − ¯ y , and z = ¯ z − z , where ( ¯ x , ¯ y , ¯ z ) is the coordinate of the CG in the V V V V V V V V V ( x , y , z ) reference frame.

V V V The vehicle geometry modeler has access to the outer mold line of the jig-shape (undeformed) aircraft geometry.

The coordinate reference frame ( x , y , z ) defines the coordinate system used in the vehicle geometry model.

V V V The aeroelastic deflections in bending and torsion result in a displacement ∆ p and rotation angle φ where φ = Θ d − W d + V d (186) 1 x 2 x 3 ∆ r = − ( W sin W + V sin V ) d + V cos V d + W cos W d (187) x x 1 x 2 x 3 The coordinate reference frame ( x , y , z ) is related to the coordinate reference frame ( x , y , z ) by the following V V V relationship:       d − sin Λ cos Γ − cos Λ cos Γ − sin Γ b 1 1       =  d   − cos Λ sin Λ 0   b  2 2 d sin Λ sin Γ cos Λ sin Γ − cos Γ b 3 3     − sin Λ cos Γ − cos Λ cos Γ − sin Γ − v     = (188)  − cos Λ sin Λ 0   v  sin Λ sin Γ cos Λ sin Γ − cos Γ − v where ( v , v , v ) are the unit vectors for the Vorview coordinate reference frame ( x , y , z ) .

1 2 3 V V V Thus, the aeroelastic deflections result in a wing twist expressed as an incremental rotation vector ( ∆ φ , ∆ φ , ∆ φ ) x y z and a displacement vector ( ∆ x , ∆ y , ∆ z ) V V V ∆ φ = Θ sin Λ cos Γ − W cos Λ − V sin Λ sin Γ (189) x x x ∆ φ = − Θ cos Λ cos Γ − W sin Λ + V cos Λ sin Γ (190) y x x ∆ φ = Θ sin Γ + V cos Γ (191) z x ∆ x = − ( W sin W + V sin V ) sin Λ cos Γ + V cos V cos Λ − W cos W sin Λ sin Γ (192) V x x x x ∆ y = ( W sin W + V sin V ) cos Λ cos Γ + V cos V sin Λ + W cos W cos Λ sin Γ (193) V x x x x ∆ z = − ( W sin W + V sin V ) sin Γ + W cos W cos Γ (194) V x x x A coordinate transformation to account for wing aeroelastic deflections is performed by rotating a wing section about its elastic axis by the incremental rotation vector ( ∆ φ , ∆ φ , ∆ φ ) and then translating the resultant coordinates x y z by the displacement vector ( ∆ x , ∆ y , ∆ z ) .

V V V To perform the coupled aeroelastic computation, the static aeroelastic model is coupled with Vorview through the automated vehicle geometry modeler. Aerodynamic force and moment coefficients as computed from Vorview are used as inputs to the static aeroelastic model. The computed aeroelastic deflections are then used to generate the aircraft deformed geometry by the automated vehicle geometry modeler. The aerodynamic solution is then recomputed with the aircraft deformed geometry in Vorview. This process is iterated until the solution is converged when errors in the computed aeroelastic deflections are within a specified tolerance. A flow chart for the coupled aeroelastic computation is shown in Fig. 13.

The static aeroelastic solution provides aerodynamic information for the deformed aircraft under a trimmed flight condition. The dynamic aeroelastic analysis is conducted to compute the unsteady contributions to wing aerodynamics.

26 of 38 American Institute of Aeronautics and Astronautics Fig. 13- Coupled Aeroelastic Vortex Lattice Computation Flow Chart VII. Simulations A coupled aeroelastic-longitudinal flight dynamic model is built by coupling the wing dynamic aeroelastic model to a linearized model of the aircraft rigid-body longitudinal flight dynamics. The wing dynamic aeroelastic model is represented by a second-order system described by Eq. (134) which is assembled using the finite-element method.

In this implementation, the wing dynamic aeroelastic model is implemented using the cantilever boundary conditions where W ( 0 ) , W ( 0 ) , V ( 0 ) , V ( 0 ) , and Θ ( 0 ) are all set to zero. Strictly speaking, the wing symmetric modes should be x x coupled with the aircraft longitudinal flight dynamics where the wing displacement at the wing root matches the aircraft displacement. By removing the wing displacement at the wing root through the cantilever boundary conditions, the coupled aeroelastic-longitudinal flight dynamic model is a reasonable approximation of the coupled wing symmetric modes with the aircraft rigid-body modes.

Let the aircraft rigid-body flight dynamics be represented by a first-order system given by M ˙ x = Sx + C u (195) r r r u [ ] > where x is the aircraft rigid-body state vector, u = is a control vector of the control surface deflections δ δ r e [ ] comprising the elevator δ and the symmetric VCCTEF deflection δ , and C = is the partitioned control C C e u δ δ e sensitivity matrix.

Note that the design concept of the VCCTEF utilizes only the third camber segments of each spanwise flap section for roll and pitch control due to the faster response of these control surfaces provided by the electric drive motors.

[ ] > Thus, δ = represents the deflections of the third camber segments of the 15 spanwise δ δ δ . . . δ 1 2 3 flap sections of the VCCTEF.

Due to the coupling between the aircraft unsteady aerodynamic coefficients with the aeroelastic deformation of the wing, the coupling matrices represented by H ( k ) , H ( k ) , and H ( k ) are introduced, where k is the reduced frequency 1 2 3 parameter. These matrices are calculated using the unsteady aircraft aerodynamic coefficients evaluated using the equations developed in the previous section. Thus, the first-order system representing the aircraft rigid-body dynamics coupled to the wing aeroelastic states is given by M ˙ x = Sx + H ( k ) x + H ( k ) ˙ x + H ( k ) ¨ x + C u (196) r r r e e e u 1 2 3 27 of 38 American Institute of Aeronautics and Astronautics The coupled equations (196) and (134) can be combined together to form a coupled first-order state space model         − 1 M 0 − H ( k ) S H ( k ) H ( k ) ˙ x x r r 3 1 2 r          ˙ x  =  0 I 0   0 0 I   x  e e ∂ F ( k ) ∂ F ( k ) ¨ x − 0 M ( k ) − K ( k ) − C ( k ) ˙ x e e ∂ ˙ s ∂ s ︸ ︷︷ ︸ ︸ ︷︷ ︸ ˙ x A     − 1 [ ] M 0 − H ( k ) C C r 3 δ δ e δ     e +  0 I 0   0 0  (197) δ ∂ F ( k ) ∂ F − 0 M ( k ) 0 ︸ ︷︷ ︸ ∂ ˙ s ∂ δ ︸ ︷︷ ︸ u B Note that the coupled state space model represents a reduced-frequency-dependent state space model. This model is generally valid for a known value of the reduced frequency parameter k and can be used to approximate a flight dy- namic model if the dominant wing aeroelastic frequency is known. Another method is to remove the reduced frequency dependency from the state space model by approximating the Theodorsen’s function C ( k ) by various methods such as 17 18 the Roger method of rational fraction approximation or the R.T. Jones method. These approximation methods can be advantageous in that the state space model is more broadly applicable to a wider range of frequencies at a cost of increasing the order of the state space model by introducing aerodynamic lag states resulting from the approximation of the Theodorsen’s function. This approach will be considered in the future work.

A. Linearized Aircraft Rigid-Body Longitudinal Flight Dynamic Model A 4-state longitudinal flight dynamic model of the GTM is implemented. The aircraft rigid-body states are given by [ ] > ∆ V ∆ V x = ∆ α q ∆ θ where is a the normalized perturbation of the aircraft forward airspeed, ∆ α is the r V V perturbation of the aircraft angle of attack, q is the pitch rate, and ∆ θ is the perturbation in the aircraft pitch angle. For the flight condition of Mach 0.8 at 35,000 ft, the matrices for the aircraft rigid-body longitudinal flight dynamic model are given by   11 . 1138 0 0 0   0 11 . 1757 0 0   M =   r   0 0 . 1310 0 . 7841 0 0 0 0 1   − 0 . 0558 − 0 . 4364 − 0 . 7480 − 0 . 4595   − 1 . 7284 − 6 . 3068 10 . 9544 − 0 . 0306   S =     − 0 . 0074 − 1 . 7648 − 0 . 3370 0 0 0 1 0 o o which represent a linearization about the rigid-body trim point α = 3 . 8142 , V = 778 . 2063 ft/sec, θ = − 3 . 8142 , o T = 5617 lbs, and δ = − 6 . 1497 .

e The eigenvalues of the aircraft rigid-body longitudinal flight dynamics are calculated to be λ = − 0 . 5779 ± 1 . 4491 i sp λ = − 0 . 0042 ± 0 . 0763 i p The eigenvalues λ correspond to the high frequency, highly damped short-period mode of the aircraft dynamics.

sp The eigenvalues λ represent the low frequency, lightly damped phugoid mode of the aircraft dynamics involving the p airspeed, pitch angle, and altitude.

The natural frequencies and the damping ratios for the short-period and the phugoid modes are calculated to be ω = 1 . 5601 rad/sec n sp ζ = 0 . 3704 sp ω = 0 . 0764 rad/sec n p ζ = 0 . 0545 p 28 of 38 American Institute of Aeronautics and Astronautics B. Flutter Analysis Two different wing models will be analyzed for the baseline wing stiffness of the GTM wings and for the reduced stiffness of the highly flexible wings of the ESAC. The baseline structural rigidities EI and GJ are estimated for the conventional stiff wing structures for the GTM. For the ESAC, the wing structural rigidities EI and GJ are purposely reduced by a factor of two to model highly flexible wing structures. The increased flexibility enables the wing shaping control actuation by the VCCTEF system for drag reduction.

In addition to the wing dry mass, the fuel mass is also accounted for. The fuel is stored in the center tank and wing main tanks. The center tank holds 20,000 lbs of fuel. Each of the main tanks holds about 15,000 lbs of fuel.

The center tank is used first until it is empty. Then, the fuel is drawn equally from the wing main tanks. The fuel mass is modeled as the combined wing mass density. As the structural rigidities are reduced, the wing dry mass also decreases. Assuming that the wing box structure is modeled as a thin-walled structure, then the mass change is related to the change in the wing structural rigidity EI can be modeled.

For the flutter analysis, the structural dynamic modes for the cantilever, symmetric, and anti-symmetric boundary conditions are computed with 80% fuel loading and no structural damping for conservatism. A trim thrust value is used in the flutter analysis based on the linearization of the aircraft rigid-body longitudinal flight dynamic model to account for aero-propulsive-elastic effects. The structural dynamic cantilever mode shapes of the stiff GTM wings are shown in Fig. 14 and their associated natural frequencies for the stiff GTM wings and flexible ESAC wings are summarized in Table 1.

Mode 1B Mode 2B Mode 1T/2B Mode 3B Mode 2T/3B Mode 2T/4B Fig. 14 - Structural Dynamic Cantilever Mode Shapes of Stiff GTM Wings Mode Natural Frequency, Hz (GTM Wings) Natural Frequency, Hz (ESAC Wings) 1B 1.4934 1.1252 2B 4.0522 2.9620 1T/2B 5.0930 3.7895 3B 9.5070 7.1341 2T/3B 10.9926 8.4271 2T/4B 17.4842 13.3635 Table 1 - Structural Dynamic Natural Frequencies of Cantilever Modes with 80% Fuel Loading 29 of 38 American Institute of Aeronautics and Astronautics The flutter analysis using a linear aeroelastic model is conducted for the cantilevered wing with and without coupling to the aircraft rigid-body longitudinal flight dynamic model. The flutter analysis of the coupled aeroelastic- longitudinal flight dynamic model allows for solving for the flutter frequency that will be used to calculate the reduced frequency parameter k needed to generate the coupled state space model. A cruise condition at 35,000 ft altitude is examined.

1. Stiff GTM Wings The frequency and damping ratios of the stiff GTM wings with the baseline stiffness were computed by sweeping over a Mach number range. Figure 15 is a plot of the aeroelastic frequencies and damping ratios for the uncoupled wing cantilever modes, while Table 2 summarizes the critical flutter mach numbers and flutter frequencies for the first two flutter modes. The critical flutter mode is observed to be due to the second bending mode (2B) at Mach 1.3792 with a flutter frequency of 2.4792 Hz.

Stiff Wing Uncoupled Cantilever Modes @ 35,000 ft w/ 80% Fuel Loading Stiff Wing Uncoupled Cantilever Modes @ 35,000 ft w/ 80% Fuel Loading 11 0.4 10 0.35 0.3 8 0.25 7 0.2 0.15 ζ , Hz ω 5 0.1 4 0.05 2 −0.05 −0.1 0.7 0.8 0.9 1 1.1 1.2 1.3 1.4 1.5 1.6 0.7 0.8 0.9 1 1.1 1.2 1.3 1.4 1.5 1.6 M M ∞ ∞ Fig. 15 - Flutter Results for Uncoupled Cantilever Modes of Stiff GTM Wings at 35,000 ft with 80% Fuel Loading Stiff Wing Coupled Cantilever Modes @ 35,000 ft w/ 80% Fuel Loading Stiff Wing Coupled Cantilever Modes @ 35,000 ft w/ 80% Fuel Loading 11 0.3 0.25 0.2 0.15 ζ , Hz ω 0.1 0.05 1 −0.05 0.7 0.8 0.9 1 1.1 1.2 1.3 1.4 1.5 1.6 0.7 0.8 0.9 1 1.1 1.2 1.3 1.4 1.5 1.6 M M ∞ ∞ Fig. 16 - Flutter Results for Coupled Cantilever Modes of Stiff GTM Wings at 35,000 ft with 80% Fuel Loading Mode Flutter Mach Flutter Frequency, Hz 2B 1.3792 2.4792 3B 1.6729 6.4578 Table 2 - First Two Flutter Speeds of Uncoupled Cantilever Modes of Stiff GTM Wings at 35,000 ft with 80% Fuel Loading 30 of 38 American Institute of Aeronautics and Astronautics The flutter analysis is repeated for the coupled aeroelastic-longitudinal flight dynamic model. The critical flutter mode is observed to be due to the third bending mode (3B) at Mach 1.5490 with a flutter frequency of 9.6130 Hz.

Note that the flutter characteristics are significantly changed as a result of the coupling with the aircraft rigid-body longitudinal flight dynamics.

2. Flexible ESAC Wings The stiffness of the flexible ESAC wings is reduced by half from the baseline stiffness of the GTM wings. The flexibility will allow the VCCTEF to be more effective in wing shaping control for drag reduction. As a result of the reduced stiffness, the flutter boundary of the flexible ESAC wing will decrease from that of the stiff GTM wings.

Without coupling to the aircraft rigid-body longitudinal flight dynamics, the flutter characteristics of the flexible ESAC wings are shown in Fig. 17 and Table 3. The critical flutter mode is observed to be due to the second bending mode (2B) at Mach 1.0115 with a flutter frequency of 3.2367 Hz.

Flexible Wing Uncoupled Cantilever Modes @ 35,000 ft w/ 80% Fuel Loading Flexible Wing Uncoupled Cantilever Modes @ 35,000 ft w/ 80% Fuel Loading 9 0.6 8 0.5 7 0.4 6 0.3 5 0.2 ζ , Hz ω 4 0.1 3 0 2 −0.1 1 −0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 1.1 1.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 1.1 1.2 M M ∞ ∞ Fig. 17 - Flutter Results for Uncoupled Cantilever Modes of Flexible ESAC Wings at 35,000 ft with 80% Fuel Loading Flexible Wing Coupled Cantilever Modes @ 35,000 ft w/ 80% Fuel Loading Flexible Wing Coupled Cantilever Modes @ 35,000 ft w/ 80% Fuel Loading 9 0.4 0.35 0.3 0.25 0.2 ζ 5 0.15 , Hz ω 0.1 0.05 −0.05 1 −0.1 0.4 0.6 0.8 1 1.2 1.4 0.4 0.6 0.8 1 1.2 1.4 M M ∞ ∞ Fig. 18 - Flutter Results for Coupled Cantilever Modes of Flexible ESAC Wings at 35,000 ft with 80% Fuel Loading Mode Flutter Mach Flutter Frequency, Hz 2B 1.0115 3.2367 1T/2B 1.1888 5.0131 Table 3 - First Two Flutter Speeds Uncoupled Cantilever Modes of Flexible ESAC Wings at 35,000 ft with 80% Fuel Loading 31 of 38 American Institute of Aeronautics and Astronautics The flutter characteristics of the coupled cantilever modes of the flexible wings are shown in Fig. 18. The critical flutter mode is observed to be due to the 2T/4B mode which occurs at Mach 1.3105 with a flutter frequency of 6.2244 Hz.

C. Modal Analysis The flutter analysis can be used to establish the most dominant aeroelastic mode at a given flight condition. This can be determined by examining the aeroelastic mode with the lowest damping. In all cases, the second bending mode is the most lightly damped mode. The frequency of this mode is then used to determine the reduced frequency parameter k for the coupled aeroelastic-longitudinal flight dynamic model.

1. Stiff GTM Wings A modal analysis of the coupled aeroelastic-longitudinal flight dynamics is performed to determine the effect the aeroelastic coupling on the aircraft rigid-body modes. The eigenvalues of the aircraft rigid-body longitudinal flight dynamics with the stiff GTM wings are calculated to be λ = − 0 . 6012 ± 1 . 3949 i sp λ = − 0 . 0051 ± 0 . 0771 i p The natural frequencies and the damping ratios for the short-period and the phugoid modes are calculated to be ω = 1 . 5189 rad/sec n sp ζ = 0 . 3958 sp ω = 0 . 0773 rad/sec n p ζ = 0 . 0655 p The frequencies and damping ratios for the coupled cantilever modes of the stiff GTM wings are shown in Table 4.

The root locus of the poles for the phugoid, short-period, and the first three cantilever modes of the stiff GTM wings are plotted in Fig. 19.

Mode Eigenvalue Frequency (rad/s) Damping Factor 1B − 0 . 7755 ± 10 . 4433 i 1.6667 0.0741 2B − 0 . 3637 ± 25 . 8867 i 4.1204 0.0140 1T/2B − 0 . 6130 ± 32 . 2080 i 5.1270 0.0190 3B − 0 . 7252 ± 59 . 5928 i 9.4852 0.0122 2T/3B − 1 . 2371 ± 66 . 9549 i 10.6580 0.0185 2T/4B − 1 . 5012 ± 108 . 3448 i 17.2453 0.0139 Table 4 - Natural Frequencies and Damping Ratios of Coupled Cantilever Modes of Stiff GTM Wings at 35,000 ft with 80% Fuel Loading 32 of 38 American Institute of Aeronautics and Astronautics Root Locus of Stiff Wing Cantilever Modes @ 35,000 ft w/ 80% Fuel Loading Uncoupled Phugoid Uncoupled Short−Period Uncoupled 1B Uncoupled 2B Uncoupled 1T/2B ω 0 i Coupled Phugoid Coupled Short−Period Coupled 1B −10 Coupled 2B Coupled 1T/2B −20 −30 −40 −0.8 −0.7 −0.6 −0.5 −0.4 −0.3 −0.2 −0.1 0 σ Fig. 19 - Root Locus of Coupled Aeroelastic-Longitudinal Flight Dynamics with Stiff GTM Wings at 35,000 ft with 80% Fuel Loading 2. Flexible ESAC Wings The eigenvalues of the aircraft rigid-body longitudinal flight dynamics with the flexible ESAC wings are calculated to be λ = − 0 . 6570 ± 1 . 2712 i sp λ = − 0 . 0063 ± 0 . 0782 i (198) p The natural frequencies and the damping ratios for the short-period and the phugoid modes are calculated to be ω = 1 . 4310 rad/sec n sp ζ = 0 . 4591 sp ω = 0 . 0785 rad/sec n p ζ = 0 . 0807 p Note that both the rigid-body modes are stable in spite of the aeroelastic coupling. This is due to the significant frequency separation between the rigid-body modes and the wing cantilever modes.

The frequencies and damping ratios for the coupled cantilever modes of the stiff GTM wings are shown in Table 5.

The root locus of poles for the phugoid, short-period, and the first three cantilever modes of the flexible ESAC wings are plotted in Fig. 20.

Mode Eigenvalue Frequency (rad/s) Damping Factor 1B − 1 . 0929 ± 9 . 3646 i 9.8499 0.1159 2B − 0 . 2881 ± 19 . 6302 i 19.6323 0.0147 1T/2B − 1 . 0147 ± 25 . 3135 i 25.3339 0.0401 3B − 0 . 8460 ± 46 . 4797 i 46.4874 0.0182 2T/3B − 5 . 1495 ± 46 . 6200 i 46.9035 0.1098 2T/4B − 6 . 7807 − 81 . 6900 i 81.9709 0.0827 Table 4 - Natural Frequencies and Damping Ratios of Coupled Cantilever Modes of Flexible ESAC Wings at 35,000 ft with 80% Fuel Loading 33 of 38 American Institute of Aeronautics and Astronautics Root Locus of Flexible Wing Cantilever Modes @ 35,000 ft w/ 80% Fuel Loading Uncoupled Phugoid 10 Uncoupled Short−Period Uncoupled 1B Uncoupled 2B Uncoupled 1T/2B ω 0 i Coupled Phugoid Coupled Short−Period Coupled 1B −10 Coupled 2B Coupled 1T/2B −20 −30 −1.4 −1.2 −1 −0.8 −0.6 −0.4 −0.2 0 σ Fig. 20 - Root Locus of Coupled Aeroelastic-Longitudinal Flight Dynamics with Flexible ESAC Wings at 35,000 ft with 80% Fuel Loading D. Dynamic Response Simulations 1. Elevator Input The dynamic response for aircraft with the stiff GTM wings with the baseline stiffness is simulated for separate inputs of the elevator and VCCTEF over a time span of 40 sec. Figure 21 shows a doublet input of the elevator.

1.5 0.5 , deg e δ −0.5 −1 −1.5 0 5 10 15 20 25 30 35 40 t, sec Fig. 21 - Elevator Doublet The dynamic responses of the aircraft with rigid-body and coupled aeroelastic-longitudinal flight dynamics with the stiff GTM wings and flexible ESAC wings are shown in Fig. 22. The dynamic response due to the stiff GTM wings matches that of the rigid-body longitudinal flight dynamics very well. The dynamic response due to the flexible ESAC wings also matches the rigid-body dynamic response reasonably well, but some small differences are noted.

34 of 38 American Institute of Aeronautics and Astronautics 20 2 Rigid Rigid Stiff Wing Stiff Wing 1.5 15 Flexible Wing Flexible Wing 0.5 V, ft , deg α −0.5 −5 −1 −10 −1.5 −15 −2 0 5 10 15 20 25 30 35 40 0 5 10 15 20 25 30 35 40 t, sec t, sec 3 3 Rigid Stiff Wing 2.5 Flexible Wing 1.5 0.5 −1 , deg θ q, deg/sec −2 −0.5 −3 −1 Rigid −4 −1.5 Stiff Wing Flexible Wing −5 −2 0 5 10 15 20 25 30 35 40 0 5 10 15 20 25 30 35 40 t, sec t, sec Fig. 22 - Aircraft Longitudinal States The dynamic responses of the aeroelastic deflections at the wing tip, denoted as W and Θ are shown in Fig.

tip tip 23. The flexible ESAC wings experience much greater aeroelastic deflections as expected since the stiffness of the ESAC wings is half that of the stiff GTM wings. It is somewhat surprising that, with the significant wing aeroelastic deflections, the aircraft longitudinal states do not seem to be much affected.

2 1.5 Rigid Rigid Stiff Wing Stiff Wing 1.5 1 Flexible Wing Flexible Wing 1 0.5 0.5 0 , ft , deg tip tip W Θ 0 −0.5 −0.5 −1 −1 −1.5 −1.5 −2 0 5 10 15 20 25 30 35 40 0 5 10 15 20 25 30 35 40 t, sec t, sec Fig. 23 - Aeroelastic Deflections at Wing Tip 35 of 38 American Institute of Aeronautics and Astronautics 2. VCCTEF Input The VCCTEF can be used for either roll or pitch control by symmetric or anti-symmetric deflections of the individual spanwise flap segments. The first inboard flap segment is designated as a high-lift flap and therefore is not used for roll or pitch control. Because of the continuous trailing edge surfaces, the flap segment adjacent to the inboard flap can only be deflected by a relative amount as permitted by the transition material.

For the dynamic response simulation of the VCCTEF, the outboard flap, designated flap 15, is commanded by the same doublet as the elevator. The commands for flaps 1 to 15 vary linearly from zero to the full doublet.

The dynamic responses of the aircraft with rigid-body and coupled aeroelastic-longitudinal flight dynamics with the stiff GTM wings and flexible ESAC wings are shown in Fig. 24. These dynamic responses are significantly different from one another. In particular, the dynamic response due to the flexible ESAC wings exhibits control reversals for α , q , and θ . The control reversals are due to a large nose-down twist caused by the flexible ESAC wings.

3 0.4 Rigid Stiff Wing 0.3 Flexible Wing 0.2 0.1 −1 , deg α V, ft/sec −2 −0.1 −3 Rigid −0.2 −4 Stiff Wing Flexible Wing −5 −0.3 0 5 10 15 20 25 30 35 40 0 5 10 15 20 25 30 35 40 t, sec t, sec 0.5 0.6 Rigid Stiff Wing 0.4 Flexible Wing 0.4 0.3 0.2 0.2 0.1 , deg θ −0.2 q, deg/sec −0.4 −0.1 −0.6 −0.2 −0.3 −0.8 0 5 10 15 20 25 30 35 40 0 5 10 15 20 25 30 35 40 t, sec t, sec Fig. 24 - Aircraft Longitudinal States The dynamic responses of the aeroelastic deflections at the wing tip, denoted as W and Θ are shown in Fig.

tip tip 25. The flexible ESAC wings experience much greater aeroelastic deflections and dynamic transients than the stiff GTM wings. In particular, the twist of the flexible ESAC wings is quite substantial, more than 2 degrees at some time instances. This could explain the sign reversal in α , q , and θ .

36 of 38 American Institute of Aeronautics and Astronautics 1 3 Rigid Rigid Stiff Wing Stiff Wing 0.8 Flexible Wing 2 Flexible Wing 0.6 0.4 0.2 , ft , deg tip tip W 0 Θ −1 −0.2 −2 −0.4 −3 −0.6 −0.8 −4 0 5 10 15 20 25 30 35 40 0 5 10 15 20 25 30 35 40 t, sec t, sec Fig. 25 - Aeroelastic Deflections at Wing Tip VIII. Conclusions This paper presents a coupled aeroelastic-longitudinal flight dynamic model of a flexible wing aircraft. This aircraft concept called Elastically Shaped Aircraft Concept (ESAC). The aircraft concept addresses the drag reduction goal in commercial aviation through an elastic wing shaping control approach for aircraft with highly flexible wing structures.

The multi-disciplinary nature of flight physics is appreciated with the recognition of the potential adverse effects of aeroelastic wing shape deflections on aerodynamic performance. By aeroelastically tailoring the wing shape with active control, a significant drag reduction benefit could be realized. To attain the potential of the elastic wing shaping control concept, a new type of aerodynamic control effector is introduced and is referred to as a Variable Camber Continuous Trailing Edge Flap (VCCTEF).

A coupled aeroelastic-flight dynamic modeling approach has been developed for this aircraft concept. The flight dynamic model is coupled with the aeroelastic states from a finite-element model of the flexible wing via the aeroe- lastic contributions to the aerodynamic coefficients. A coupled aeroelastic-longitudinal flight dynamic model has been developed for both the stiff wing GTM and flexible wing ESAC. Initial simulations show that the short-period and phugoid modes remain stable for the flexible wing ESAC. An open-loop response simulation is conducted to demonstrate the coupled dynamic response. The wing flexibility results in a significant deflection as compared to a conventional stiff wing transport aircraft which could cause some issues of control reversal. Future work will further develop the coupled aeroelastic lateral-directional dynamic model and ultimately a fully coupled 6-dof flight dynamic model.

References Nguyen, N., “Elastically Shaped Future Air Vehicle Concept,” NASA Innovation Fund Award 2010 Report, October 2010, Submitted to NASA Innovative Partnerships Program, http://ntrs.nasa.gov/archive/nasa/casi.ntrs.nasa.gov/20110023698_2011024909.pdf Nguyen, N., Trinh, K., Reynolds, K., Kless, J., Aftosmis, M., Urnes, J., and Ippolito, C., “Elastically Shaped Wing Optimization and Aircraft Concept for Improved Cruise Efficiency,” AIAA Aerospace Sciences Meeting, AIAA-2013-0141, January 2013.

Boeing Report No. 2012X0015, “Development of Variable Camber Continuous Trailing Edge Flap System,” October 4, 2012.

Urnes, J., Nguyen, N., Ippolito, C., Totah, J., Trinh, K., and Ting, E., “A Mission Adaptive Variable Camber Flap Control System to Optimize High Lift and Cruise Lift to Drag Ratios of Future N+3 Transport Aircraft,” AIAA Aerospace Sciences Meeting, AIAA-2013-0214, January 2013.

Jordan, T. L., Langford, W. M., Belcastro, C. M., Foster, J. M., Shah, G. H., Howland, G., and Kidd, R., “Development of a Dynamically Scaled Generic Transport Model Testbed for Flight Research Experiments,” AUVSI Unmanned Unlimited, Arlington, VA, 2004.

Nguyen, N. and Urnes, J., “Aeroelastic Modeling of Elastically Shaped Aircraft Concept via Wing Shaping Control for Drag Reduction," AIAA Atmospheric Flight Mechanics Conference, AIAA-2012-4642, August 2012.

Nguyen, N., Nelson, A., and Pulliam, T., “Damage Adaptive Control System Research Report,” Internal NASA Report, April 2006.

Stanford University, AA241, http://adg.stanford.edu/aa241/drag/lsformfactor.html Hodges, D.H. and Pierce, G.A., Introduction to Structural Dynamics and Aeroelasticity , Cambridge University Press, 2002.

Nguyen, N., “Integrated Flight Dynamics Modeling of Flexible Aircraft with Inertial Force-Propulsion - Aeroelastic Couplings,” 46th AIAA Aerospace Sciences Meeting and Exhibit, AIAA-2008-194, January 2008.

37 of 38 American Institute of Aeronautics and Astronautics Houbolt, J. C. and Brooks, G. W., “Differential Equations of Motion for Combined Flapwise Bending, Chordwise Bending, and Torsion of Twisted Nonuniform Rotor Blades,” NACA Technical Note 3905, February 1957.

Nguyen, N., Trinh, K., Nguyen, D., and Tuzcu, I., “Nonlinear Aeroelasticity of Flexible Wing Structure Coupled with Aircraft Flight Dynamics,” AIAA Structures, Structural Dynamics, and Materials Conference, AIAA-2012-1792, April 2012.

Theodorsen, T. and Garrick, I.E., “Mechanism of Flutter - a Theoretical and Experimental Investigation of the Flutter Problem”, NACA Report 685, 1940.Theodorsen, T, “General Theory of Aerodynamic Instability and the Mechanism of Flutter”, NACA Report 496, 1935.

Hughes, T., The Finite Element Method Linear Static and Dynamic Finite Element Analysis , Prentice-Hall, Inc, 1987.

Miranda, L.R., Elliot, R.D., and Baker, W.M., “A Generalized Vortex Lattice Method for Subsonic and Supersonic Flow Applica- tions,” NASA CR-2865, 1977.

Aftosmis, M. J., Berger, M. J., and Melton, J. E., “Robust and Efficient CartesianMesh Generation for Component- Based Geometry,” AIAA Journal, Vol. 36, No. 6, 1998, pp. 953-960.

Abel, I., “An Analytical Technique for Predicting the Characteristics of a Flexible Wing Equipped with an Active Flutter-Suppression System and Comparison with Wind-Tunnel Data,” NASA TP 1367, 1979.

Brunton, S. and Rowley, C., “Empirical State-Space Representations for Theodorsen’s Lift Model,” Elsevier Journal of Fluids and Structures, Vol. 38, April 2013, pp. 174–186.

38 of 38 American Institute of Aeronautics and Astronautics