Static Aeroelastic Scaling and Analysis of a Sub-Scale Flexible Wing Wind Tunnel Model

Static Aeroelastic Scaling and Analysis of a Sub-Scale Flexible

Wing Wind Tunnel Model

∗

Eric Ting

Stinger Ghaffarian Technologies, Inc., Moffett Field, CA 94035

†

Sonia Lebofsky

Stinger Ghaffarian Technologies, Inc., Moffett Field, CA 94035

‡

Nhan Nguyen

NASA Ames Research Center, Moffett Field, CA 94035

§

Khanh Trinh

Stinger Ghaffarian Technologies, Inc., Moffett Field, CA 94035

This paper presents an approach to the development of a scaled wind tunnel model for static aeroelastic similarity with a full-scale wing model. The full-scale aircraft model is based on the NASA Generic Transport Model (GTM) with flexible wing structures referred to as the Elastically Shaped Aircraft Concept (ESAC).

The baseline stiffness of the ESAC wing represents a conventionally stiff wing model. Static aeroelastic scaling is conducted on the stiff wing configuration to develop the wind tunnel model, but additional tailoring is also conducted such that the wind tunnel model achieves a 10% wing tip deflection at the wind tunnel test condition.

An aeroelastic scaling procedure and analysis is conducted, and a sub-scale flexible wind tunnel model based on the full-scale’s undeformed jig-shape is developed. Optimization of the flexible wind tunnel model’s undeflected twist along the span, or pre-twist or wash-out, is then conducted for the design test condition. The resulting wind tunnel model is an aeroelastic model designed for the wind tunnel test condition.

I. Introduction

Due to recent strides in the development of light-weight materials, the aircraft industry has been investigating the possibility of reducing airframe weight to increase energy efficiency. Reduction of the aircraft weight translates into a lower lift requirement and in turn, reduces induced drag, thrust requirements, fuel burn, and cost. These modern materials, such as advanced composites, are able to maintain the same load-carrying capacity of conventional air- frame material selections. The provided structural rigidity of these materials, however, can be reduced. It becomes increasingly important for these modern designs to take into account the aeroelastic interactions of flight aerodynamics and the flexible aircraft structures within flight. These aeroelastic interactions can potentially degrade aerodynamic efficiency, and thus must be accurately modeled and analyzed.

A NASA conceptual study titled “Elastically Shaped Future Air Vehicle Concept” was conducted in 2010 to investigate the potential benefits of several advanced concepts of a transport aircraft. The study showed that there exists potential benefits in shaping wing surface aeroelastic deformation actively in flight with control. In designs where structural flexibility is lessened, active wing shaping control can be used to tailor a wing’s aeroelastic shape.

The results of the study, however, also showed that conventional flap and slat devices are not aerodynamically ideal as control surfaces for active wing shaping control.

A novel control surface known as the Variable Camber Continuous Trailing Edge Flap (VCCTEF) system was proposed as a new control surface candidate. Under the Fixed Wing project Active Aeroelastic Shape Control (AASC) element, NASA and Boeing are currently conducting a joint study to investigate the application of the VCCTEF 2, 3 system on a commercial transport class aircraft. The goal of the VCCTEF study is to investigate the applicability ∗ Engineer, Intelligent Systems Division, eric.b.ting@nasa.gov † Engineer, Intelligent Systems Division, sonia.lebofsky@nasa.gov ‡ Research Scientist, Intelligent Systems Division, nhan.t.nguyen@nasa.gov, AIAA Associate Fellow § Engineer, Intelligent Systems Division, khanh.v.trinh@nasa.gov 1 of 41 American Institute of Aeronautics and Astronautics of the VCCTEF as a method to optimize the wing’s spanwise twist shape to establish the best lift-to-drag (L/D) ratio during any point within the flight envelope. This offers a significant advantage over the majority of conventional commercial aircraft wing designs which are twisted for a set cruise condition and cannot be retailored within flight.

The VCCTEF is implemented on a model of the GTM with structural flexibility of the wing considered, herein referred to as the Elastically Shaped Aircraft Concept (ESAC). As investigation of the VCCTEF system continues, wind tunnel testing as a method to gauge the potential of the new control surface has been proposed. In the current efforts, NASA and Boeing have joined together with the University of Washington Aeronautical Laboratory (UWAL) to evaluate the VCCTEF in a subsonic wind tunnel test. A procedure for modeling the development of a wind tunnel model configuration using software and numerical tools is developed in order to facilitate decision making with regard to wind tunnel testing.

This paper describes the approach for analyzing and scaling the full-scale ESAC wing to model a wind tunnel model configuration prior to inclusion of the VCCTEF. A static aeroelastic model is developed based on the ESAC wing’s jig-shape for the candidate wind tunnel model. The model is based on mimicking the aeroelastic behavior of full-scale ESAC wing, but higher wing tip deflection is desired of approximately 10% of the wing semi-span. The 5, 6 model utilizes a one-dimensional structural model of the the wing structure as a beam in coupled bending-torsion.

6–9 Finite-element method (FEM) will be used to formulate a discretized weak-form solution to the structural equations.

Vortex-lattice will be used to conduct the aerodynamic modeling and determine the loads over the wing surface. FEM and vortex-lattice are coupled together in structural-aerodynamic loops to generate the aeroelastic model.

Design of the wind tunnel model configuration is completed with optimization of the wing twist of the undeformed wind tunnel model, or wing pre-twist, along the span. This optimization is conducted so that the wind tunnel model will experience maximum L/D or minimum induced drag at the wind tunnel test condition. A gradient-based optimization using the conjugate directions method and one-dimensional line searches is applied. The resulting wind tunnel model is thus a scaled static aeroelastic clean wing model designed for the wind tunnel conditions.

II. Elastically Shaped Aircraft Concept

The elastically shaped aircraft concept (ESAC) is modeled as a notional single-aisle, mid-size, 200-passenger aircraft based on the NASA Generic Transport Model configuration. The GTM is a research platform that includes a wind tunnel model and a remotely piloted vehicle, and the geometry of the ESAC is obtained by scaling up the GTM wind tunnel model geometry by a scale of 200:11. Figure 1 is an illustration of the GTM planform. The reason for selecting the GTM as a starting point is that an extensive wind tunnel aerodynamic database exists that can be used in analysis. The benchmark configuration represents one of the most common types of transport aircraft in the commercial aviation section that provides short-to-medium range passenger carrying capacities.

Figure 1. GTM Planform In the aeroelastic model of the ESAC, the wing is allowed to freely deform based on reference B757 wing stiffness 2 of 41 American Institute of Aeronautics and Astronautics values and the GTM jig-shape planform.

As an active wing shaping control surface, the VCCTEF is implemented on the ESAC jig-shape wing. The VC- CTEF is divided into 14 sections attached to the outer wing and three sections attached to the inner wing for each side, as shown in Fig. 2. Each 24-inch section has three camber flap segments that can be individually commanded, as shown in Fig. 3. These camber flaps are joined to the neighboring sections by using a flexible and supported material (shown in blue), which deforms and provides smooth transitions between flap sections without drag producing gaps present on most control surfaces of existing aircraft. The flexible skin materials that cover the spanwise camber flap sections also constrain the flap deflections such that the relative flap deflections between any two adjacent spanwise 2, 3 flap sections are limited. More information on the VCCTEF concept can be found in references.

Figure 2. ESAC with VCCTEF Figure 3. Variable Camber Flap A. UWAL Wind Tunnel Test By early 2013, UWAL had begun investigation into the construction of a wind tunnel wing concept equipped with the VCCTEF configuration. This tunnel test was aimed to analyze the behavior of a highly flexible sub-scale model and gauge the usage of the VCCTEF as a control surface. A subsonic wind tunnel test at the UWAL facilities was planned. A notional diagram representing how the proposed test will be conducted is shown in Fig. 4. The proposed wind tunnel test is a floor-mounted test where the semi-span wing model is limited based on the height of the wind tunnel test chamber at approximately six feet. For a model of this size, the wind tunnel dynamic pressure is chosen to lb f be q = 20 .

∞ , w ft 3 of 41 American Institute of Aeronautics and Astronautics Figure 4. UWAL Wind Tunnel Test Concept As joint efforts between NASA, Boeing, and UWAL proceeded, candidate wind tunnel planforms needed to be designed. A preliminary wind tunnel model concept with five sections of VCCTEF is represented in Fig. 5. Because construction of the 17 VCCTEF sections (14 on outer wing, three on inner wing) would be intensive on a sub-scale model, a representation of the VCCTEF would be used instead. Instead of construction of 17 segments, the number is reduced and the inner high lift sections are not included.

Figure 5. UWAL Wind Tunnel Model Concept 4 of 41 American Institute of Aeronautics and Astronautics B. Wing Alone Model As a response to the need to design a wing model to meet the UWAL wind tunnel test requirements, an aeroelastic model of a full-scale “wing alone” design is initially generated. Because a candidate wind tunnel model configuration does not possess any airframe structures other than the wing, a symmetric model based on GTM geometry that has no fuselage, tails, engine nacelles, or pylons is created. The symmetry of the model is only maintained to facilitate aerodynamic modeling.

A water-tight geometry is generated by creating a new airfoil section at the fuselage-wing intersection station and shifting the jig-shape wing root to the aircraft centerline. The average location of the intersection between the fuselage wing box and the wing on the ESAC is estimated to be 6 . 1708 ft along the Body Butt Line (BBL) from the aircraft ◦ centerline. The dihedral of the ESAC wing is preserved, and the new wing alone model still possesses a Γ = 5 dihedral.

Figure 6. Vortex-Lattice Model of the Wing Alone Model As an idealization, the curved wing tip of the original jig-shape geometry is also removed and replaced with an “idealized straight wing tip”. This is not expected to have much impact on the overall wing alone model’s aerodynam- ics and is conducted only to prevent any issues with the vortex-lattice modeling and analysis.

Figure 7. Wing Alone Model Idealized Straight Wing Tip The wing alone reference area is obtained by integrating the chord over the semi-span and multiplying by a factor of 2. Let y represent the coordinate on the wing along the aircraft pitch axis running from the aircraft centerline B b outwards towards the right wing tip (seen in Fig. 8 as the b b direction). The subscript f is used to refer to quantities related to the full-scale wing alone model.

b ˆ f S = 2 c ( y ) dy ≈ 1640 . 8 ft (1) B B re f , f In the wing alone model, the wing reference area does not include the “fictitious” area that would have been covered by the fuselage of the full-scale wing generally included in trapezoidal area estimate of a wing.

The mean aerodynamic chord is also obtained through integration.

b f ˆ ¯ c = c ( y ) dy ≈ 17 . 0991 ft (2) f B B S re f , f 0 The wing aspect ratio is determined using the span and reference area.

b f AR = = 7 . 6000 (3) f S re f , f 5 of 41 American Institute of Aeronautics and Astronautics The taper ratio is obtained using the root and tip chord values.

c 5 . 5990 ft t , f λ = = ≈ 0 . 1950 (4) f c 28 . 7122 ft r , f

III. Wing Structural Modeling

A structural model of the wing using beam theory is developed which is later incorporated into a fully coupled structural-aerodynamic aeroelasticity model.

A. Reference Frames Figure 8. Aircraft Reference Frames Figure 8 illustrates three orthogonal views for a typical aircraft and several associated reference frames. These refer- ence frames are useful in developing the structural models of the lifting surfaces of an aircraft, although the coordinate frames associated with the aircraft wings are primarily used in this analysis. The aircraft body-fixed reference frame B is defined by the unit vectors b b b , b b b , and b b b which are aligned with the aircraft roll, pitch, and yaw axes, respectively.

1 2 3 The reference frame C is aligned with the right wing’s elastic axis and is defined by the unit vectors c c c , c c c , and c c c .

1 2 3 Let Λ be the sweep of the elastic axis. The B frame can be related to C through three successive rotations: 1) the first ′ π ′ rotation about b b b by an angle of + Λ to generate an intermediate reference frame B defined by the unit vectors b b b , ′ ′ ′ b b b , and b b b (not shown), 2) the second rotation about b b b by the dihedral angle Γ of the elastic axis that results in the 2 3 2 ′ ′ ′ ′ intermediate reference frame C defined by the unit vectors c c c , c c c , and c c c (not shown), and 3) the third rotation about 1 2 3 ′ c c c by an angle of π to result in the reference frame C . The transformation can be represented by a series of coordinate rotations expressed as ⎡ ⎤ ⎡ ⎤ ⎡ ⎤ ⎡ ⎤ ⎡ ⎤ b b b − sin Λ − cos Λ 0 cos Γ 0 sin Γ 1 0 0 c c c 1 1 ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ = ⎣ b b b ⎦ ⎣ cos Λ − sin Λ 0 ⎦ ⎣ 0 1 0 ⎦ ⎣ 0 − 1 0 ⎦ ⎣ c c c ⎦ 2 2 b b b 0 0 1 − sin Γ 0 cos Γ 0 0 − 1 c c c 3 3 ⎡ ⎤ ⎡ ⎤ − sin Λ cos Γ cos Λ sin Λ sin Γ c c c ⎢ ⎥ ⎢ ⎥ = (5) ⎣ cos Λ cos Γ sin Λ − cos Λ sin Γ ⎦ ⎣ c c c ⎦ − sin Γ 0 − cos Γ c c c The analysis can be repeated for the left wing. The reference frame D is aligned with the left wing’s elastic axis and is defined by the unit vectors d d d , d d d , and d d d . The B frame can be related to D through three successive rotations: 1 2 3 6 of 41 American Institute of Aeronautics and Astronautics π ′′ 1) the first rotation about − b b b by an angle of + Λ to generate an intermediate reference frame B defined by the ′′ ′′ ′′ ′′ b b b b unit vectors b b , b b , and b b (not shown), 2) the second rotation about b b by the dihedral angle Γ of the elastic axis that 1 2 3 2 ′ ′ ′ ′ d d d results in the intermediate reference frame D defined by the unit vectors d d , d d , and d d (not shown), and 3) the third 1 2 3 ′ rotation about d d d by an angle of π to result in the reference frame D . The relationship can be expressed as ⎡ ⎤ ⎡ ⎤ ⎡ ⎤ ⎡ ⎤ ⎡ ⎤ b b b − sin Λ cos Λ 0 cos Γ 0 sin Γ 1 0 0 d d d 1 1 ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ = ⎣ b b b ⎦ ⎣ − cos Λ − sin Λ 0 ⎦ ⎣ 0 1 0 ⎦ ⎣ 0 − 1 0 ⎦ ⎣ d d d ⎦ 2 2 b b b 0 0 1 − sin Γ 0 cos Γ 0 0 − 1 d d d 3 3 ⎡ ⎤ ⎡ ⎤ − sin Λ cos Γ − cos Λ sin Λ sin Γ d d d ⎢ ⎥ ⎢ ⎥ = (6) ⎣ − cos Λ cos Γ sin Λ cos Λ sin Γ ⎦ ⎣ d d d ⎦ − sin Γ 0 − cos Γ d d d B. Elastic Axis An analysis of the combined motion of the left wing is conducted in the present section, and the motion of the right wing is considered to be equivalent for symmetric flight. This analysis is equivalent to that in a previous study and is included for completeness.

Let x represent the coordinate along the elastic axis of a wing running from root to tip. The wing pre-twist angle γ ( x ) thus represents the incidence of the airfoil section at the corresponding elastic axis coordinate. A typical wing pre- twist varies from nose-up at the wing root to nose-down at the wing tip and is commonly referred to as a “wash-out” twist distribution.

The internal structure of a wing is typically composed of a complex arrangement of load carrying spars and wing boxes that carry the stresses and strains introduced by aerodynamic forces and aeroelastic deflections. For this analysis, an equivalent beam approach is used which models the wing’s elastic behavior using equivalent stiffness properties. It is a common approach in analyzing aeroelastic deflections and can be used to analyze high aspect ratio wings with good accuracy. The effect of wing curvature is ignored and straight beam theory is used to model the wing deflection.

The axial or extensional deflection of a wing is also generally very small and is neglected.

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

7 of 41 American Institute of Aeronautics and Astronautics Differentiating with respect to x gives [ ] [ ] [ ] [ ] ′ ′ y − sin γ − cos γ η − z γ x = γ = (8) ′ z cos γ − sin γ ξ y γ x 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 d d − W d d d + V d d d (9) 1 x 2 x 3 where the subscript x denotes 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. Then the coordinates ( x , y , z ) 1 1 1 1 1 1 are computed using the small angle approximation as ⎡ ⎤ ⎡ ⎤ ⎡ ⎤ ⎡ ⎤ φ d d d x ( x , t ) x φ φ × ( yd d + zd d ) . d d x − yV − zW 1 2 3 1 x x ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ = + φ d d d = (10) ⎣ y ( x , t ) ⎦ ⎣ y + V ⎦ ⎣ φ φ × ( yd d + zd d ) . d d ⎦ ⎣ y + V − z Θ ⎦ 1 2 3 2 φ d d d z ( x , t ) z + w φ φ × ( yd d + zd d ) . 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 ⎢ ⎥ ⎢ ′ ′ ⎥ = (11) ⎣ 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 , x ε = = − 1 (12) ds s x where √ √ ′ 2 2 2 2 2 s = 1 + y + z = 1 + ( y + z )( γ ) (13) 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 γ ) (14) xx xx x 1 , x 1 , x 1 , x x Ignoring the second-order terms and using Taylor series expansion, s is approximated as 1 , x ′ 2 2 − yV − zW + ( y + z ) γ Θ xx xx x s ≈ s + 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 ) γ y 1 + ( y + z )( γ ) Θ (15) xx xx x ′ For a small wing twist angle γ , ( γ ) ≈ 0. Then longitudinal strain can be expressed as ′ 2 2 ε = − yV − zW + ( y + z ) γ Θ (16) xx xx x The moments acting on the wing are then obtained as ⎡ ⎤ ⎡ ⎤ ⎡ ⎤ ′ 2 2 ˆ ˆ M GJ Θ ( y + z )( γ + Θ ) x x x ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ = + E ε dydz ⎣ M ⎦ ⎣ 0 ⎦ ⎣ − z ⎦ y M 0 − y z ⎡ ⎤ ⎡ ⎤ ′ ′ ′ GJ + EB ( γ ) − EB γ − EB γ Θ 1 2 3 x ⎢ ′ ⎥ ⎢ ⎥ = ⎣ − EB γ EI − EI ⎦ ⎣ W ⎦ (17) 2 yy yz xx ′ − EB γ − EI E V 3 yz Izz xx 8 of 41 American Institute of Aeronautics and Astronautics ′ 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 (18) ⎣ 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 as in highly twisted wings such as turbomachinery blades.

C. Aeroelastic Angle of Attack The aeroelastic angle of attack is the effective angle of attack of a flexible wing section that is undergoing aeroelastic deformation defined by elastic axis twist Θ , flapwise bending W , and chordwise bending V . It can be calculated by solving for the relative velocity of air as it approaches a wing section perpendicular to the elastic axis. The aeroelastic angle of attack encompasses a wing section’s rigid local angle of attack and the contribution due to wing elastic deformation, and also governs the aerodynamic forces and moments on a local wing section.

The local angle of attack depends on the relative approaching air velocity, the rotation angle φ φ φ , and 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 v = ¯ v v v + ω × r r r = ( ub b b + vb b b + wb b b ) + ( pb b b + qb b b + rb b b ) × ( − x b b b − y b b b − z b b b ) Q 1 2 3 1 2 3 a 1 a 2 a 3 = ( u + ry − qz ) b b b + ( v − rx + pz ) b b b + ( w + qx − py ) b b b (19) a a 1 a a 2 a a 3 where ( x , y , z ) are the coordinates of point Q in the aircraft body B frame relative to the aircraft center of gravity a a a (C.G.) such that x is positive when point Q is aft of the aircraft C.G., y is positive when point Q is towards the left a a wing from the aircraft C.G., and z is positive when point Q is above the aircraft C.G. The aircraft velocity is ( u , v , w ) a in the aircraft body axes, and ( p , q , r ) are the aircraft angular velocity components.

This can be expressed in the left wing frame D as v v v = x d d d + y d d d + z d d d , the local velocity due to aircraft Q t 1 t 2 t 3 rigid-body dynamics.

⎡ ⎤ ⎡ ⎤ x − ( u + ry − qz ) sin Λ cos Γ − ( v − rx + pz ) cos Λ cos Γ − ( w + qx − py ) sin Γ t a a a a a a ⎢ ⎥ ⎢ ⎥ = (20) ⎣ 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 For a trim case where β = 0, p = q = r = 0, then ⎡ ⎤ ⎡ ⎤ x − u sin Λ cos Γ − w sin Γ t ⎢ ⎥ ⎢ ⎥ = (21) ⎣ y ⎦ ⎣ − u cos Λ ⎦ t z u sin Λ sin Γ − w cos Γ t The local velocity at point Q due to both 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 v = v v v + φ φ φ × p p p = v d d d + v d d d + v d d d = (22) d d d d d d d d d ⎣ y + V − ( yV + zW ) V − ( z + W + y Θ ) Θ ⎦ Q x 1 y 2 z 3 1 2 3 t t x x xt t z + W − ( yV + zW ) W + ( y + V − z Θ ) Θ t t x x xt t where ( x , y , z ) are the coordinates for the point Q in the reference frame D without any aeroelastic deflection.

For static aeroelasticity, all the velocity components of the aeroelastic deflections are set to zero. Thus, v = x , x t v = y , and v = z .

y t z 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 ( μ , η , ξ ) . The transformation can be performed 9 of 41 American Institute of Aeronautics and Astronautics using successive rotation matrix multiplication operations as ⎡ ⎤ ⎡ ⎤ ⎡ ⎤ ⎡ ⎤ ⎡ ⎤ v 1 0 0 cos V sin V 0 cos W 0 sin W x μ x x x x t ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ v = ⎣ ⎦ ⎣ 0 cos ( Θ + γ ) sin ( Θ + γ ) ⎦ ⎣ − sin V cos V 0 ⎦ ⎣ 0 1 0 ⎦ ⎣ y ⎦ η x x t v 0 − sin ( Θ + γ ) cos ( Θ + γ ) 0 0 1 − sin W 0 cos W z x x t ξ ⎡ ⎤ 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 ) z x x x z x y x x x z x t ⎤ ⎡ v + v V + v W x y x z x ⎢ ⎥ ≈ ⎣ − v [ V + W ( Θ + γ )] + v + v [( Θ + γ ) − V W ] ⎦ (23) 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 c η ξ reference frame D is computed as v + w ¯ v + Δ v + w v + w ( ¯ v + w ) Δ v i i i i η ξ ξ ξ ξ ξ α = = = − (24) c v ¯ v + Δ v ¯ v ¯ v η η η η η where w is the downwash due to the three-dimensional lift distribution over a finite-aspect ratio wing.

i The velocity components are ¯ v = u sin Λ sin Γ − w cos Γ (25) ξ Δ v = v [ − W + V ( Θ + γ )] − v ( Θ + γ ) (26) ξ x x x y ¯ v = − u cos Λ (27) η Δ v = − v [ V + W ( Θ + γ )] + v [( Θ + γ ) − V W ] (28) η x x x z x x The local aeroelastic angle of attack is evaluated as u sin Λ sin Γ − w cos Γ + ( − u sin Λ cos Γ − w sin Γ ) [ − W + V ( Θ + γ )] + u cos Λ ( Θ + γ ) + w x x i α = − c u cos Λ { u sin Λ sin Γ − w cos Γ + w i − − ( − u sin Λ cos Γ − w sin Γ ) [ V + W ( Θ + γ )] x x 2 2 u cos Λ } + ( u sin Λ sin Γ − w cos Γ ) [( Θ + γ ) − V W ] (29) x x Assuming a trim case, let u ≈ V and w ≈ V α . Neglecting chordwise bending components, V , V , and also ∞ ∞ x neglecting the three-dimensional finite-wing effect, w , allows α be expressed as i c sin Λ sin Γ α cos Γ sin Λ cos Γ + α sin Γ α = − + − W − Θ − γ c x cos Λ cos Λ cos Λ { } sin Λ sin Γ − α cos Γ − ( u sin Λ cos Γ + w sin Γ ) ( W Θ + W γ ) + ( u sin Λ sin Γ − w cos Γ ) ( Θ + γ ) (30) x x u cos Λ sin Λ sin Γ α cos Γ sin Λ cos Γ + α sin Γ α = − + − W − Θ − γ c x cos Λ cos Λ cos Λ ( ) 2 2 2 − sin Λ sin Γ cos Γ + α sin Λ cos Γ − α sin Λ sin Γ + α sin Γ cos Γ + + ( W Θ + W γ ) x x 2 2 cos Λ cos Λ ( ) 2 2 2 2 − sin Λ sin Γ + α sin Λ sin Γ cos Γ α sin Λ sin Γ cos Γ − α cos Γ + + ( Θ + γ ) (31) 2 2 cos Λ cos Λ Eliminating higher order terms results in the aeroelastic angle of attack expressed as cos Γ α = − γ − tan Λ sin Γ + α − Θ − W tan Λ cos Γ (32) c x cos Λ 10 of 41 American Institute of Aeronautics and Astronautics This can be re-expressed after applying small angle approximation in terms of partial derivatives as ∂ α ∂ α ∂ α ∂ α c c c c α ( x , y , z ) = + α + W + Θ (33) c x ∂ 1 ∂ α ∂ W ∂ Θ x ∂ α c = − γ − tan ΛΓ (34) ∂ 1 ∂ α 1 c = (35) ∂ α cos Λ ∂ α c = − tan Λ (36) ∂ W x ∂ α c = − 1 (37) ∂ Θ The aeroelastic deflections terms W and Θ contribute to aerodynamic stiffness.

x D. Coupled Bending-Torsion Equations Without considering chordwise bending of the wing, the equilibrium conditions for bending and torsion are expressed as ∂ M x = − m (38) x ∂ x ∂ M ∂ m y y = f − (39) z ∂ x ∂ x where m is the pitching moment per unit span about the elastic axis, f is the lift force per unit span, and m is the x z y bending moment per unit span about the flapwise axis of the wing.

Because the structural modeling is intended for use in a static aeroelasticity model, a steady-state aerodynamics model is used. Aerodynamic information can be obtained through vortex-lattice modeling to develop the forces and moments for coupled bending-torsion of a flexible wing.

Figure 10. Airfoil Forces and Moments Neglecting the effect of downwash that is caused due to lift generation over a three-dimensional finite-wing, the lift coefficient over the span of a clean wing assuming linear aerodynamics is as follows: c ( x ) = c ( x ) α ( x ) (40) L L c α where α is the aeroelastic angle of attack as shown in Fig. 10, assumed to be constant for airfoil cross sections c perpendicular to the elastic axis and only a function along the wing such that α = α ( x ) .

c c The aeroelastic angle of attack α can be expressed as contributions due to aircraft rigid-body angle of attack α as c well as the contribution due to aeroelastic deformation. Let α represent the contribution the aeroelastic angle of attack r due to rigid-body considerations including airfoil shape, and α represent the effect on the local aeroelastic angle of e attack due to aeroelastic deformation at the aerodynamic center of the airfoil section. Based on Eq. 33 and neglecting chordwise bending and assuming dihedral is small, α and α can be represented as r e ( ) ∂ α ∂ α c c α ( x ) = ( x ) + ( x ) α ( x ) cos Λ (41) r ∂ 1 ∂ α 11 of 41 American Institute of Aeronautics and Astronautics ( ) ∂ α ∂ α c c α ( x ) = ( x ) W ( x ) + ( x ) Θ ( x ) cos Λ (42) e x ∂ W ∂ Θ x α ( x ) + α ( x ) r e α ( x ) = (43) c cos Λ where both α and α are about the pitch axis direction, positive nose-up.

r e The rigid and elastic lift coefficient contributions to the local sectional lift coefficient of the elastic axis airfoil cross sections can be expressed as α ( x ) r c ( x ) = c ( x ) (44) L L r α cos Λ α ( x ) e c ( x ) = c ( x ) (45) L L e α cos Λ It is also important to note that the elastic contribution α to the local aeroelastic angle of attack α can be repre- e c sented based on the partial derivatives calculated in Eqs. 36 and 37. Given a deformation characterized by elastic axis twist Θ and vertical bending slope W , the elastic contribution to the aeroelastic angle of attack can be calculated as x α ( x ) = − Θ ( x ) cos Λ − W ( x ) sin Λ (46) e x where α is about the aircraft pitch axis. This applies for the static case using the assumptions and simplifications e applied in derivation of α .

c The steady-state drag coefficient can be modeled by a parabolic drag polar as c ( x ) = c ( x ) + k ( x ) c ( x ) (47) D D L where c is the section parasitic drag coefficient and k is the section drag polar parameter.

D Likewise, the pitching moment coefficient about the aircraft pitch axis can be represented as e ( x ) c ( x ) = c ( x ) + c ( x ) cos Λ (48) m m L ac c ( x ) where e is the location of the aerodynamic center relative to the elastic axis along the body axis, positive when the aerodynamic center is forward of the elastic axis, and c is defined about the pitch axis, positive nose-up.

m ac Expanding Eq. 48 using the aeroelastic angle of attack α definition in Eq. 33 produces e ( ) e ( x ) ∂ α ∂ α ∂ α ∂ α c c c c c ( x ) = c ( x ) + c ( x ) + α + Θ + W cos Λ (49) m m L x ac α c ( x ) ∂ 1 ∂ α ∂ Θ ∂ W x which allows us to define the following quantity ( ) e ( x ) ∂ α ∂ α c c c = c ( x ) + c ( x ) + α cos Λ (50) m m L ac α c ( x ) ∂ 1 ∂ α The lift force, drag force, and pitching moment about the aircraft pitch axis are expressed as l = c q cos Λ c (51) L ∞ d = c q cos Λ c (52) D ∞ m = c q c (53) m ∞ where cos Λ takes into account the correction due to the elastic axis sweep but is not needed in the pitch moment calculation since Eq. 48 is already about the pitch axis.

The forces and moments in the local coordinate reference frame are obtained as a f = ( l cos α + d sin α ) Γ + ( d cos α − l sin α ) sin Λ (54) x a f = ( d cos α − l sin α ) cos Λ (55) y a f = l cos α + d sin α − ( d cos α − l sin α ) sin ΛΓ (56) z a m = − m cos Λ (57) x 12 of 41 American Institute of Aeronautics and Astronautics a m = m sin Λ (58) y a m = m cos ΛΓ (59) z For a model with only flapwise bending and torsion considered, the beam deflection analysis the only aerodynamic a a a force and moment terms that have an effect are the terms f , m , and m . For this analysis, the aerodynamic force and z x y moment terms are thus considered to be a 2 f ≈ c q cos Λ c (60) L ∞ z a 2 2 m ≈ − c q cos Λ c (61) m ∞ x a ∂ m ∂ c y m ≈ q sin Λ cos Λ c (62) ∞ ∂ x ∂ x where an additional cos Λ term considers the change in the direction of q over the wing due to sweep.

∞ Inserting Eq. 17 and the force and moment terms Eqs. 60-62 into the governing equilibrium equations Eqs. 38 and 39, the following equations can be used to describe the coupled bending and torsion motion of the wing: ∂ ′ ∂ c m 2 2 2 ( − EB γ Θ + EI W ) = − mW + me Θ + c q cos Λ c − q tan Λ cos Λ c (63) 2 x yy xx tt cg tt L ∞ ∞ ∂ x ∂ x {[ ] } ∂ ′ ′ 2 2 2 2 GJ + EB ( γ ) Θ − EB γ W = mr Θ − me W + c q cos Λ c (64) 1 x 2 xx tt cg tt m ∞ k ∂ x Although the wing alone model is modeled as a symmetric wing in a horizontal plane, the actual wing tunnel model will be only a semi-span wing mounted in a vertical plane. Gravitational forces on the wing alone and candidate wind tunnel model are thus ignored. There are also no engines on the model, and thus the only forces and moments being considered in the aeroelastic model are from aerodynamic and inertial sources.

IV. Finite-Element Modeling

The development of the coupled bending-torsion partial differential equations describing the wing allows for wing bending and torsional deflections to be solved. FEM is used as a numerical technique that uses locally-defined basis functions to numerically approximate the solution of the governing partial differential equations. The FEM is used to discretize the wing structure into n equally spaced one-dimensional elements. The bending and torsional deflections can be approximated as n Θ ( x , t ) = Θ ( x , t ) (65) i

∑

i = 1 n W ( x , t ) = W ( x , t ) (66)

∑ i

i = 1 where i refers to the i -th element.

For each element, the bending and torsional deflections are approximated as [ ] [ ] θ ( t ) i Θ ( x , t ) = ψ θ ( t ) + ψ ( x ) θ ( t ) = = N ( x ) θ ( t ) (67) ψ ( x ) ψ ( x ) i i 1 2 2 θ i i i 1 2 θ ( t ) i [ ] ′ ′ W ( x , t ) = φ ( x ) w ( t ) + φ ( x ) w ( t ) + φ ( x ) w ( t ) + φ ( x ) w ( t ) i 1 1 2 3 2 4 i 1 i 2 i i ⎡ ⎤ w ( t ) i ′ [ ] ⎢ ⎥ w ( t ) ⎢ ⎥ i = = N ( x ) w ( t ) (68) φ ( x ) φ ( x ) φ ( x ) φ ( x ) ⎢ ⎥ w i 1 2 3 4 ⎣ ⎦ w ( t ) i ′ w ( 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 − (69) l 13 of 41 American Institute of Aeronautics and Astronautics x ψ ( x ) = (70) l x x 2 3 φ ( x ) = 1 − 3 ( ) + 2 ( ) (71) l l [ ] x x x 2 3 φ ( x ) = l − 2 ( ) + ( ) (72) l l l x x 2 3 φ ( x ) = 3 ( ) − 2 ( ) (73) l l [ ] x x 2 3 φ ( x ) = l − ( ) + ( ) (74) l l L where x ∈ [ 0 , 1 ] is the local coordinate and l = is the element length.

n The weak-form integral expressions of the coupled bending-torsion partial differential equations are obtained by T T multiplying the equations by N ( x ) and N ( x ) and then integrating over the wing span. The aerodynamic coefficients w θ are expanded here based on the aeroelastic angle of attack representation in Eq. 33 and using Eq. 49. This yields ˆ {[ ] } n l d ′ ′ ′ ′′ T 2 N GJ + EB ( γ ) N θ − EB γ N w dx =

∑ 1 i 2 i

θ θ w dx i = 0 [ ( )] ˆ ˆ n l n l e ∂ α ∂ α ′ c c T 2 T 2 2 ¨ N mr N θ − me N ¨ w ) dx + N c + cos Λ c N θ + N w q cos Λ c dx (75) θ i cg w i m L θ i i ∞

∑ θ k ∑ θ 0 α w

c ∂ Θ ∂ W 0 0 x i = 1 i = 1 ˆ n l 2 d ′ ′ ′′ T N ( − EB γ N θ + EI N w ) dx = 2 i yy i

∑ w θ w

dx i = 0 ˆ [ ( )] n n l ∂ α ∂ α ∂ α ∂ α ′ c c c c T T 2 ¨ N ( ρ AN ¨ w + ρ Ae N θ ) dx + N c + + N θ + N w q cos Λ cdx w i cg θ i L θ i i ∞

∑ w ∑ w α w

∂ 1 ∂ α ∂ Θ ∂ W x i = 1 i = 1 ˆ [ ( )] n l d e ∂ α ∂ α ′ c c T 2 2 − N c + cos Λ c N θ + N w tan Λ q cos Λ c dx (76) m L θ i i ∞

∑ w α w

dx c ∂ Θ ∂ W x i = 1 The expressions of the left hand sides can be integrated by parts upon enforcing the boundary conditions, resulting in ˆ ˆ {[ ] } {[ ] } l l d ′ ′ ′ ′′ ′ ′ ′ ′ ′′ T 2 T 2 N GJ + EB ( γ ) N θ − EB γ N w dx = − N GJ + EB ( γ ) N θ − EB γ N w dx (77) 1 i 2 i 1 i 2 i θ θ w θ θ w dx 0 0 ˆ ˆ l 2 l d ′ ′ ′′ ′′ ′ ′ ′′ T T N ( − EB γ N θ + EI N w ) dx = N ( − EB γ N θ + EI N w ) dx (78) 2 i yy i 2 i yy i w θ w w θ w dx 0 0 The elemental mass matrix, stiffness matrices, and force vector are then established as [ ] ˆ l 2 T T r N N − e N N θ cg w k θ θ M = m dx (79) s i T T − e N N N N 0 cg θ w w w ] [ [ ] ˆ ′ ′ ′ ′ ′ ′′ l 2 T T GJ + EB ( γ ) N N − EB γ N N 1 2 w θ θ θ K = dx (80) s i ′ ′′ ′ ′′ ′′ T T − EB γ N N EI N N 2 yy w w w θ [ ] ˆ ′ ∂ α ∂ α l c T c T ec cos Λ N N ec cos Λ N N L θ L α α w ∂ Θ θ ∂ W θ 2 x K = q cos Λ cdx (81) a ′ ∞ i ∂ α ∂ α c T c T − c N N − c N N 0 L θ L α w α w w ∂ Θ ∂ W x ([ ] [ ]) ˆ l T − c cN 0 m 0 θ F = q cos Λ c + dx r ∞ ′ i T T c N − c c tan Λ N 0 L m w w [ ]∣ l ∣ ∣ + q cos Λ c (82) ∣ ∞ T ∣ − c c tan Λ N m w 14 of 41 American Institute of Aeronautics and Astronautics where K is the structural stiffness matrix and K is the stiffness matrix due to the result of aerodynamics. The moment s a component due to lift eccentricity that contributes to bending terms in K are neglected.

a The globally assembled system is described by the matrix equation M ¨ x + K x = F − K x (83) s e s e r a e [ ] T ′ ′ ′ ′ where x = .

e θ w w θ w w . . . θ w w . . . θ w w 1 1 2 2 i i n + 1 n + 1 1 2 i n + 1 Equation 83 represents the governing equation for solving the structural deflection of a flexible wing given aero- dynamic force and moment inputs. By setting ˙ x = 0 0 0, the equilibrium solution for x can be obtained through inverting e e the stiffness matrix and pre-multiplying the force matrix. Information can then be extracted from the solution including the wing’s deflection along the elastic axis of the wing Θ and W . It can also be used to calculate the elastic contribution to the aeroelastic angle of attack α in Eq. 33.

e Without considering the effect of aeroelastic coupling that results in the aerodynamic stiffness matrix K , a wing a structural deflection can be solved by − 1 x = K F (84) e r s Aeroelastic deflection of the wing can be calculated by including the aerodynamic stiffness. The aeroelastic solu- tion is represented by the term ¯ x to differentiate it from x , which is structural deflection calculated without using the e e aerodynamic stiffness matrix.

− 1 ¯ x = ( K + K ) F (85) e s a r The term F represents the force matrix in the FEM that is constructed solely based on the aerodynamic character- r istics of the undeformed wing. Thus, Eq. 85 represents a method in which the aeroelastic deformation of a wing can be solved solely using the rigid wing aerodynamic loads and properties. While this is a powerful framework, another approach exists that involves updating the force matrix in the FEM instead of utilizing the rigid wing properties.

Let F represent the load on the wing due to the rigid planform and the incremental effect due to aeroelastic deformation. The force F can be represented as F = F − K x (86) r a e where x is the deflection of the wing. Aeroelastic deformations can be calculated by an iterative technique where e k + 1 − 1 k k + 1 x = K F , x → ¯ x as k increases (87) e e s e k where the force vector F is estimated using computational aerodynamic modeling tools to calculate the loads on the deformed geometry. This approach, which is used in this study, requires an aerodynamic modeling tool to be run coupled with the FEM model. One of the advantages of using this method is that the aircraft and sectional characteristics are estimated for the deformed geometry, instead of using the assumption that they remain constant at the rigid wing values. This can also provide a better estimate of the effect of K than the analytical model developed.

a

V. Vortex-Lattice Aerodynamic Modeling

Vorview is a computational tool used for aerodynamic modeling of aircraft configurations using vortex-lattice method. Based on lifting line aerodynamic theory, Vorview provides a rapid method for estimating aerodynamic force and moment coefficients. Geometric input vehicle configuration are constructed within Vorview by discretizing the surface into a series of panels, which are then represented by placement of spanwise and chordwise locations of bound or horseshoe vortices. Vorview computes the vehicle aerodynamics in both the longitudinal and lateral directions independently, and these can be combined to produce the overall aerodynamic characteristics of the vehicle at any arbitrary angle of attack and angle of sideslip.

Vorview is considered a medium fidelity tool, and limitations associated with vortex-lattice modeling in general apply to Vorview aerodynamic analysis. The drag prediction by Vorview is most reliable only for induced drag prediction due to the inviscid nature of any vortex-lattice method. Prediction of viscous drag due to boundary layer separation and wave drag due to shock-induced boundary layer separation are generally not conducted by vortex- lattice, and viscous drag must be estimated using other methods.

In addition to force and moment analysis, Vorview can provide a rapid estimation of aerodynamic derivatives including dynamic derivatives due to angular rates. These aerodynamic stability and control derivatives are useful in analyzing the stability and handling characteristics of an aircraft configuration. Owing to the computationally 15 of 41 American Institute of Aeronautics and Astronautics 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 using the results from these stability 4 13 and handling analyses. 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.

Figure 11. GTM Aircraft Model in Vorview In this study, Vorview will be utilized as the primary tool for conducting aerodynamic modeling for the aircraft configurations. Total aircraft characteristics as well as sectional data along the aircraft wing surfaces can be post- processed from Vorview.

VI. Automated Geometry Modeling Tool

An automated geometry generation tool is developed in Matlab that is used to close the loop between the structural and aerodynamic modeling needed to generate an aeroelastic model. The geometry generation tool uses structural deflection data that is computed by the FEM model and applies it to the undeformed aircraft wing geometry to reflect static aeroelastic deflections. The vehicle geometry modeler directly outputs a geometry input file that can be read by Vorview when computing an aeroelastic solution.

Figure 12. GTM Coordinate Systems Consider the reference frames in Fig. 12. The coordinate reference frame ( x , y , z ) defines the Body Station A A A (BS), the Body Butt Line (BBL), and the Body Water Line (BWL) of the aircraft, respectively. The coordinate reference frame ( x , y , z ) is the translated coordinate system attached to the nose of the aircraft such that x = V V V V 16 of 41 American Institute of Aeronautics and Astronautics x − 13 . 25 ft, y = y , and z = z − 15 . 8333 ft. This reference frame is used by vortex-lattice aerodynamic modeling B V B V B tool. The aircraft body reference frame ( x , y , z ) is the same B coordinate system defined earlier in Fig. 8 by the unit B B B vectors b b b , b b b , and b b b . The B coordinate frame is attached to the aircraft center of gravity (CG) such that x = ¯ x − x , 1 2 3 B V V y = y − ¯ y , and z = ¯ z − z , where ( ¯ x , ¯ y ¯ z ) is the coordinate of the CG in the ( x , y z ) reference frame.

B V V B V V V V , V V V , V The vehicle geometry modeler has access to the outer mold line of the aircraft geometry. It is capable of applying geometric transformations onto the outer mold coordinates of the wing’s jig-shape to simulate aeroelastic deflection.

Neglecting chordwise bending deflection and utilizing the coordinate system of the left wing developed earlier (coor- dinate frame D ), the aeroelastic deflections in bending and torsion are in expressed in a vector form as φ d d φ φ = Θ d d − W d d (88) 1 x 2 Δ Δ Δ r r r = − W sin W d d d + W cos W d d d (89) x 1 x 3 The coordinate reference frame ( x , y , z ) of the left wing is related to the coordinate reference frame ( x , y , z ) by V V V the following relationship ⎡ ⎤ ⎡ ⎤ ⎡ ⎤ d d d − sin Λ cos Γ − cos Λ sin Γ − sin Γ b b b 1 1 ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ d = b ⎣ d d ⎦ ⎣ − cos Λ sin Λ 0 ⎦ ⎣ b b ⎦ 2 2 d b d d sin Λ sin γ cos Λ sin Γ − cos Γ b b 3 3 ⎡ ⎤ ⎡ ⎤ v − sin Λ cos Γ − cos Λ cos Γ − sin Γ − v v ⎢ ⎥ ⎢ ⎥ = v (90) ⎣ − cos Λ sin Λ 0 ⎦ ⎣ v v ⎦ v sin Λ sin Γ cos Λ sin Γ − cos Γ − v v v v v where ( v v , v 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 angle of attack Δ α (positive nose-up), a horizontal deflection Δ y (positive deflection towards wing tip), and a vertical deflection Δ z (positive V V displacement upward) as follows: Δ α = − Θ cos Λ cos Γ − W sin Λ (91) x Δ y = − W sin W cos Λ cos Γ − W cos W cos Λ sin Γ (92) V x x Δ z = − W sin W sin Γ + W cos W cos Γ (93) v 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 angle of attack Δ α and then translating the resultant coordinates by the horizontal deflection Δ y and the vertical deflection Δ z .

V V Note that the transformation for Δ α is equivalent to the value of α , the local change in the angle of attack for a e wing section due to aeroelastic deformation represented by Eq. 46, when dihedral Γ is small.

VII. Static Aeroelastic Model

In a standard static aeroelastic model, it is understood that the modeling effort needs to take into account that structural deformations during flight will alter the aircraft aerodynamics, and changing the aerodynamics will thus change the structural deformations. This realizes an aeroelastic model where coupling exists between the structural modeling and aerodynamic modeling approaches. Previous studies have analytically constructed fully coupled aeroe- 5, 6 lastic finite-element models that utilize rigid wing lift-curve slopes as an aerodynamic model. This study employs a static aeroelastic model that is constructed by utilizing a structural FEM model coupled with a vortex-lattice solution.

A static aeroelastic code is developed by utilizing the automated geometry generation modeling tool to close the loop between the FEM model and the vortex-lattice model using Eq. 87. For a model considering only flapwise ¯ ¯ bending and axial torsion, the aeroelastic deflection can be summarized by the quantities of ¯ Θ ( x ) , W ( x ) , W ( x ) , the x aeroelastic elastic axis twist, aeroelastic vertical (flapwise) bending, and aeroelastic vertical bending slope respectively.

These quantities are emphasized to be aeroelastic deflections, while the terms Θ ( x ) , W ( x ) , W ( x ) are considered x structural deflection terms which may or may not be the aeroelastic solution for a given flight condition. Closing the static aeroelastic loop causes the the structural deflections Θ ( x ) , W ( x ) , W ( x ) to converge to the aeroelastic solution x ¯ Θ ( x ) , ¯ W ( x ) , ¯ W ( x ) as iterations are conducted. The structural and aeroelastic deformations can also be represented by x the elastic contribution to the aeroelastic angle of attack in Eq. 46, or α ( x ) and ¯ α ( x ) .

e e 17 of 41 American Institute of Aeronautics and Astronautics Figure 13. Static Aeroelastic Model Concept The static aeroelastic model maps an input desired ¯ C , Mach number M , and altitude h into the respective static L ¯ ¯ aeroelastic deflection ( ¯ Θ , W , W ) and the angle of attack ¯ α that the flexible wing aircraft would experience when x leveled at the desired ¯ C . The following procedure is followed: L 1. Vortex-lattice modeling is conducted on an input geometry at an input flight condition α , M to determine the aircraft total aerodynamic quantities, as well as sectional aerodynamic distributions of c ( x ) , c ( x ) , k ( x ) , L m ac c ( x ) , and x ( x ) or the location of the section aerodynamic centers.

L ac α 2. The structural FEM model uses the sectional aerodynamic inputs to calculate the wing’s structural deflection Θ ( x ) and W ( x ) .

3. The geometry generation tool converts Θ ( x ) and W ( x ) into the series of deformations α ( x ) , Δ y ( x ) , and Δ z ( x ) , e V V and generates a new aircraft geometry with the deformed wing.

4. A lift curve is generated based on the deflected wing aircraft geometry. The angle of attack α for the value ¯ C L is determined and selected for the next iteration.

5. Steps 1-4 are repeated until Δ α between iterations is within a criteria.

This converged solution is represented by the angle of attack ¯ α that the flexible wing model would need to have a lift coefficient of ¯ C . The wing shape at the converged flight condition is the converged aeroelastic deflection ¯ Θ , ¯ W , ¯ W , L x and ¯ α .

e The model is used to determine the full-scale wing alone model’s static aeroelastic deflection using the baseline stiffness values. A cruise flight condition for the full-scale ESAC is considered to be at Mach = 0.797, altitude W 210 , 000 lbs h = 36 , 000 ft, with a wing loading of = corresponding to a design ¯ C = 0 . 510.

2 L S re f 1951 ft The code is first run restricting any coupled structural-aerodynamic loops. This solution thus represents the case where the structural deflections do not affect aerodynamics experienced on the flexible wing, or a model where the aerodynamics correspond to the rigid planform only. The deflection results are presented as W and Θ , where W tip tip tip is the wing tip vertical deflection (positive upwards), and Θ is the wing tip twist about the elastic axis (positive tip nose-down). These results are summarized in Table 1.

18 of 41 American Institute of Aeronautics and Astronautics ¯ C 0 . 510 L ¯ α , deg 2 . 279 W , ft 2 . 986 tip 2 W tip 100 ( ) , % 5 . 348 b Θ , deg − 0 . 351 tip 2 Θ deg tip − 3 ( ) , − 6 . 286 × 10 b ft Table 1. Structural Deflection Results for Full-Scale Wing Alone Model The model is then used to determine the aeroelastic deflection allowing the structural-aerodynamic loops. The results are summarized in Table 2.

¯ C 0 . 510 L ¯ α , deg 3 . 064 ¯ W , ft 2 . 740 tip 2 ¯ W tip 100 ( ) , % 4 . 907 b ¯ Θ , deg − 0 . 231 tip 2 ¯ Θ deg tip − 3 ( ) , − 4 . 135 × 10 b ft Table 2. Aeroelastic Deflection Results for Full-Scale Wing Alone Model

VIII. Wind Tunnel Model Scaling

A candidate wind tunnel model is generated from conducting scaling of the full-scale wing alone model. Scaling must be conducted considering several factors and desired characteristics.

A. Geometric Scaling Geometric scaling of the wing alone model is conducted so that the wind tunnel model can fit within the wind tunnel.

Given a desired wind tunnel model, a geometric scaling factor n can be determined based on the span of the full- scale scale wing alone model b and the desired span of the wind tunnel model b . The subscript w hereinafter refers to the f s sub-scale wind tunnel model characteristics.

Suppose a wind tunnel height is given to be 6 ft and a desired wind tunnel model semi-span is 5 . 4219 ft. A geometric scaling factor can be determined by: b 2 ( 56 . 1625 ft ) f n = = = 10 . 3585 (94) scale b 2 ( 5 . 4219 ft ) w The geometric scaling factor can be used to scale all the coordinates of the full-scale wing alone model that is used by the geometry generation tool. The reference values for the sub-scale wind tunnel model can also be obtained through the scaling factor.

S = S = 15 . 2919 ft (95) re f , w re f , f n scale ¯ c f ¯ c = = 1 . 6507 ft (96) w n scale The aspect ratio and the taper ratio remain unchanged.

AR = 7 . 6000 (97) w λ = 0 . 1950 (98) w The geometric scaling is not expected to affect the aerodynamics of the models given that all the reference values are correctly computed. To verify this, the lift and drag curves for the full-scale model are generated and compared to that of the scaled down wind tunnel model. Note that the drag polar only includes vortex-lattice computed drag.

19 of 41 American Institute of Aeronautics and Astronautics 1 0.12 0.9 0.1 0.8 0.7 0.08 0.6 L D 0.5 0.06 C C 0.4 0.04 0.3 0.2 0.02 0.1 Full − Scale Full − Scale Sub − Scale Sub − Scale 0 0 − 2 − 1 0 1 2 3 4 5 6 7 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 C α , angle of attack, deg L Figure 14. Lift and Drag Curve Verification for Full-Scale Wing Alone and Sub-Scale Wind Tunnel Model, Undeformed and Rigid The lift curve and drag polar for the geometrically scaled down wind tunnel model rests almost exactly on top of that of the full-scale wing alone model. This indicates that the scaling was done properly such that the total aerodynamics of the two models are preserved with geometric scaling.

B. Aeroelastic Scaling Aeroelastic scaling is conducted to determine the scaling factors on the wing’s stiffness properties. From the results in Tables 1 and 2, it can be seen that for a flexible wing model, structural-aerodynamic coupling changes the deformation result. It can even be modeled that aerodynamic considerations make the wing model stiffer in bending and softer in torsion while coupling the two motions, and this can be deduced from observing the terms in Eq. 81. Because the magnitudes of the deformations are generally similar when the total load over the wing is maintained, an uncoupled model will be used for preliminary sizing. This approach simplifies the scaling process and allows the torsional and bending stiffness of the wing to be analyzed separately due to the fact that the B and B terms in Eq. 18 are considered 1 2 to be negligible, removing any coupling between the two deformations.

1. Torsional Stiffness Scaling It is known that the structural deformation is calculated using the system equation given by Eq. 83. In the spirit of the finite-element analysis previously presented, if a single element modeled with constant structural properties is considered with a fixed end, the static deflection equation can approximated by − 1 θ = ( k ) f (99) s , torsion a , torsion GJ where the value of k is related to to the torsional constant J , k ∝ where l is the length of the beam s , torsion s , torsion l element.

In designing the wind tunnel model from the full-scale wing alone, it is desired that the relative elastic axis twist of the wing alone and the wind tunnel model are equal. That is, for each beam element: θ θ f w = (100) l l f w Let k = k and f = f , t s , torsion t a , torsion θ θ f w = (101) l l f w − 1 − 1 ( k ) f t t ( k ) f f f t t w w = (102) l l f w − 1 − 1 ( GJ / l ) f f f t ( GJ / l ) f f w w t w = (103) l l f w 20 of 41 American Institute of Aeronautics and Astronautics f t f f t w = (104) GJ GJ f w The force terms are expanded into the aerodynamic contributions, and S is the reference area for an element based e on Eq. 82.

e f = − ( c c + c cos Λ ) q cos Λ S (105) a , torsion m L ∞ e ac c The values of the aerodynamic coefficients for the wing alone and wind tunnel model as shown in Fig. 14 are equivalent: c = c = c and c = c = c .

m , f m , w m L , f L , w L ac ac ac ( ) 2 2 − c c + c e cos Λ q cos Λ S − ( c c + c e cos Λ ) q cos Λ S m f L f ∞ e , f m w L w ∞ e , w ac ac = (106) GJ GJ f w It is also known that the geometric parameters of the flying wing and the wind tunnel model are related through the value n thus allowing us to formulate a relationship between the torsional constants.

scale ( ) ( [ ] [ ] ) ( ) 2 2 2 − c c + c e cos Λ q cos Λ S − c c / n + c e / n cos Λ q cos Λ S / n m f L f ∞ e , f m f scale L f scale ∞ e , f ac ac scale = (107) GJ GJ f w ( ) ( ) 2 2 − c c + c e cos Λ q cos Λ S − c c + c e cos Λ q cos Λ S m f L f ∞ e , f m f L f ∞ e , f ac ac = (108) GJ GJ n f w scale GJ = GJ (109) w f n scale ︸ ︷︷ ︸ r torsion Thus the value r is defined as a scaling factor on the GJ of the flying wing model when scaling for the wind torsion tunnel model where r = .

torsion n scale The static structural deflection for the full-scale and the scaled down model with the application of the scaling factor r are tabulated in Table 3. It is shown that the application of the torsional scaling factor is able to match torsion the relative elastic axis twist relative to the span for both models when examined using an uncoupled structural- aerodynamic model.

¯ C = 0 . 510 Full Scale Wing Alone Wind Tunnel Model L ¯ α , deg 2 . 279 2 . 249 − 1 Θ , deg − 0 . 351 − 0 . 308 × 10 tip 2 Θ tip deg − 3 − 3 ( ) , − 6 . 286 × 10 − 5 . 707 × 10 b ft Table 3. Structural Deflection Results for Wind Tunnel Model with Torsional Stiffness Scaling In this analysis, the relative twist per semi-span was preserved in Eq. 100, and this is one of two ways in which torsional stiffness scaling can be conducted. Another approach for conducting the scaling would be instead to preserve the magnitude of the twist, or θ = θ (110) f w Following the same derivation process shown before, Eq. 104 becomes f l t f f l f t w w = (111) GJ GJ f w Substituting in Eq. 105, then ( ) 2 2 − c c + c e cos Λ q cos Λ S l − ( c c + c e cos Λ ) q cos Λ S l m f L f ∞ e , f f ac m w L w ∞ e , w w ac = (112) GJ GJ f w ( ) ( [ ] [ ] ) ( ) 2 2 2 − c c + c e cos Λ q cos Λ S l − c c / n + c e / n cos Λ q cos Λ S / n ( l / n ) m L ∞ m L ∞ f f e , f f f scale f scale e , f f scale ac ac scale = GJ GJ f w (113) 21 of 41 American Institute of Aeronautics and Astronautics ( ) ( ) 2 2 − c c + c e cos Λ q cos Λ S l − c c + c e cos Λ q cos Λ S l m f L f ∞ e , f f m f L f ∞ e , f f ac ac = (114) GJ GJ n f w scale If this scaling method is used, the result is GJ = GJ (115) w f n scale ︸ ︷︷ ︸ r torsion Thus two possible scaling factors for torsional deflection are developed. For scaling that preserves the relative twist relative to the span of the model, a scaling factor r = can be used. For scaling that preserves the magnitude torsion n scale of twist, a scaling factor r = can be used.

torsion n scale 2. Vertical Bending Stiffness Scaling Scaling the vertical bending stiffness EI is conducted similar to the torsional stiffness scaling. In the development yy for this wind tunnel model, however, the bending stiffness is desired to be scaled such that a 10% tip deflection relative to half of the span is achieved. For a beam element undergoing bending, the deflection can be given by − 1 w = ( k ) f (116) s , bending a , bending EI yy The value of k is given by the relationship k ∝ where l is the length of the beam element.

s , bending s , bending l Let k = k and f = f . Let an additional scaling factor n be defined to increase the tip deflec- b s , bending b a , bending bending tion to 10% relative to the model’s semi-span.

w w f w n = (117) bending l l f w − 1 − 1 ( k ) f b b ( k ) f f f b b w w n = (118) bending l l f w [ ] [ ] 3 − 1 3 − 1 ( EI ) / ( l ) f yy ( EI ) / ( l ) f f f b yy w w f b w n = (119) bending l l f w f ( l ) b f f ( l ) f b w w n = (120) bending ( EI ) ( EI ) yy f yy w Vertical bending from aerodynamic sources is due mainly from the sectional lift coefficient. Ignoring the contri- bution to the vertical bending force due to the component of pitching moment, the bending force on a beam element can be represented as f = c q cos Λ S (121) a , bending L ∞ e 2 2 2 2 c q cos Λ S ( l ) c q cos Λ S ( l ) L ∞ e , f f L ∞ e , w w n = (122) bending ( EI ) ( EI ) yy f yy w 2 2 2 2 c q cos Λ S ( l ) c q cos Λ S ( l ) L ∞ e , f f L ∞ e , f f n = (123) bending ( EI ) ( EI ) n yy f yy w scale ( EI ) = ( EI ) (124) yy w yy f n n bending scale ︸ ︷︷ ︸ r bending Thus the value r is defined as a scaling factor on the baseline EI when scaling for the wind tunnel model bending yy where r = .

bending n n bending scale Table 4 shows the results for the wind tunnel model’s static structural deflection when computed without any structural-aerodynamic coupling and with the application of the scaling factors derived. It can be seen that when the scaling factors are applied, a 10% wing tip deflection is achieved on the wind tunnel model.

22 of 41 American Institute of Aeronautics and Astronautics ¯ C = 0 . 510 Full Scale Wing Alone Wind Tunnel Model Wind Tunnel Model L ( n = 1 ) ( n = 0 . 537 ) bending bending ¯ α , deg 2 . 279 2 . 249 2 . 249 W , ft 2 . 986 0 . 289 0 . 539 tip 2 W tip 100 ( ) , % 5 . 348 5 . 370 10 . 000 b Table 4. Structural Deflection Results for Wind Tunnel Model with Vertical Bending Stiffness Scaling 3. Dynamic Pressure Effects A limitation in wind tunnel tests is that wind tunnel facilities may not be equipped to run at Mach numbers as high as the cruise condition of a full-scale aircraft. While use of a transonic wind tunnel can be used to test a model at high Mach number and high dynamic pressure q , operational and usage times at these facilities are very costly. For the ∞ VCCTEF study, it is not necessary that the wind tunnel model will need to run at high dynamic pressure provided that the elastic stiffness is scaled such that the lower dynamic pressure at the wind tunnel test speed is accounted for.

Let the value q represent the dynamic pressure of the wing alone model at a cruise Mach number and altitude.

∞ , f Let the value q represent the dynamic pressure of the wind tunnel model, which is restricted based on the wind ∞ , w tunnel test configuration.

If the wind tunnel dynamic pressure is the same as the wing alone model’s, the load on the wind tunnel model must be scaled to be the value W using the relationship: w W L w = = C q (125) L ∞ , f S S re f , w re f , w However, the wing loading on the wind tunnel model must be adjusted to take into account the change in dynamic ( ) q ∞ , w pressure to preserve the lift coefficient C . Multiplying Eq. 125 by the dynamic pressure ratio , a new value for L q ∞ , f ′ the load W for a design ¯ C can be determined.

L w ( ) ( ) ′ W W q q w ∞ , w ∞ , w w ¯ ¯ = = C q = C q (126) L L ∞ , w ∞ , f S S q q re f , w re f , w ∞ , f ∞ , f With respect to the torsional stiffness, the analysis beginning with Eq. 100 is still valid. However, the assumption that the dynamic pressure value q remains constant in Eq. 106 is no longer valid. Instead, the analysis needs to be ∞ adjusted as follows: ( ) 2 2 − c c + c e cos Λ q cos Λ S − ( c c + c e cos Λ ) q cos Λ S m f L f ∞ , f e , f m w L w ∞ , w e , w ac ac = (127) GJ GJ f w ( ) ( [ ] [ ] ) 2 2 2 − c c + c e cos Λ q cos Λ S − c c / n + c e / n cos Λ q cos Λ S / n m f L f ∞ , f e , f m f scale L f scale ∞ , w e , f ac ac scale = (128) GJ GJ f w ( ) ( ) ( ) 2 2 − c c + c e cos Λ q cos Λ S − c c + c e cos Λ q cos Λ S q m f L f ∞ , f e , f m f L f ∞ , f e , f ∞ , w ac ac = (129) GJ GJ n q f ∞ , f w scale ( ) 1 q ∞ , w GJ = GJ (130) w f q n ∞ , f scale The bending analysis is also altered from Eq. 122.

2 2 2 2 c q cos Λ S ( l ) c q cos Λ S ( l ) L ∞ , f e , f f L ∞ , w e , w w n = (131) bending ( EI ) ( EI ) yy f yy w ( ) 2 2 2 2 2 c q cos Λ S / n ( l / n ) c q cos Λ S ( l ) L ∞ , f e , f f L ∞ , w e , f f scale scale n = (132) bending ( EI ) ( EI ) yy f yy w ( ) 2 2 2 2 c cos Λ S ( l ) c cos Λ S ( l ) q L e , f f L e , f f ∞ , w n = (133) bending ( EI ) q ( EI ) n yy f yy w ∞ , f scale 23 of 41 American Institute of Aeronautics and Astronautics ( ) 1 q ∞ , w ( EI ) = ( EI ) (134) yy w yy f q n n bending ∞ , f scale It can be seen that in test situations where the wind tunnel dynamic pressure cannot be the same as the full-scale ( ) q ∞ , w flight condition, an additional factor equal to the dynamic pressure ratio must be added to the stiffness scaling q ∞ , f factors r and r such that torsion bending ( ) 1 q ∞ , w r = (135) torsion n q ∞ , f scale ( ) 1 q ∞ , w r = (136) bending q n n bending ∞ , f scale 4. Mach Number Effects The above analysis considers the case where the change in the dynamic pressure from q to q does not affect the ∞ , f ∞ , w aerodynamic coefficients. That is, c = c = c and c = c = c . This assumption allows for the clean m , f m , w m L , f L , w L ac ac ac derivation of the analyses above. In reality, the change in the dynamic pressure affects the Mach number of the flight condition, which affects the aerodynamic coefficients.

If the wind tunnel is operating at a different dynamic pressure q and at ambient sea level altitude, the Mach ∞ , w number of the wind tunnel model can be calculated.

( ) 1 2 q ∞ , w M = (137) w γ p SL This Mach number M is different than the Mach number M of the full scale wing alone model and is expected w f to be lower. This affects the aerodynamic coefficients.

The aerodynamic data for the wind tunnel model is obtained for two different flight conditions, assuming a rigid undeformed wing planform.

• The ESAC’s cruise condition at M = 0 . 797, h = 36 , 000 ft. For ¯ C = 0 . 510 the angle of attack is taken to be f L ◦ ¯ α = 2 . 249 .

lb f • An estimated wind tunnel flight condition of q = 20 , sea-level flight, corresponding to a Mach number ∞ , f ft ◦ M = 0 . 116. For ¯ C = 0 . 510 the angle of attack is taken to be ¯ α = 3 . 993 .

w L The lift curve of the wind tunnel model at both flight conditions are compared. The results show that the total aircraft C is different as a result of the Mach number effect, which is to be expected due to compressibility. A higher angle L α of attack for level flight is required to achieve the same ¯ C = 0 . 510 at the lower Mach number as well.

L 0.9 0.8 0.7 0.6 L 0.5 C 0.4 0.3 0.2 0.1 Mach = 0.797 Mach = 0.116 − 2 − 1 0 1 2 3 4 5 6 7 α , angle of attack, deg Figure 15. Mach Number Effect on Lift Curve 24 of 41 American Institute of Aeronautics and Astronautics The spanwise lift c and moment c distributions are shown in Fig. 16.

m l ac 0.7 − 0.01 Mach = 0.797 Mach = 0.797 − 0.02 Mach = 0.116 Mach = 0.116 0.65 − 0.03 0.6 − 0.04 − 0.05 0.55 l ac c − 0.06 m c 0.5 − 0.07 − 0.08 0.45 − 0.09 0.4 − 0.1 0.35 − 0.11 0 1 2 3 4 5 6 0 1 2 3 4 5 6 y, BBL, ft y, BBL, ft Figure 16. Mach Number Effect on Spanwise Aerodynamic Coefficients The spanwise lift distributions are similar in magnitude due to the fact that the overall load is maintained at ¯ C = L 0 . 510, but it is observed that the load shifts slightly outboard at the higher Mach number. The spanwise distribution of c is drastically different between the flight conditions. At the lower Mach number, the values of c are much m m ac ac lower. Thus, the aeroelastic scaling for the wind tunnel model will require consideration that Mach number effects can alter the spanwise c .

m ac The scaling analysis conducted for bending still holds because c does not affect the derivation of Eq. 136.

m ac However, additional analysis is required for aeroelastic scaling of torsion. The analysis is modified from Eq. 127.

While the simplification c = c = c is still made, c  = c .

L , w L m , w L , f m , f ac ac ( ) − c c + c e cos Λ q cos Λ S − ( c c + c e cos Λ ) q cos Λ S m , f f L f ∞ , f e , f ac m , w w L w ∞ , w e , w ac = (138) GJ GJ f w ( ) ( [ ] [ ] ) ( ) 2 2 2 − c c + c e cos Λ q cos Λ S − c c / n + c e / n cos Λ q cos Λ S / n m , f f L f ∞ , f e , f m , w f scale L f scale ∞ , w e , w ac ac scale = GJ GJ f w (139) ( ) ( ) ( ) 2 2 − c c + c e cos Λ q cos Λ S − c c + c e cos Λ q cos Λ S q m , f f L f ∞ , f e , f m , w f L f ∞ , f e , f ac ac ∞ , w = (140) GJ GJ n q f ∞ , f w scale ( ) ( ) ( ) 1 c q mr , w ∞ , w GJ = GJ q cos Λ S c (141) w f ∞ , f e , f f c q n mr , f ∞ , f scale where ( ) c e c e L f L f c = c + cos Λ = c + cos Λ (142) mr , f m , f m ac ac c c f f f and ( ) c e c e L f L f c = c + cos Λ = c + cos Λ (143) mr , w m , w m ac ac c c f f w c e L f For the flight condition, c is generally a negative number and is generally a positive number. However, the m ac c f Mach number effect introduces the possibility that c and c are opposing in sign due to the relative magnitudes mr , f mr , w ( ) c mr , w of the terms. This causes the ratio < 0. In effect, reducing the Mach number from M to M can cause the f w c mr , f wing aerodynamics to attempt to twist in an opposite direction. If this results, then GJ becomes a negative number, w which makes no physical sense. The interpretation is that there exists value of M or wind tunnel dynamic pressure w ( ) c mr , w q at which the ratio is negative and scaling the torsional stiffness of the wind tunnel model cannot achieve ∞ , w c mr , f the same amount of twist as the full scale flight condition.

25 of 41 American Institute of Aeronautics and Astronautics The values c and c are related to the torsional force experienced by a wing element about the elastic axis.

mr , f mr , w lb f For the ESAC’s cruise condition at M = 0 . 797 and the wind tunnel model’s Mach number for q = 20 , sea-level, f ∞ , f ft the values of c are compared.

mr 0.06 0.05 0.04 0.03 mr c 0.02 0.01 Mach = 0.797 Mach = 0.116 − 0.01 0 1 2 3 4 5 6 y, BBL, ft Figure 17. Mach Number Effect on Torsional Force c At C = 0 . 510 mr L It can be seen that c and c are the same in sign for most locations along the span. For the purpose of scaling mr , f mr , w torsional stiffness, this allows it to be possible for the wind tunnel model’s tip twist to be scaled to reflect the full-scale wing alone model’s at cruise.

For the purpose of illustration, the values of c for two other flight conditions are plotted in Fig. 18. The two mr flight conditions correspond to a ¯ C = 0 . 346. The full scale wing alone corresponds to M = 0 . 8 at a h = 30 , 000 ft L f lb f altitude. The wind tunnel flight condition corresponds to q = 20 at sea-level flight corresponding to M = 0 . 116.

∞ , f 2 w ft 0.015 0.01 0.005 − 0.005 mr c − 0.01 − 0.015 − 0.02 − 0.025 Mach = 0.8 Mach = 0.116 − 0.03 0 1 2 3 4 5 6 y, BBL, ft Figure 18. Mach Number Effect on Torsional Force c At C = 0 . 346 mr L It can be seen that c and c are opposite in sign for all locations along the span when scaling from M = 0 . 8.

mr , f mr , w f c mr , f This indicates that the ratio is generally negative along the span. The Mach number effects prevent the wind c mr , w 26 of 41 American Institute of Aeronautics and Astronautics tunnel model’s tip twist at the wind tunnel condition to be scaled to the cruise M = 0 . 8 condition.

5. Coupled Model Considerations The above analysis is intended for models that do not possess a coupled aerodynamic-structural nature and can be readily applied to problems where bending and torsion do not exhibit strong coupling. The actual static aeroelastic problem for the wind tunnel model possesses both coupling with an aerodynamic model, and the structural bending and torsion modes are also coupled together as a result of aerodynamic stiffening and softening.

Development of an analytical estimate for a scaling factor for the coupled model can be extremely intensive.

Instead, additional factors r and r representative of scaling factors are added into the relationships that should be 1 2 tuned heuristically using the static aeroelastic framework in Fig. 13 for the design ¯ C = 0 . 510. For stiff wing problems, L the coupled results and the uncoupled results are generally close and the scaling factors determined for the uncoupled model can be used. As even softer wing models are developed with higher vertical tip deflections, the usage of r and r factors become more necessary.

( ) 1 q ∞ , w GJ = r GJ = r GJ (144) w torsion f 1 f q n ∞ , f scale ︸ ︷︷ ︸ r torsion ( ) 1 q ∞ , w ( EI ) = r ( EI ) = r ( EI ) (145) yy w yy f 2 yy f bending n n q ∞ , f bending scale ︸ ︷︷ ︸ r bending The factor r is also used to represent a constant scaling factor to correct for mach number effects in Eq. 141.

2 ¯ W tip The relative percentage of aeroelastic wing tip deflection ( ) and the relative aeroelastic twist about the elastic axis b 2 ¯ Θ tip ( ) are determined using the static aeroelastic model of scaled wind tunnel model using the scaling factors r torsion b and r in Eqs. 144 and 145. The results are summarized in Table 5 and 6.

bending 2 ¯ W tip 100 ( ) , % r = 0 . 50 r = 0 . 80 r = 1 . 00 r = 1 . 50 2 2 2 2 b r = 1 . 00 14 . 474 10 . 032 8 . 343 5 . 871 r = 1 . 50 14 . 465 10 . 025 8 . 335 5 . 868 r = 1 . 75 14 . 466 10 . 020 8 . 336 5 . 868 r = 2 . 00 14 . 468 10 . 020 8 . 333 5 . 867 Table 5. Aeroelastic Relative Wing Tip Deflection for Wind Tunnel Model with r and r Scaling Factors 1 2 2 ¯ Θ deg tip ( ) , r = 0 . 50 r = 0 . 80 r = 1 . 00 r = 1 . 50 2 2 2 2 b ft − 3 − 3 − 3 − 3 r = 1 . 00 − 5 . 577 × 10 − 7 . 251 × 10 − 7 . 923 × 10 − 8 . 958 × 10 − 3 − 3 − 3 − 3 r = 1 . 50 − 3 . 717 × 10 − 4 . 828 × 10 − 5 . 264 × 10 − 5 . 965 × 10 − 3 − 3 − 3 − 3 r = 1 . 75 − 3 . 195 × 10 − 4 . 129 × 10 − 4 . 515 × 10 − 5 . 113 × 10 − 3 − 3 − 3 − 3 r = 2 . 00 − 2 . 789 × 10 − 3 . 613 × 10 − 3 . 946 × 10 − 4 . 472 × 10 Table 6. Aeroelastic Relative Wing Tip Elastic Axis Twist for Wind Tunnel Model with r and r Scaling Factors 1 2 6. Summary of Aeroelastic Scaling The scaling equations represented by Eqs. 144 and 145 are applied to the stiffness distributions of GJ and EI . While yy r and r can be defined as functions of span, this would result in a more complicated scaling procedure.

torsion bending Instead, the values of r and r are selected as single values applied to the entire baseline GJ and EI torsion bending yy distributions. This simplification can be done because the wind tunnel model’s aeroelastic behavior is scaled based on its tip deflection and tip twist.

27 of 41 American Institute of Aeronautics and Astronautics Given a value of n = 10 . 3585, aeroelastic scaling of the wind tunnel model from the cruise condition of Mach scale lb f number M = 0 . 797, h = 36 , 000 ft down to q = 20 ( M = 0 . 116), sea-level is conducted, where the design lift f ∞ , w w ft coefficient is ¯ C = 0 . 510. The scaling factors that were determined are summarized in Table 7.

L ¯ C 0 . 510 L n 0 . 537 bending r 1 . 75 r 0 . 80 q ∞ , w − 2 9 . 47 × 10 q ∞ , f − 4 r 1 . 28 × 10 bending − 6 r 4 . 42 × 10 torsion Table 7. Aeroelastic Scaling Factors For Wind Tunnel Model The application of the determined scaling factors results in aeroelastic deflections for the flexible wind tunnel model summarized in Table 8. The scaling factors are able to scale the wind tunnel model such that the wing tip ( ) 2 ¯ Θ tip deflection is 10 . 02% relative to the wind tunnel model’s semi-span, and the relative elastic axis twist is = b w ( ) 2 ¯ Θ − 3 deg tip − 3 deg − 4 . 129 × 10 , which compares to the full-scale model’s elastic axis twist of = − 4 . 135 × 10 .

ft b ft f ¯ C = 0 . 510 Full Scale Wing Alone Wind Tunnel Model L ¯ α , deg 3 . 064 5 . 773 ¯ W , ft 2 . 740 0 . 540 tip 2 ¯ W tip 100 ( ) , % 4 . 907 10 . 020 b − 1 ¯ Θ , deg − 0 . 231 − 0 . 223 × 10 tip 2 ¯ Θ tip deg − 3 − 3 ( ) , − 4 . 135 × 10 − 4 . 129 × 10 b ft Table 8. Aeroelastic Deflection Results for Wind Tunnel Model with Aeroelastic Stiffness Scaling C. Static Divergence An analysis of the scaled down wind tunnel model is conducted to determine the divergence dynamic pressure q , or d the dynamic pressure in which the wind tunnel model will experience static divergence. Determining the divergence dynamic pressure places a restriction on the wind tunnel test condition and is important to analyze to ensure that the model will be able to be properly utilized at the wind tunnel test conditions. If the divergence dynamic pressure is significantly larger than the test condition of the tunnel, static instability does not pose a problem. This, however, does not preclude the possibility of the wind tunnel model experiencing dynamic instability due to aeroelasticity, or flutter.

Flutter is not investigated within the scope of this study.

1. Torsional Divergence Initially, a preliminary analysis of static divergence can be performed on the wind tunnel model focusing only on torsional divergence. Torsional divergence is a classically examined phenomenon due to the basic aeroelastic coupling 7, 15 in twist. Aerodynamic forces can cause a wing to twist nose-up (negative Θ ). This nose-up twist increases the aeroelastic angle of attack α on the wing sections. Since lift force is proportional to aeroelastic angle of attack, this c can cause the wing to twist even more nose-up. Thus, a positive feedback loop exists between twist and angle of attack that can cause the wing to exceed its structural limitations and experience torsional divergence.

The structural stiffness and aerodynamic stiffness matrices for the globally assembled finite-element system were previous represented as K and K in Eq. 85 and assembled by the element matrices in Eqs. 80 and 81. Let the global a s stiffness matrices be partitioned as follows: [ ] K K s , t s , tb K = (146) s K K s , bt s , b 28 of 41 American Institute of Aeronautics and Astronautics [ ] K K a , t a , tb K = (147) a K K a , bt a , b where K and K are sub-matrices of the stiffness matrix elements corresponding to the torsional degrees-of-freedom s , t a , t [ ] , K and K are the sub-matrices of the stiffness matrix elements corresponding to the bend- θ θ . . . θ 1 2 n + 1 s , b a , b [ ] ′ ′ ′ ′ ing degrees-of-freedom , and K , K , K , and K are w w w w . . . w w . . . w w 1 2 i n + 1 s , tb s , bt a , tb a , bt 1 2 i n + 1 the coupling matrices.

A torsional divergence analysis involves examining the matrices K and K . Four cases are examined using s , t a , t different scaling factors r : torsion ( ) q q ∞ , w ∞ , w 1 − 2 • A first scaling case where r = r , r = 1 . 00, n = 10 . 3585, and = 9 . 47 × 10 .

torsion 1 1 scale q q n ∞ , f ∞ , f scale ( ) q q ∞ , w ∞ , w 1 − 2 • A second scaling case where r = r , r = 1 . 00, n = 10 . 3585, and = 9 . 47 × 10 .

torsion 1 1 scale q q n ∞ , f ∞ , f scale ( ) q q ∞ , w ∞ , w 1 − 2 • A third scaling case where r = r , r = 1 . 75, n = 10 . 3585, and = 9 . 47 × 10 .

torsion 1 1 scale q q n ∞ , f ∞ , f scale ( ) q q 1 ∞ , w ∞ , w − 2 • A fourth scaling case where r = r , r = 1 . 75, n = 10 . 3585, and = 9 . 47 × 10 .

torsion 1 3 1 scale q q n ∞ , f ∞ , f scale The aerodynamic stiffness matrix K is reliant on dynamic pressure as seen in Eq. 81, and static divergence occurs a , t when the term K + K is non-invertible or singular. In order to evaluate this, let the determinant of the total stiffness s , t a , t matrix normalized to determinant of the zero-speed structural stiffness matrix be defined as det ( K + K ) s , t a , t Δ = (148) t det ( K ) s , t lb f The wind tunnel test facility is limited at operating dynamic pressure q ≤ 60 . For this analysis, however, ∞ 2 ft lb f the value of Δ is plotted versus q ranging up to q = 400 for illustration purposes. The value in which Δ = 0 t ∞ ∞ t ft represents the divergence speed q .

d 1.2 0.8 0.6 − 4 r ∝ n , r =1.00 torsion scale 1 t Δ − 3 r ∝ n , r =1.00 torsion scale 1 0.4 − 4 r ∝ n , r =1.75 torsion scale 1 − 3 0.2 r ∝ n , r =1.75 torsion scale 1 − 0.2 0 50 100 150 200 250 300 350 400 q , dynamic presure, lbf/ft ∞ Figure 19. Δ versus Dynamic Pressure q , Torsional Divergence Analysis t ∞ lb f A torsional divergence dynamic pressure of q = 162 is observed for the first scaling case where r is the d torsion ft lb f smallest of the four scaling cases. The third scaling case has a torsional divergence dynamic pressure of q = 274 , d 2 ft while the second and fourth scaling cases where r ∝ have divergence dynamic pressures beyond q = torsion ∞ n xcale 29 of 41 American Institute of Aeronautics and Astronautics lb f 400 . It is clear that increasing the value of r , and thus the torsional stiffness of the model, causes the torsional torsion ft divergence dynamic pressure to increase due to the fact that the additional structural stiffness helps to prevent onset of structural instability. For all four scaling cases, the torsional divergence dynamic pressure is far beyond the desired lb f wind tunnel test speed of q = 20 , and thus, torsional divergence is not expected to be a problem for the test ∞ , w ft condition of the sub-scale wind tunnel model.

2. Coupled Static Divergence An analysis of torsional divergence represents an uncoupled analysis of the static instability problem for aeroelasticity.

The full FEM model developed in this study can actually be used to conduct a more refined analysis by examining the full K and K matrices in Eq. 85. Based on the sign of the terms in Eq. 81, aeroelasticity is expected to result in s a softening in torsion and stiffening in bending. The bending slope contributes to the aeroelastic angle attack through Eq. 46 and relieves the angle of attack, actually improving the divergence properties of the aeroelastic model and indicating that sole analysis of torsional divergence can actually be more conservative than the real system. Static divergence analyses for four scaling cases are examined: ( ) ( ) q q 1 ∞ , w 1 ∞ , w • A first scaling case where r = r , r = , r = 1 . 00, r = 1 . 00, n = torsion 1 4 bending 4 1 2 scale q q n ∞ , f n n ∞ , f bending scale scale q ∞ , w − 2 10 . 3585, and = 9 . 47 × 10 .

q ∞ , f ( ) ( ) q q 1 ∞ , w 1 ∞ , w • A second scaling case where r = r , r = , r = 1 . 00, r = 1 . 00, n = torsion 1 1 2 3 bending 4 scale q q n ∞ , f n n ∞ , f bending scale scale q ∞ , w − 2 10 . 3585, and = 9 . 47 × 10 .

q ∞ , f ( ) ( ) q q ∞ , w ∞ , w 1 1 • A third scaling case where r = r , r = , r = 1 . 75, r = 0 . 80, n = torsion 1 bending 1 2 scale 4 4 q q n n n ∞ , f ∞ , f bending scale scale q ∞ , w − 2 10 . 3585, and = 9 . 47 × 10 .

q ∞ , f ( ) ( ) q q 1 ∞ , w 1 ∞ , w • A fourth scaling case where r = r , r = , r = 1 . 75, r = 0 . 80, n = torsion 1 bending 1 2 scale 3 4 q q n ∞ , f n n ∞ , f bending scale scale q ∞ , w − 2 10 . 3585, and = 9 . 47 × 10 .

q ∞ , f Let the determinant of the total stiffness matrix normalized to the determinant of the zero-speed structural stiffness matrix be defined as det ( K + K ) s a Δ = (149) det ( K ) s The value of Δ is plotted versus q to determine the value q when K + K becomes singular.

∞ d s a − 4 r ∝ n , r =1.00, r =1.00 torsion scale 1 2 − 3 r ∝ n , r =1.00, r =1.00 torsion scale 1 2 − 4 r ∝ n , r =1.75, r =0.80 torsion scale 1 2 − 3 r ∝ n , r =1.75, r =0.80 torsion scale 1 2 50 Δ 0 50 100 150 200 250 300 350 400 q , dynamic presure, lbf/ft ∞ Figure 20. Δ versus Dynamic Pressure q , Coupled Static Divergence Analysis ∞ 30 of 41 American Institute of Aeronautics and Astronautics The results in Fig. 20 demonstrate the effect of modeling coupled bending-torsion on static divergence. The term Δ is positive for all the scaling cases within the dynamic pressure range examined. In fact, addition of the coupled bending-torsion consideration actually causes the term Δ to increase as dynamic pressure increases, indicating that static instability is not an issue for the model.

Let the total stiffness matrix be represented as K = K + K (150) s a This can be expanded using Eqs. 146 and 147 as [ ] [ ] [ ] K K K K K K t tb s , t s , tb a , t a , tb K = = + (151) K K K K K K bt b s , bt s , b a , bt a , b The terms K and K are negligible because B and B in Eq. 18 are considered to be zero. Thus Eq. 151 s , tb s , bt 1 2 becomes [ ] [ ] K K K + K K t s , t a , t tb a , tb K = = (152) K K K K + K bt b a , bt s , b a , b The matrix in Eq. 152 can be expressed as [ ] [ ] [ ] ( ) − 1 K + K K I K K + K − K K + K K 0 s , t a , t a , tb a , tb s , t a , t a , tb s , b a , b a , bt = ( ) (153) − 1 K K + K 0 K + K K + K K I a , bt s , b a , b s , b a , b s , b a , b a , bt which allows the determinant to be calculated as ( ) − 1 det ( K ) = det ( K + K ) det ( K + K − K K + K K ) s , b a , b s , t a , t a , tb s , b a , b a , bt − 1 − 1 = det ( K + K ) det ( K + K − K K K K ) (154) s , t a , t s , b a , b a , tb a , bt a , b s , b Because bending stiffness does not experience static instability, the determinant det ( K + K ) is positive. While s , b a , b it is difficult to make any generalizations about the sign of the second determinant term in Eq. 154, it can be concluded, however, that static divergence occurs only if − 1 − 1 det ( K + K − K K K K ) ≤ 0 (155) s , t a , t a , tb a , bt a , b s , b − 1 − 1 It is possible that the final term in Eq. 155 can be always positive, det ( K + K − K K K K ) > 0, based s , t a , t a , tb a , bt a , b s , b on the values of the elements in the stiffness matrices. This is the case for the aeroelastically scaled wind tunnel model whose results are in Fig. 20 where static divergence does not occur for the model. It can be seen that the torsional − 1 − 1 divergence problem is alleviated based on the − K K K K term in Eq. 155.

a , tb a , bt a , b s , b

IX. Wing Twist Optimization

To complete the development of the wind tunnel model, the geometrically and aeroelastically scaled model needs to be re-twisted for the wind tunnel test condition. A new unloaded shape is developed such that when the flexible wing model is operating at wind tunnel test condition, it aeroelastically deforms to a deflected shape that has minimum induced drag or maximum L/D ratio. This tailors the model such that it becomes ideal for conducting trade studies and drag analysis. An optimization procedure is developed and applied to the sub-scale wind tunnel model.

A. Optimization Method Optimization is achieved using an unconstrained gradient-based optimization algorithm. The foundation of a gradient- based optimization method is the determination of an optimal search direction which sufficiently minimizes the objec- tive function, calculated using the function gradient information. For this particular problem, the objective function is not an explicit analytical function, but instead the static aeroelastic mapping developed in Fig. 13. The input flight condition for the wind tunnel model is fixed, where ¯ C = 0 . 510, Mach number M = M = 0 . 116, and altitude h = 0 L w corresponds to sea-level testing conditions. However, a new design input is added which allows a user to add addi- tional twist onto the sub-scale model’s existing pre-twist distribution. Let the design input be expressed as a x which i contains η individual variables that specify an additional twist distribution on the sub-scale model.

31 of 41 American Institute of Aeronautics and Astronautics The models utilizes x to control the wing twist shape, and the static aeroelastic model determines the aeroelastic i ¯ ¯ shape of the ( ¯ Θ , W , W ) and the angle of attack ¯ α required for the static aeroelastic model to have a lift coefficient x equivalent to ¯ C . Aerodynamic modeling of the deformed shape also allows for the drag coefficient of the model C , L D to be determined for the static aeroelastic model. For all intents and purposes, then, the static aeroelastic model can be seen as an equivalent functional mapping such that C = J ( x ) (156) D c i where the objective function J is accomplished by utilizing the static aeroelastic mapping.

c Because J is not expressed analytically, the gradient of the objective function cannot be explicitly calculated and c [ ] η 1 2 must be approximated. Let x = . A forward finite-difference method is used to approximate the i x x . . . x i i i gradient about a known design point, x using the following: i , 0 [ ] T ∂ J ∂ J ∂ J c c c . . .

∇ J ( x ) = η (157) c i 1 2 ∂ x ∂ x ∂ x i i i ∂ J ( x ) J ( x + Δ x ) − J ( x ) c i , 0 c i , 0 i c i , 0 ≈ (158) ∂ x Δ x i , 0 i The gradient of a function points in the direction of greatest increase, so a possible search direction is one in the exact opposite direction of the gradient itself. This is called the direction of steepest descent. However, while choosing steepest descent will result in convergence on a minimum, it is known to be slow and inefficient. Therefore, several other methods have been developed to determine a more efficient search direction, such as the method of conjugate directions. The conjugate direction method uses the following property of conjugate vectors S [ H ] S = 0 , i  = j (159) i j where [ H ] is a symmetric and positive-definite matrix. For a quadratic function, where [ H ] is the Hessian of the function, if conjugate search directions S are used, then the conjugate property produces a complete decoupling which results in the ability to optimize the function in exactly η line searches, where η is the number of problem design variables. While this method is the most efficient for exactly quadratic problems, it is also effective for non-quadratic functions.

To start the conjugate direction method, the first search direction, S , is calculated using steepest descent since there is no prior gradient information. For following iterations, each updated search direction is found using information about the gradient at the current design point and the search direction from the previous iteration. The formula for updating the search direction at each iteration is k k k − 1 S = − ∇ J + β S (160) c where k T k k − 1 ( ∇ J ) ( ∇ J − ∇ J ) c c c β = (161) k − 1 k − 1 T ( ∇ J ) ( ∇ J ) c c using the Polak-Ribière method, one of several different possible β formulations.

It is possible for the search direction to become ill-conditioned due to numerical imprecision or if the objective function is particularly non-quadratic. Therefore, it is necessary to monitor whether or not the conjugate search direction does not result in a sufficient decrease in the objective function and reset to steepest descent if necessary.

There are two particular scenarios in which the search direction should be reset. The first situation occurs when the line search does not produce an improvement in the objective function. Second, each time a conjugate search direction is calculated, it should be compared with the gradient of the function at that design point. The closer the search direction is to the direction of the gradient, the less likely it is to result in sufficient decrease. In particular, the angle between ◦ the conjugate search direction vector and the negative of the gradient should be less than 90 .

With the search direction determined, a one-dimensional optimization in the search direction, or line search, is conducted. For the line search, the design variable is the search step in the particular search direction. The formal statement of the problem is as follows: k k min J ( x + σ S ) ; t > 0 (162) c i t 32 of 41 American Institute of Aeronautics and Astronautics k k where x are the actual problem design variables at iteration step k , S is the search direction at iteration step k , and σ i is the search step variable for the line search.

The line search consists of several steps in order to determine a minimum of the objective function in the search direction. First, basic bracketing is used to determine a lower and upper bound between which a minimum in the objective function exists. With upper and lower bounds defined, the bracketed interval is further refined using the Golden-Ratio Search method. This method uses the Golden-Section ratio value to reduce the interval around the minimum to a desired tolerance, in the fewest number of function evaluations.

For this method, two interior points in the interval are calculated using the following formulas: σ = σ + τ ( σ − σ ) (163) a lwr upr lwr σ = σ + ( 1 − τ )( σ − σ ) (164) b lwr upr lwr where τ is derived from the Golden-Section ratio, and is given by √ 3 − 5 τ = = 0 . 3819 (165) The objective function is then calculated at the two additional interior points, so that all four points, [ ] [ ] , and their corresponding function values, , are known. The σ σ σ σ J J J J lwr a b upr c , lwr c , a c , b c , upr interval refinement algorithm then proceeds as follows: • If J > J , then σ becomes the new upper bound, σ , and σ becomes the new σ interior point. A new c , b c , a b upr a b interior point, σ , is calculated using Eq. 163, along with it’s corresponding function evaluation, J .

a c , a • If J < J , then σ becomes the new lower bound, σ , and σ becomes the new σ interior point. A new c , b c , a a lwr b a interior point, σ , is calculated using Eq. 164, along with it’s corresponding function evaluation, J .

b c , b This process is continued until the interval has been reduced to a desired level of accuracy relative to the original interval.

Once the final interval bracketing the minimum has been calculated, the minimum of the objective function in the particular search direction can be approximated using a polynomial fit. In this case, since there are four points available from the Golden-Search method with four known function values, a cubic polynomial approximation of the objective function can be determined from which the approximate minimum can be calculated by finding the roots of the derivative of the polynomial approximation.

The result of the line search is the minimum of the objective function in the particular search direction and the ∗ corresponding minimum search point, σ , which is then used to find the next design point as follows: k + 1 k ∗ k x = x + σ S (166) i i This design point is then used as the initial design point for the next iteration of the optimization.

The method of calculating the search direction, performing a line search, and updating the search direction at the new design point is continued until a convergence criteria is met. For this problem, convergence is assumed when the absolute value of the objective function has not changed over several iterations.

B. Design Variable Distributions The design variables for this particular wing twist optimization problem control the values of the additional twist to be added to the already existing jig-shape twist of the flexible wing. Let the wing pre-twist ¯ γ ( x ) be the existing pre-twist on the jig-shape of sub-scale wind tunnel model, positive nose-down, about the pitch axis of the wind tunnel model.

Let the total wing pre-twist ˜ γ be represented as ˜ ¯ γ ( x ) = γ ( x ) + Δ γ ( x ) (167) where Δ γ ( x ) represents an additional pre-twist, positive nose-down, applied to the jig-shape on the wing about the pitch axis. The design variable x controls the distribution of Δ γ such that i Δ γ = f ( x ) (168) i Two different design variable cases are considered for the wing twist optimization: discrete point, and polynomial shape function coefficients.

33 of 41 American Institute of Aeronautics and Astronautics 1. Discrete-Point Design Variables For the first case, the design variables to be input into the optimization method are the actual additional twist values at two discrete points on the wing. In particular, the points are at the wing break point due to the wing trailing edge extension x , and the wing tip x . That is, for a discrete-point design variable optimization break tip [ ] x = (169) i Δ γ ( x ) Δ γ ( x ) break tip The additional twist is also constrained to be zero at the wing root Δ γ ( 0 ) = 0. The total twist distribution ¯ γ ( x ) is determined by linearly interpolating Δ γ ( x ) across the wing span and adding it to γ ( x ) .

2. Shape Function Design Variables For the second case, the additional twist along the wing span is represented by Chebyshev polynomial functions. In particular, the following shape function is initially considered: Δ γ ( x ) = a T ( x ) + a T ( x ) + a T ( x ) + a T ( x ) + a T ( x ) (170) 0 0 1 1 2 2 3 3 4 4 where T = 1 (171) T = x (172) T = 2 x − 1 (173) T = 4 x − 3 x (174) 4 2 T = 8 x − 8 x + 1 (175) However, in order to compare directly with the discrete optimization and prevent a under-constrained optimization problem, the functions needs to be modified such that the additional twist is always fixed to be zero at the wing root, Δ γ ( 0 ) = 0. This is done by subtracting the root value from Equation (170). For this particular model, the constant terms in the equation are eliminated.

Additional scaling of the shape function is done so that the polynomial coefficients stay nearby in order of magni- tude. The location along the wing x is scaled by the length of the wing L . Therefore, the final shape function that is used to describe the wing twist distribution for the model is given as ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) 2 3 4 2 x x x x x x x Δ γ = a + a 2 + a 4 − 3 + a 8 − 8 (176) 1 2 3 4 L L L L L L L The design variables for the shape function optimization are the four polynomial coefficients of the above shape function [ ] x = (177) a a a a i 1 2 3 4 C. Optimization Results The optimization framework is utilized to determine the new wind tunnel undeformed shape, such that when aeroelas- tically deformed at C = 0 . 510, C is minimized. This corresponds also to an optimized L/D for the wind tunnel test L D condition.

1. Discrete-Point Optimization A series of optimization runs are conducted to minimize C using wing pre-twist specified at discrete points along the D wing. The value of the additional pre-twist is fixed such that no additional wash-out is added to the root Δ γ ( 0 ) = 0, but the pre-twist at two locations on the wing are prescribed as design variables and linearly interpolated for stations in between. The two locations selected are the wing tip located at y = 5 . 388 ft, and y = 1 . 793 ft along pitch axis tip break or the y − axis in Fig. 12 or the b b b − direction in Fig. 8.

B 2 Without adding any additional pre-twist, the wind tunnel model has a C = 0 . 03228. Four optimization runs are D conducted using the discrete-point design input, each initialized at different starting values. The total amount of design variables for the discrete-point optimization is η = 2 corresponding to Δ γ at the wing extension break and the wing tip.

34 of 41 American Institute of Aeronautics and Astronautics The results of the value of the cost function J or C of the four optimization runs which were initialized at different c D starting values, are plotted in Fig. 21. For this case with only η = 2, the optimization requires only a minimal amount of iterations before converging to the minimum C value.

D 0.038 Run 1 Run 2 0.037 Run 3 Run 4 0.036 0.035 D 0.034 C 0.033 0.032 0.031 0.03 1 2 3 4 5 6 Iteration Number Figure 21. Cost Function ( C ) per Iteration of Discrete-Point Optimization D The optimization result corresponded to C = 0 . 03020, representing a 6 . 444% decrease in C relative to the un- D D optimized wind tunnel model. The additional pre-twist distribution that is added to the wind tunnel model planform results are summarized in Table 9, and a plot of the pre-twist distribution along the wing shown in Fig. 22. It can be seen that a nose-up additional pre-twist Δ γ is imposed on the wind tunnel model in order to reduce the C value at the D wind tunnel test condition.

y , BBL, ft Δ γ , deg, positive nose-down 0 0 1 . 793 − 4 . 220 5 . 388 − 6 . 203 Table 9. Optimization Result for Discrete-Point Optimization − 1 − 2 down − − 3 − 4 , deg, positive nose − 5 Δγ − 6 − 7 0 1 2 3 4 5 6 y, BBL, ft Figure 22. Optimized Δ γ Result for Discrete-Point Optimization 35 of 41 American Institute of Aeronautics and Astronautics Deflection information for the static aeroelastic model for the un-optimized and the discrete-point optimized results are summarized in Table 10.

¯ C = 0 . 510 Un-optimized Wind Tunnel Model Optimized Wind Tunnel Model L ¯ α , deg 5 . 773 2 . 411 C , counts 322 . 8 302 . 0 D ¯ W , ft 0 . 540 0 . 649 tip 2 ¯ W tip 100 ( ) , % 10 . 020 11 . 872 b − 1 − 1 ¯ Θ , deg − 0 . 223 × 10 − 0 . 343 × 10 tip 2 ¯ Θ deg tip − 3 − 3 ( ) , − 4 . 129 × 10 − 6 . 371 × 10 b ft Table 10. Aeroelastic Deflection Results for Discrete-Point Optimized Wind Tunnel Model The lift distribution of the un-optimized and discrete-point optimized wind tunnel models undergoing aeroelastic deformation are also shown in Fig. 23.

1.4 Un − optimized Discrete − Point Optimized 1.2 0.8 cl*c 0.6 0.4 0.2 0 1 2 3 4 5 6 y, BBL, ft Figure 23. Lift Distribution of Wind Tunnel Model Using Discrete-Point Optimization Results The un-optimized model’s lift distribution is generally triangular in shape, due in part to the aeroelastic deformation which causes the wing tip to twist nose-down. This nose-down aeroelastic deformation causes the lift distribution to shift towards the wing root. The addition of the optimized Δ γ corrects this, however, and a new pre-twist is prescribed that twists the wing tip more nose-up as shown in Fig. 22. The resulting lift distribution becomes more elliptical in shape.

2. Shape Function Optimization A second series of optimization runs are conducted using the shape function in Eq. 176 to prescribe additional pre- twist on the wing Δ γ . In this case, the design input variable x has increased to η = 4 degrees-of-freedom, where the i input design variable represented in Eq. 177, corresponds to the shape function coefficients in Eq. 176. The shape function imposes no additional wash-out added to the root, or Δ γ ( 0 ) = 0.

A total of seven optimization runs are conducted from different initial values, and the evolution and decrease in the value of C as iterations are conducted is plotted in Fig. 24. While it requires more iterations before the cost D function decreases to a minimum value, Fig. 24 shows that the optimization is able to drive C of the wind tunnel D model at the test conditions down. The independent optimization runs also generally converge to similar minimum C values. The lowest minimum C value obtained in the optimization study was C = 0 . 03018 corresponding to a D D D 6 . 506% improvement from the un-optimized wind tunnel model shape. The shape function optimization C result is D lower than that of the discrete-point optimization, but is to be expected due to the fact that Δ γ is parametrized by more degrees-of-freedom in the shape function optimization.

36 of 41 American Institute of Aeronautics and Astronautics 0.055 Run 1 Run 2 Run 3 Run 4 0.05 Run 5 Run 6 Run 7 0.045 D C 0.04 0.035 0.03 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 Iteration Number Figure 24. Cost Function ( C ) per Iteration of Shape Function Optimization D The optimization result is summarized in Table 11 for the parameters of the shape function, where it is seen that normalization of the independent variable in Eq. 176 is able to return coefficient results of similar order.

Parameter Value a − 0 . 3379 a − 0 . 9700 a 0 . 5701 a 1 . 3101 Table 11. Optimization Result for Shape Function Optimization The optimization results translate into an additional pre-twist distribution Δ γ ( x ) added to the wind tunnel model planform, and the distribution is shown in Fig. 25. By using polynomial basis functions, the pre-twist distribution is a smooth continuous curve.

− 1 − 2 down − − 3 − 4 − 5 , deg, positive nose Δγ − 6 − 7 − 8 0 1 2 3 4 5 6 y, BBL, ft Figure 25. Optimized Δ γ Result for Shape Function Optimization The static aeroelastic model is used to obtain the deflection results for the un-optimized and the shape function optimized results, and the values are summarized in Table 12.

37 of 41 American Institute of Aeronautics and Astronautics ¯ C = 0 . 510 Un-optimized Wind Tunnel Model Optimized Wind Tunnel Model L ¯ α , deg 5 . 773 2 . 345 C , counts 322 . 8 301 . 8 D ¯ W , ft 0 . 540 0 . 644 tip 2 ¯ W tip 100 ( ) , % 10 . 020 11 . 949 b − 1 − 1 ¯ Θ , deg − 0 . 223 × 10 − 0 . 355 × 10 tip 2 ¯ Θ tip deg − 3 − 3 ( ) , − 4 . 129 × 10 − 6 . 594 × 10 b ft Table 12. Aeroelastic Deflection Results for Shape Function Optimized Wind Tunnel Mode The lift distribution of the un-optimized and shape function optimized wind tunnel models are also shown in Fig.

26. The lift distribution of the wind tunnel model with wing pre-twist optimized using shape functions is very similar to the lift distribution of the model with wing pre-twist optimized using the discrete-points. Both optimization results impose a nose-up twist onto the wing to counteract the nose-down twist due to aeroelastic deformation.

1.4 Un − optimized Shape Function Optimized 1.2 0.8 cl*c 0.6 0.4 0.2 0 1 2 3 4 5 6 y, BBL, ft Figure 26. Lift Distribution of Wind Tunnel Model Using Shape Function Optimization Results 3. Summary The results of the discrete-point optimization and the shape function optimization are used to analyze the static aeroe- lastic model of the wind tunnel model. Lift curves and drag polars are plotted in Fig. 27 representing the flexible wind tunnel model, and the curves provide insight into the optimization results of the model.

Both the discrete-point optimization result and the shape function optimization result produce lift curves which are very similar to each other, and the lift curves for the optimized models are shifted upwards of the un-optimized lift curve. This means that the additional optimized pre-twist Δ γ helps to recover the loss of lift due to the nose-down aeroelastic deformation of the flexible wind tunnel model.

The drag polar of the flexible sub-scale wind tunnel model is also affected when the optimized pre-twist Δ γ is applied to the static aeroelastic model. Though slight, the drag polars of the optimized models are shifted, and the drag polar at the design ¯ C value is lower than that of the unoptimized model. To recall, the discrete-point optimized model L observed a 6 . 444% decrease in C at the test condition relative to the un-optimized model and the shape function D optimized model observed a 6 . 506% decrease.

38 of 41 American Institute of Aeronautics and Astronautics 1 0.16 Un − optimized Un − optimized Discrete Point Optimized Discrete Point Optimized 0.9 0.14 Shape Function Optimized Shape Function Optimized 0.8 0.12 0.7 0.1 0.6 L D 0.5 0.08 C C 0.4 0.06 0.3 0.04 0.2 0.02 0.1 0 0 − 2 0 2 4 6 8 10 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 C α , angle of attack, deg L Figure 27. Lift Curves and Drag Polars for Aeroelastic Wind Tunnel Models The aeroelastic deflections for the optimized and un-optimized aeroelastic models are plotted in Fig. 28.

1.4 0.06 Un − optimized Un − optimized Discrete Point Optimized Discrete Point Optimized 0.04 Shape Function Optimized Shape Function Optimized 1.2 down − 0.02 0.8 − 0.02 0.6 − 0.04 0.4 − 0.06 , tip elastic axis twist, positive nose W, vertical wing tip deflection, positive up Θ 0.2 − 0.08 0 − 0.1 − 2 0 2 4 6 8 10 − 2 0 2 4 6 8 10 α , angle of attack, deg α , angle of attack, deg Figure 28. Aeroelastic Deflections for Wind Tunnel Models Because of the nose-up pre-twist applied that resulted in the shift of the lift curves of the optimized models, the vertical bending deflection W for the optimized models is higher. The higher lift also contributes to the moment tip about the elastic axis of the wing, thus driving Θ more negative, or nose-up.

tip The aircraft rigid body angle of attack is α , the aeroelastic deformation effect on the angle of attack is α , and Δ γ e represents the additional prescribed wash-out determined through optimization–all about the aircraft pitch axis where α and α are positive nose-up, and Δ γ is positive nose-down. Let a new quantity α be defined such that e p α = α + α − Δ γ (178) p e where α represents the angle of attack of a local section of the wing perpendicular to the pitch axis, positive nose-up, p relative to the wind tunnel un-optimized jig-shape existing pre-twist ¯ γ . Thus, α represents the effective angle of attack p of a wing section relative to the un-optimized existing pre-twist ¯ γ . The value of α is plotted versus the aircraft p , tip angle of attack α for the aeroelastic wind tunnel models in Fig. 29.

39 of 41 American Institute of Aeronautics and Astronautics up − Un − optimized − 2 Discrete Point Optimized , physical angle of attack, deg, positive nose p α Shape Function Optimized − 4 − 2 0 2 4 6 8 10 α , angle of attack, deg Figure 29. Physical Angle of Attack at Wing Tip for Aeroelastic Wind Tunnel Models The plot of α shows the physical angle of attack of the wing tip relative to the jig-shape for the un-optimized p , tip and optimized models, where it is the largest for the model optimized using the shape function. For an un-optimized model with a rigid planform, α = α , and thus the plot of α for the un-optimized model demonstrates the effect of p p , tip aeroelastic deformation only. Each plot of the optimized models, thus, represents the effect of re-twisting of the wing and the aeroelastic deformation due to the re-twisting.

X. Conclusion

This study presents the development of a static aeroelastic model of a flexible sub-scale wind tunnel model based upon scaling of a full-scale transport aircraft wing of the NASA GTM. A static aeroelastic framework is developed by coupling a structural model of the flexible wing and the aerodynamic model of the aircraft. The structural model is constructed using finite-element modeling of an equivalent one-dimensional simple beam model of the wing. Aero- dynamic modeling is conducted using a vortex-lattice solution. The resulting static aeroelastic model is a coupled finite-element vortex-lattice model capable of converging aeroelastic solutions for the flexible wing model and devel- oping flexible aircraft lift curves and drag polars.

The static aeroelastic model is implemented on a full-scale wing model of the ESAC or GTM. In order to develop the sub-scale model, a scaling procedure is developed first by geometrically scaling the full-scale model to sub-scale, then by conducting aeroelastic scaling. Aeroelastic scaling is conducted to scale the torsional stiffness of the sub-scale model such that the relative twist to span ratio matches that of the full-scale model, and the bending stiffness is scaled such that the sub-scale model has 10% wing tip deflection. Additional aeroelastic scaling is conducted to take into account the lower dynamic pressure and mach number of the wind tunnel test relative to the full-scale model’s cruise flight condition. Heuristic scaling factors are added to the scaling to account for coupled bending-torsion effects due to aerodynamic stiffening/softening. A static divergence analysis is conducted on the final, fully scaled flexible wind tunnel model to evaluate the risk of static instability of the model in wind tunnel testing.

A final design of the wind tunnel model is developed by re-twisting the model through optimization targeted at minimizing induced drag or maximizing L/D at the wind tunnel design test condition. The gradient-based optimization approach is based on utilizing one-dimensional line searches in conjugate search directions. Two optimizations are conducted: one where additional pre-twist is applied to the wind tunnel model by linearly interpolating between pre- twist values specified at discrete points along the wing corresponding to the wing trailing edge extensions break and tip, and one where additional pre-twist is applied based on a shape function inspired by Chebyshev polynomials. The resulting optimized pre-twist of the wing is able to reduce the drag coefficient at the design test condition by 6 . 444% when optimized by specifying twist at discrete points, and 6 . 506% when optimized by specifying twist through a shape function.

The final result of the study is a flexible sub-scale wind tunnel model configuration with static aeroelastic similarity to the full-scale ESAC wing but with increased wing tip deflection and tailored for the design test condition. An aeroelastic model accompanies the developed configuration that can be used in future validation against wind tunnel 40 of 41 American Institute of Aeronautics and Astronautics testing. A clean wing model is analyzed, but investigation of the VCCTEF control surface can further use the model developed in future studies.

XI. Acknowledgments

The authors would like to thank the NASA Aeronautics Research Mission Directorate (ARMD) Fixed Wing Project under the Fundamental Aeronautics Program for providing the funding to support this work.

References

Nguyen, N., “NASA Innovation Fund 2010 Project: Elastically Shaped Future Air Vehicle Concept,” NASA Internal Report Submitted to NASA Innovative Partnerships Program Office, October 8, 2010.

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

Urnes, Sr., J., Nguyen, N., Ippolito, C., Totah, J., Trinh, K., 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,” 51st 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 Testbest for Flight Research Experiments,” AUVSI Unmanned Unlimited, Arlington, VA, 2004.

Nguyen, N., Trinh, K., Frost S., and Reynolds, K., “Coupled Aeroelastic Vortex Lattice Modeling of Flexible Aircraft,” AIAA Applied Aerodynamics Conference, AIAA-2011-3021, June 2011.

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

Hodges, D.H. and Pierce, G.A., Introduction to Structural Dynamics and Aeroelasticity , Cambridge University Press, 2002.

Craig, Jr., R.R. and Kurdila, A.J., Fundamentals of Structural Dynamics , Second Edition, John Wiley & Sons, Inc., 2006.

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

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.

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.

Miranda, L.R., Elliot, R.D., and Baker, W.M., “A Generalized Vortex Lattice Method for Subsonic and Supersonic Flow Applications,” 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.

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

Nguyen, N., Ting, E., Swei, S., Ishihara, A., “Distributed Parameter Optimal Control by Adjoint Aeroelastic Differential Operators for Mode Suppression Control,” AIAA Guidance, Navigation, and Control Conference, AIAA-2013-4859, August 2013.

Vanderplaats, G. N., Multidiscipline Design Optimization , Vanderplaats Research & Development, Inc., Monterey, CA, 2007.

Fletcher, R. and Reeves, C. M., “Function Minimization by Conjugate Gradients,” The Computer Journal , Vol. 7, No. 2, 1964, pp. 149-154.

Polak, E., and Ribiere, G., “Note sur la convergence de méthodes de directions conjuguées,” Revue française d’informatique et de recherche opérationnelle , Vol. 3, No. 1, 1969, pp. 35-43.

41 of 41 American Institute of Aeronautics and Astronautics