Aeroelasticity of Axially Loaded Aerodynamic Structures for Truss-Braced Wing Aircraft

Aeroelasticity of Axially Loaded Aerodynamic Structures for

Truss-Braced Wing Aircraft

∗ Nhan Nguyen NASA Ames Research Center, Moffett Field, CA 94035 † Eric Ting Stinger Ghaffarian Technologies, Inc., Moffett Field, CA 94035 ‡ Sonia Lebofsky Stinger Ghaffarian Technologies, Inc., Moffett Field, CA 94035 This paper presents an aeroelastic finite-element formulation for axially loaded aerodynamic structures.

The presence of axial loading causes the bending and torsional sitffnesses to change. For aircraft with axially loaded structures such as the truss-braced wing aircraft, the aeroelastic behaviors of such structures are non- linear and depend on the aerodynamic loading exerted on these structures. Under axial strain, a tensile force is created which can influence the stiffness of the overall aircraft structure. This tension stiffening is a geometric nonlinear effect that needs to be captured in aeroelastic analyses to better understand the behaviors of these types of aircraft structures. A frequency analysis of a rotating blade structure is performed to demonstrate the analytical method. A flutter analysis of a truss-braced wing aircraft is performed to analyze the effect of geo- metric nonlinear effect of tension stiffening on the flutter speed. The results show that the geometric nonlinear tension stiffening effect can have a significant impact on the flutter speed prediction. In general, increased wing loading results in an increase in the flutter speed. The study illustrates the importance of accounting for the geometric nonlinear tension stiffening effect in analyzing the truss-braced wing aircraft.

I. Introduction Future aircraft design trends are employing high-aspect ratio wing technology to improve aerodynamic efficiency for reduced fuel burn. As a wing aspect ratio increases, the need for maintaining sufficient structural load carrying capacity is becoming more important. Cantilever wing design can only accommodate up to a certain aspect ratio beyond which the wing root bending moment can become too large that will impose structural limitations on a wing design. At the same time, wing deflection can also become too excessive that could cause aerodynamic penalty.

Truss-braced wing aircraft concepts provide a structural solution to high aspect ratio wing aircraft designs. The long slender wing would employ structural bracing via the use of axially loaded strut members to provide intermediate span supports in addition to the wing root attachment. These struts generally support a portion of the spanload carried by the wing and are loaded in tension. Under a negative-g flight condition such as during a dive, a load reversal could occur ∗ NASA Ames Research Center, Research Scientist, AIAA Associate Fellow, nhan.t.nguyen@nasa.gov † Stinger Ghaffarian Technologies Inc., NASA Ames Research Center, Research Engineer, eric.b.ting@nasa.gov ‡ Stinger Ghaffarian Technologies Inc., NASA Ames Research Center, Research Engineer, eric.b.ting@nasa.gov 1 of 26 American Institute of Aeronautics and Astronautics that could put the struts in compression. The compressive loading would require design considerations for buckling strength.

Under aerodynamic loading, an axially loaded member also experiences the normal bending and torsion generated by aerodynamic lift force and pitching moment. Aeroelasticity of an axially loaded structure undergoing transverse bending can be significantly different from that with transverse bending alone.

The SUGAR Truss-Braced Wing (TBW) aircraft concept is a Boeing-developed N+3 aircraft configuration funded 1, 2 by NASA ARMD Fixed Wing project. The TBW aircraft concept is designed to be aerodynamically efficient by employing a high-aspect ratio wing design. The aspect ratio of the TBW is in the order of 14 which is significantly greater than those of conventional aircraft wings. As a result, intermediate structural supports are required. The main wings are braced at approximated at mid-span by two main struts. In addition, two jury struts; one on each wing, provide additional reinforcement. Figure 1 is an illustration of the TBW aircraft.

Fig. 1 - Boeing SUGAR Truss-Braced Wing (TBW) Aircraft Concept The main struts and jury struts can impact aerodynamics of the main wings since they create aerodynamic inter- ference drag that can offset the aerodynamic benefit of high-aspect ratio wing design. Thus, the TBW is designed to strive a balance between aerodynamic efficiency and structural efficiency.

In this paper, we will develop an aeroelastic analysis for axially loaded aerodynamic structures such as the TBW main struts. The presence of axial loading causes the bending and torsional sitffnesses to change. A tensile loading will result in an increase in bending and torsional stiffnesses. This tension stiffening is a geometric nonlinear effect which should be captured in aeroelastic analyses of aircraft structures that carry significant axial loads. A finite- element formulation is developed to model the axially loaded structures in the presence of the bending-torsion motion of a typical aircraft wing. A flutter analysis is performed to analyze the geometric nonlinear effect of tension stiffening on the flutter speed.

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

Consider an airfoil section on the left wing as shown in Fig. 2 undergoing bending and torsional deflections. Let ( x , y , z ) be the undeformed coordinates of point Q on a wing airfoil section in the reference frame D. Let p = x d be 0 1 a position vector along the elastic axis. Then, point Q is defined by a position vector p = p + q where q = y d + z d 0 2 3 defines point Q in the y − z plane from the elastic axis. Then the undeformed local airfoil coordinates of point Q are       y cos γ − sin γ η       = (1) z sin γ cos γ ξ where η and ξ are local airfoil coordinates, and γ is the wing section pre-twist angle, positive nose-down.

Differentiating y and z with respect to x gives         ′ y ′ − sin γ − cos γ η − z γ x         = γ = (2) ′ z cos γ − sin γ ξ y γ x Fig. 2 - Left Wing Reference Frame of Wing in Combined Bending-Torsion Let Θ be a torsional twist angle about the x -axis, positive nose-down. Let W and V be flapwise and chordwise bending deflections of point Q, respectively. Let U be the axial displacement of point Q. Then, the displacement and rotation vectors due to the elastic deformation can be expressed as r = U d + V d + W d (3) 1 2 3 3 of 26 American Institute of Aeronautics and Astronautics φ = Θ d − W d + V d (4) 1 x 2 x 3 where the subscripts x and t denote the partial derivatives of Θ , W , and V .

Let ( x , y , z ) be the deformed coordinates of point Q on the airfoil in the reference frame D and p = x d + 1 1 1 1 1 1 y d + z d be its position vector. Then the coordinates ( x , y , z ) are computed as 1 2 1 3 1 1 1 p = p + r + φ × q (5) where     x x + U − yV − zW 1 x x         = (6) y y + V − z Θ     z z + W + y Θ Differentiating x , y , and z with respect to x yields 1 1 1     ′ ′ x 1 + U − yV + z γ V − zW − y γ W 1 , x x xx x xx x     ′ ′     = (7) y − z γ + V − z Θ − y γ Θ 1 , x x x     ′ ′ z y γ + W + y Θ − z γ Θ 1 , x x x Neglecting the transverse shear effect, the longitudinal strain is computed as ds − ds s 1 1 , x ε = = − 1 (8) ds s x where √ √ ( ) ′ 2 2 2 2 2 s = 1 + y + z = 1 + ( y + z ) γ (9) 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 U − 2 yV − 2 zW + 2 ( y + z ) γ Θ + ( x − 1 ) + y + z γ + z − y γ (10) x xx xx x 1 , x 1 , x 1 , x x Using the Taylor series expansion, s is approximated as 1 , x ( ) ( ) 2 2 ′ ′ ( ) ′ 2 2 ( x − 1 ) + y + z γ + z − y γ 1 , x 1 , x 1 , x U − yV − zW + y + z γ Θ x xx xx x s ≈ s + + (11) 1 , x x s 2 s x x The longitudinal strain is then obtained as ( ) ( ) 2 2 ′ ′ ( ) ′ 2 2 ( x − 1 ) + y + z γ + z − y γ 1 , x 1 , x 1 , x U − yV − zW + y + z γ Θ x xx xx x ε = + (12) 2 2 s 2 s x x ( ) ′ For a small wing twist angle γ , γ ≈ 0. Then ( ) ( ) 2 2 ′ ′ ( x − 1 ) + y + z γ + z − y γ 1 , x 1 , x 1 , x ( ) ′ 2 2 ε = U − yV − zW + y + z γ Θ + (13) x xx xx x 4 of 26 American Institute of Aeronautics and Astronautics The axial force and moments acting on a wing are then obtained as     ∫∫ E ε dydz P x ( )     ∫∫ ( ) ′    2 2  GJ Θ + E ε y + z γ + Θ dydz M x x  x      =   ∫∫     M y − E ε zdydz     ∫∫ M z − E ε ydydz       ′ EA EI γ − EAe − EAe xx z y U ∆ P [ ] x x   ( ) 2     ′ ′ ′ ′        EI γ GJ + EB γ + EI U − EB γ − EB γ  xx 1 xx x 2 3 Θ ∆ M  x   x    = + (14)       ′       W ∆ M xx y     − EAe − EB γ EI − EI z 2 yy yz   ′ V ∆ M xx z − EAe − EB γ − EI EI y 3 yz zz where ∆ P and ∆ M are the nonlinear terms; E is the Young’s modulus; G is the shear modulus; A is the tensile x ( x , y , z ) area; e and e are the offsets of the centroid from the elastic axis; I , I , and I are the section area moments of inertia y z yy yz zz about the flapwise axis; J is the torsional constant; and B , B , and B are the bending-torsion coupling constants which 1 2 3 are defined as     2 2 B y + z ∫ ∫     ( ) 2 2     = y + z dydz (15) B z     B y ′ The slope of the twist angle γ can play a significant role in structures with large twists such as turbomachinery blades. For aircraft wings, this effect is negligible and therefore can be neglected. Also, in general, the offset of the centroid from the elastic axis in the z -direction, e , is small compared to that in the y − direction e . Therefore, it may z y be neglected. Then, the moments acting a wing are reduced to       P EAU − EAe V ∆ P x x y xx x             M ( GJ + EI U ) Θ ∆ M  x   xx x x   x    =   +   (16)       M EI W − EI V ∆ M y yy xx yz xx y       M − EAe U − EI W + EI V ∆ M z y x yz xx zz xx z Note that the signs of the moments are defined in the positive deflection sense such that M = M d − M d + M d (17) x 1 y 2 z 3 5 of 26 American Institute of Aeronautics and Astronautics The nonlinear terms are given by     ∆ P 1 x     ( ) ∫ ∫     [ ] 2 2 ∆ M 1 y + z Θ  x   x  2 2 2   = E ( U − yV − zW ) + ( V − z Θ ) + ( W + y Θ )   dydz x xx xx x x x x     ∆ M − z y     ∆ M − y z         EA EI EI EI zz yy xx                 ( ) 1 EI Θ 1 EB Θ 1 EB Θ 1 EB Θ  xx x   4 x   5 x   1 x  2 2 2 2 2 2 = U + V + W   + V   + W   + Θ   x x x xx xx x 2   2   2   2   − EAe − EB − EB − EB z 7 8 2         − EAe − EB − EB − EB y 10 9 3         − EAe − EAe − EI − EAe y z yz z                 − EB Θ − EB Θ EB Θ − EB Θ  3 x   2 x   6 x   2 x  + U V + U W + V W + V Θ         x xx x xx xx xx x x         − EI EI − EB EI yz yy 9 yy         EI − EI − EB − EI zz yz 7 yz   EAe y     EB Θ  3 x  + W Θ (18)   x x   EI yz   − EI zz where     ( ) 2 2 2 B y y + z     ( )     2 2 2  B   z y + z      ( )     2 2 B yz y + z  6    ∫ ∫        2  = dydz (19) B y z             B z             B yz     B y Consider two cases: 1. If the nonlinear strain component due to axial displacement ( x − 1 ) ≈ 0, then the nonlinear terms are simpli- 1 , x fied to           ∆ P EA EI − EAe EAe x xx z y                     ( ) ∆ M EI Θ EB Θ − EB Θ EB Θ 1 1  x   xx x   1 x   2 x   3 x  2 2 2 = V + W + Θ + V Θ + W Θ (20)           x x x x x x x   2   2       ∆ M − EAe − EB EI EI y z 2 yy yz           ∆ M − EAe − EB − EI − EI z y 3 yz zz ( ) ′ 2. If, in addition, the nonlinear strain component due to chordwise deflection y + z γ ≈ 0, then the nonlinear 1 , x 6 of 26 American Institute of Aeronautics and Astronautics terms are simplified to         ∆ P EA EI EAe x yy y                 ∆ M 1 EI Θ 1 EB Θ EB Θ  x   xx x   4 x   3 x  2 2   = W   + Θ   + W Θ   (21) x x x x         2 2 ∆ M − EAe − EB EI y z 7 yz         ∆ M − EAe − EB − EI z y 10 zz III. Aeroelastic Angle of Attack The relative velocity of the air approaching a wing section includes the contribution from the wing elastic deflection that results in changes in the local angle of attack. Since aerodynamic forces and moments are dependent on the local angle of attack, the wing aeroelastic deflections will generate additional elastic forces and moments. The local angle of attack depends on the relative approaching air velocity as well as the rotation angle φ from Eq. (4). The relative air velocity in turn also depends on the deflection-induced velocity. The velocity at point Q due to the aircraft velocity and angular velocity in the reference frame D is computed as v = ¯ v + ω × r = ( u b + v b + w b ) + ( p b + q b + r b ) × ( − x b − y b − z b ) Q 1 2 3 1 2 3 a 1 a 2 a 3 = ( u + ry − qz ) b + ( v − rx + pz ) b + ( w + qx − py ) b a a 1 a a 2 a a 3 = x d + y d + z d (22) t 1 t 2 t 3 where       b − sin Λ cos Γ − cos Λ sin Λ sin Γ d 1 1             = (23) b − cos Λ cos Γ sin Λ cos Λ sin Γ d  2     2  b − sin Γ 0 − cos Γ d 3 3     x − ( u + ry − qz ) sin Λ cos Γ − ( v − rx + pz ) cos Λ cos Γ − ( w + qx − py ) sin Γ t a a a a a a         = (24) y − ( u + ry − qz ) cos Λ + ( v − rx + pz ) sin Λ t a a a a     z ( u + ry − qz ) sin Λ sin Γ + ( v − rx + pz ) cos Λ sin Γ − ( w + qx − py ) cos Γ t a a a a a a and ( p , q , r ) are aircraft angular velocity components in the roll, pitch, and yaw axes, and ( x , y , z ) is the coordinate a a a of point Q in the aircraft body-fixed reference frame B relative to the aircraft C.G. (center of gravity) such that x is a positive when point Q is aft of the aircraft CG, y is positive when point Q is toward the left wing from the aircraft a C.G., and z is positive when point Q is above the aircraft C.G.

a Consider a trim problem when β = 0, p = q = r = 0. Then     x − u sin Λ cos Γ − w sin Γ t         = (25) y − u cos Λ t     z u sin Λ sin Γ − w cos Γ t The local velocity at point Q due to aircraft rigid-body dynamics and aeroelastic deflections in the reference frame D is obtained as ∂ r ∂ φ v = v + + × ( r + φ × q ) = v d + v d + v d (26) Q x 1 y 2 z 3 ∂ t ∂ t 7 of 26 American Institute of Aeronautics and Astronautics where     v x + U − ( z + W + y Θ ) W − ( y + V − z Θ ) V x t t xt xt         = (27) v y + V + ( U − yV − zW ) V − ( z + W + y Θ ) Θ y t t x x xt t     v z + W + ( U − yV − zW ) W + ( y + V − z Θ ) Θ z t t x x xt t In order to compute the aeroelastic forces and moments, the velocity must be transformed from the reference frame D to the airfoil local coordinate reference frame defined by ( μ , η , ξ ) as shown in Figure 2. Then the transformation can be performed using successive rotation matrix multiplication operations as           v 1 0 0 cos V sin V 0 cos W 0 sin W v μ x x x x x                     = v 0 cos ( Θ + γ ) sin ( Θ + γ ) − sin V cos V 0 0 1 0 v η x x y           v 0 − sin ( Θ + γ ) cos ( Θ + γ ) 0 0 1 − sin W 0 cos W v x x z ξ   cos V ( v cos W + v sin W ) + v sin V x x x z x y x     = cos ( Θ + γ ) [ − sin V ( v cos W + v sin W ) + v cos V ] + sin ( Θ + γ ) ( − v sin W + v cos W ) x x x z x y x x x z x   − sin ( Θ + γ ) [ − sin V ( v cos W + v sin W ) + v cos V ] + cos ( Θ + γ ) ( − v sin W + v cos W ) x x x z x y x x x z x   v + v V + v W x y x z x     ≈ (28) − v [ V + W ( Θ + γ )] + v + v [( Θ + γ ) − V W ] x x x y z x x   v [ − W + V ( Θ + γ )] − v ( Θ + γ ) + v [ 1 + V W ( Θ + γ )] x x x y z x x for small deflections.

The local aeroelastic angle of attack on the airfoil section due to the velocity components v and v in the reference η ξ frame D, as shown in Fig. 8, is computed as ( ) v + w ¯ v + ∆ v + w v + w ¯ v + w ∆ v ξ i ξ ξ i ξ i ξ i η α = = = − (29) c v ¯ v + ∆ v ¯ v ¯ v η η η η η where w is the downwash due to the three-dimensional lift distribution over a finite-aspect ratio wing and i ¯ v = u sin Λ sin Γ − w cos Γ (30) ξ ∆ v = W + ( U − yV − zW ) W + ( y + V − z Θ ) Θ + v [ − W + V ( Θ + γ )] − v ( Θ + γ ) + v V W ( Θ + γ ) (31) t x x xt t x x x y z x x ξ ¯ v = − u cos Λ (32) η ∆ v = V + ( U − yV − zW ) V − ( z + W + y Θ ) Θ − v [ V + W ( Θ + γ )] + v [( Θ + γ ) − V W ] (33) η t x x xt t x x x z x x Assuming that z ≈ 0 and neglecting any third-order terms or higher, the local aeroelastic angle of attack is evaluated 8 of 26 American Institute of Aeronautics and Astronautics as u sin Λ sin Γ − w cos Γ + w W + ( U − yV ) W + ( y + V ) Θ i t x xt t α = − − c u cos Λ u cos Λ [ − u sin Λ cos Γ − w sin Γ + U − ( W + y Θ ) W − ( y + V ) V ] [ − W + V ( Θ + γ )] t xt xt x x − u cos Λ [ − u cos Λ + V + ( U − yV ) V − ( W + y Θ ) Θ ] ( Θ + γ ) t x xt t + u cos Λ { u sin Λ sin Γ − w cos Γ + w i − V + ( U − yV ) V − ( W + y Θ ) Θ t x xt t 2 2 u cos Λ − [ − u sin Λ cos Γ − w sin Γ + U − ( W + y Θ ) W − ( y + V ) V ] [ V + W ( Θ + γ )] t xt xt x x } + [ u sin Λ sin Γ − w cos Γ + W + ( U − yV ) W + ( y + V ) Θ ] [( Θ + γ ) − V W ] (34) t x xt t x x Let u ≈ V , w ≈ V α , and w = V α , where α is the induced angle of attack along the wing span. Then the ∞ ∞ i ∞ i i nonlinear partial derivatives of the local aeroelastic angle of attack are evaluated as ( sin Λ sin Γ + α ) ( 1 + γ tan Λ sin Γ ) i α = − γ − (35) cos Λ [ ] ∂ α cos Γ γ ( 2 sin Λ sin Γ + α ) c i = 1 + (36) ∂ α cos Λ cos Λ ∂ α γ cos Γ c = − (37) 2 2 ∂ α cos Λ ∂ α ( sin Λ sin Γ − α cos Γ ) ( sin Λ sin Γ − α cos Γ + α ) c i = − 1 − (38) ∂ Θ cos Λ [ ] ∂ α sin Λ cos Γ + α sin Γ γ ( sin Λ sin Γ − α cos Γ + α ) c i = − 1 + (39) ∂ W cos Λ cos Λ x ( ) ∂ α sin Λ cos Γ + α sin Γ sin Λ sin Γ − α cos Γ + α c i = γ − (40) ∂ V cos Λ cos Λ x ∂ α W V ( sin Λ sin Γ − α cos Γ + α ) c x x i = + (41) ∂ U V cos Λ V cos Λ t ∞ ∞ ∂ α y + V ( y γ − W ) ( sin Λ sin Γ − α cos Γ + α ) c i = − − (42) ∂ Θ V cos Λ V cos Λ t ∞ ∞ ∂ α 1 ( Θ + γ ) ( sin Λ sin Γ − α cos Γ + α ) c i = − − (43) ∂ W V cos Λ V cos Λ t ∞ ∞ ∂ α Θ + γ sin Λ sin Γ − α cos Γ + α c i = − (44) ∂ V V cos Λ V cos Λ t ∞ ∞ ∂ α ( sin Λ cos Γ + α sin Γ ) ( sin Λ sin Γ − α cos Γ + α ) c i = − (45) ∂ Θ W cos Λ x ∂ α sin Λ cos Γ + α sin Γ c = (46) ∂ Θ V cos Λ x ∂ α ( sin Λ sin Γ − α cos Γ ) ( sin Λ sin Γ − α cos Γ + α ) c i = (47) ∂ W V cos Λ x x ∂ α U − yV c x = − (48) ∂ W V cos Λ xt ∞ ∂ α yW U ( sin Λ sin Γ − α cos Γ + α ) c x i = − − (49) ∂ V V cos Λ V cos Λ xt ∞ ∞ We note that the time derivative of the axial displacement contributes to the aeroelastic angle of attack as a damping 9 of 26 American Institute of Aeronautics and Astronautics source. Otherwise, the axial displacement U is weakly coupled to the wing aeroelasticity.

The nonlinear angle of attack is then expressed as α ( x , y ) = α ( x ) + α ( x , y ) (50) c r e where α ( x ) is the rigid angle of attack and α ( x , y ) is the elastic angle of attack r e ∂ α ∂ α c c α ( x , y ) = α + α + α (51) r 0 ∂ α ∂ α ∂ α ∂ α ∂ α ∂ α ∂ α ∂ α c c c c c c α ( x , y ) = Θ + W + V + Θ W + Θ V + W V e x x x x x x ∂ Θ ∂ W ∂ V ∂ Θ W ∂ Θ V ∂ W V x x x x x x ∂ α ∂ α ∂ α ∂ α ∂ α ∂ α c c c c c c + U + Θ + W + V + W + V (52) t t t t xt xt ∂ U ∂ Θ ∂ W ∂ V ∂ W ∂ V t t t t xt xt The elastic angle of attack can be linearized as a function of the static deflection as ( ) ( ) ∂ α ∂ α ∂ α ∂ α ∂ α ∂ α c c c c c c ¯ ¯ ¯ ¯ α ( x , y ) = + W + V Θ + + Θ + V W e x x x x ∂ Θ ∂ Θ W ∂ Θ V ∂ W ∂ Θ W ∂ W V x x x x x x ( ) ∂ α ∂ α ∂ α ∂ α ∂ α ∂ α ∂ α ∂ α ∂ α c c c c c c c c c ¯ ¯ + + Θ + W V + U + Θ + V + W + V + W (53) x x t t t t xt xt ∂ V ∂ Θ V ∂ W V ∂ U ∂ Θ ∂ V ∂ W ∂ V ∂ W x x x x t t t t xt xt Finally, if we ignore all the deflection-dependent terms, then the linearized elastic angle of attack becomes ∂ α ∂ α ∂ α ∂ α ∂ α ∂ α c c c c c c α ( x , y ) = Θ + W + V + Θ + W + V (54) e x x t t t ∂ Θ ∂ W ∂ V ∂ Θ ∂ W ∂ V x x t t t IV. Aerodynamic Forces and Moments The wing section unsteady lift coefficient is given by π ˙ α cos Λ c ac c = c [ α + C ( k ) α ] cos Λ + (55) L L r ac α 2 V ∞ ω c where k = is the reduced frequency parameter, ω is the frequency of wing oscillations, c is the section chord in the 2 V ∞ streamwise direction, c is the section lift curve slope, α is the effective elastic angle of attack at the three-quarter- L ac α 5 8 chord point, and C ( k ) is the Theodorsen’s complex-valued function.

The section pitching moment coefficient about the aircraft pitch axis is evaluated as e e π ˙ α c u ac c = c + c [ α + C ( k ) α ] − (56) m m L r ac ac α c c 2 V ∞ where c is the section pitching moment coefficient at the aerodynamic center, positive nose up, e is the offset of the m ac aerodynamic center forward of the elastic axis, and e is the offset of the three-quarter chord point aft of the elastic u axis.

The section drag coefficient is expressed in a parabolic drag polar form as c = c + Kc (57) D D L 10 of 26 American Institute of Aeronautics and Astronautics where c is the section parasitic drag coefficient and K is a parabolic drag parameter.

D The force and moment coefficients can be written as ∂ c ∂ c ∂ c ∂ c ∂ c ∂ c ∂ c ∂ c L L L L L L L L c = c + Θ + W + V + Θ + W + V + W + V L L x x t t t xt xt r ∂ Θ ∂ W ∂ V ∂ Θ ∂ W ∂ V ∂ W ∂ V x x t t t xt xt ∂ c ∂ c ∂ c L L L + Θ + W + V (58) tt tt tt ∂ Θ ∂ W ∂ V tt tt tt ∂ c ∂ c ∂ c ∂ c ∂ c ∂ c ∂ c ∂ c m m m m m m m m c = c + Θ + W + V + Θ + W + V + W + V m m x x t t t xt xt r ∂ Θ ∂ W ∂ V ∂ Θ ∂ W ∂ V ∂ W ∂ V x x t t t xt xt ∂ c ∂ c ∂ c m m m + Θ + W + V (59) tt tt tt ∂ Θ ∂ W ∂ V tt tt tt ∂ c D c = c + ( c − c ) = c + 2 Kc ( c − c ) (60) D D L L D L L L r r r r r ∂ c L where c , c , and c are the rigid lift, pitching moment, and drag coefficient contributions.

L m D r r r The lift force, drag force, and pitching moment in the aircraft coordinate reference frame along the elastic axis are expressed as l = c q cos Λ c (61) L ∞ d = c q cos Λ c (62) D ∞ m = c q cos Λ c (63) m ∞ The axial and normal force coefficients are defined as c = c cos α + c sin α ≈ c (64) n L D L c = − c sin α + c cos α (65) a L D The forces and moments in the local coordinate reference frame are obtained as a f = c q c cos Λ sin Γ + c q c cos Λ sin Λ cos Γ ≈ 0 (66) n ∞ a ∞ x a 2 f = c q c cos Λ (67) a ∞ y a f ≈ c q c cos Λ cos Γ − c q c cos Λ sin Λ sin Γ ≈ c q c cos Λ cos Γ (68) n ∞ a ∞ L ∞ z a 2 2 m = − c q c cos Λ cos Γ (69) m ∞ x a 2 m = c q c sin Λ cos Λ (70) m ∞ y a 2 2 m = c q c cos Λ sin Γ (71) m ∞ z A. Propulsive Forces and Moments In addition, the propulsive effects of the aircraft engines must be accounted for in the analysis. Both the engine mass and thrust can contribute to the aeroelasticity. The propulsive force and moment vector are expressed in the reference frame D as ( ) e e e e f = δ ( x − x ) F d + F d + F d e 1 2 3 x y z = δ ( x − x ) [( − T sin Λ cos Γ − m g sin Γ ) d − T cos Λ d + ( T sin Λ sin Γ − m g cos Γ ) d ] (72) e e 1 2 e 3 e m = r × f = ( x d − y d − z d ) × f (73) e e e 1 e 2 e 3 e 11 of 26 American Institute of Aeronautics and Astronautics where T is the engine thrust, m is the engine mass, ( x , y , z ) is the coordinate of the engine thrust center such that e e e e y > 0 forward of the elastic axis and z > 0 below the elastic axis, and δ ( x − x ) is the Dirac delta function such that e e e ∫ δ ( x − x ) f ( x ) dx = f ( x ) (74) e e Transforming into the local coordinate reference frame and neglecting nonlinear contributions, the propulsive forces and moments are given by ( ) e e e e f = δ ( x − x ) F + F W + F V (75) e x x x x z y [ ] e e e e f = δ ( x − x ) F + F ( Θ + γ ) − F V (76) e x y y z x [ ] e e e e f = δ ( x − x ) F − F ( Θ + γ ) − F W (77) e x z z y x [ ] e e e m = δ ( x − x ) F z − F y + ( F y + F x ) W − ( F z + F x ) V (78) e e e x e y e x x e z e x x y z [ ] e e e m = δ ( x − x ) − F z − F x + ( F y + F x ) ( Θ + γ ) − ( F z − F y ) V (79) e e e x e y e y e z e x y x z e m = δ ( x − x ) [ F y + F x + ( F z + F x ) ( Θ + γ ) − ( F z − F y ) W ] (80) e x e y e x e z e y e z e x z B. Inertial Forces and Moments The inertial forces and moments are due to acceleration experienced by the wing, the air volume surrounding it which gives rise to the apparent mass acting at mid-chord, as well as the engine mass and half of the mass of the fuselage and tail empennage. The acceleration results in the inertial forces and moments as follows: π 1 i 2 f = − mg sin Γ − mU − ρ c cos Λ U − δ ( x ) m g sin Γ (81) tt ∞ tt f x 4 2 ( ) π 1 1 i 2 f = − mV − ρ c cos Λ V + δ ( x − x ) ( − m z Θ − m V ) + δ ( x ) m z Θ − m V (82) tt ∞ tt e e e tt e tt f f tt f tt y 4 2 2 π π i 2 2 f = − mg cos Γ + me Θ + ρ c cos Λ e Θ − mW − ρ c cos Λ W + δ ( x − x ) ( m y Θ − m W ) cg tt ∞ m tt tt ∞ tt e e e tt e tt z 4 4 ( ) 1 1 1 + δ ( x ) − m g cos Γ − m y Θ − m W (83) f f f tt f tt 2 2 2 π π π i 2 2 2 4 2 m = mge cos Γ − mr Θ − ρ c e cos ΛΘ − ρ c cos ΛΘ + me W + ρ c e cos Λ W cg tt ∞ tt ∞ tt cg tt ∞ m tt x k m 4 128 4 ( ) ( ) 1 1 1 1 2 f + δ ( x − x ) − m r Θ + m y W − m z V + δ ( x ) − m gy cos Γ − I Θ − m y W + m z V (84) e e tt e e tt e e tt f f tt f f tt f f tt e xx 2 2 2 2 ( ) 1 1 i f m = δ ( x ) − m gz sin Γ + I W (85) f f xtt y yy 2 2 ( ) 1 1 i f m = − mge sin Γ + δ ( x ) m gy sin Γ − I V (86) cg f f xtt z zz 2 2 where e is the offset of the mid-chord aft of the elastic axis; e is the offset of the center of mass forward of the m cg ( ) elastic axis; 0 , y , z is the coordinate of the aircraft CG such that y > 0 aft of elastic axis and z > 0 above the f f f f f f f elastic axis; m is the mass of the fuselage and tail empennage; and I , I , and I are the inertias about the left wing f xx yy zz reference frame of the fuselage and tail empennage.

The aircraft CG is located in the aircraft body-fixed reference frame B at a distance x aft of the elastic axis and z B B 12 of 26 American Institute of Aeronautics and Astronautics above the elastic axis at x = 0. Then by coordinate transformation into the reference frame D y = x cos Λ (87) f B z = − x sin Λ sin Γ + z cos Γ (88) f B B The roll, pitch, yaw, and roll-yaw inertias of the fuselage and the tail empennage at the elastic axis at x = 0 are computed as f f ¯ I = I (89) XX XX ( ) f f 2 2 ¯ I = I + m x + z (90) f B B YY YY f f ¯ I = I + m x (91) f B ZZ ZZ f f ¯ I = I + m x z (92) f B B XZ XZ f f f f where ¯ I , ¯ I , ¯ I , and ¯ I are the roll, pitch, yaw, and roll-yaw inertias at the aircraft CG.

XX YY ZZ XZ The inertias can be transformed into the wing reference frame D as     f f I 0 − I − sin Λ cos Γ XX XZ [ ]     f f     I = − sin Λ cos Γ − cos Λ cos Γ − sin Γ 0 I 0 − cos Λ cos Γ xx  YY    f − I 0 I − sin Γ ZZ XZ f f f 2 2 2 2 2 = I sin Λ cos Γ + I cos Λ cos Γ + I sin Γ − I sin Λ sin 2 Γ (93) ZZ XX YY XZ     f f I 0 − I − cos Λ XX XZ [ ]     f f f 2 2 f     I = = I cos Λ + I sin Λ (94) − cos Λ sin Λ 0 0 I 0 sin Λ yy XX YY     YY f − I 0 I 0 ZZ XZ     f f I 0 − I sin Λ sin Γ XX XZ [ ]     f f     I = sin Λ sin Γ cos Λ sin Γ − cos Γ 0 I 0 cos Λ sin Γ zz     YY f − I 0 I − cos Γ ZZ XZ f 2 2 f 2 2 2 f = I sin Λ sin Γ + I cos Λ sin Γ + I cos Γ + I sin Λ sin 2 Γ (95) ZZ XX YY XZ V. Aeroelastic Equations for Combined Extensional, Bending, and Torsion Motion A. Equilibrium Conditions Consider the equilibrium of forces and moments acting on a wing section. The resulting equilibrium equations are given by ∂ P x + f = 0 (96) x ∂ x ∂ P y + f = 0 (97) y ∂ x ∂ P y + f = 0 (98) y ∂ x ∂ M x − P W + P V + m = 0 (99) y x z x x ∂ x 13 of 26 American Institute of Aeronautics and Astronautics ∂ M y − P W + P − m = 0 (100) x x z y ∂ x ∂ M z − P V + P + m = 0 (101) x x y z ∂ x These equations become ( ) ( ) ∂ M ∂ M ∂ M x z y + − P V + m W − − P W − m V + m = 0 (102) x x z x x x y x x ∂ x ∂ x ∂ x ∂ M ∂ ( P W ) ∂ m y y x x − − − f = 0 (103) z ∂ x ∂ x ∂ x ∂ M ∂ ( P V ) ∂ m z x x z − + − f = 0 (104) y ∂ x ∂ x ∂ x ¯ ¯ ¯ ¯ Let U ( x , t ) = U ( x ) + u ( x , t ) , Θ ( x , t ) = Θ ( x ) + θ ( x , t ) , W ( x , t ) = W ( x ) + w ( x , t ) , and V ( x , t ) = V ( x ) + v ( x , t ) where the bar symbols denote the static deflection. The static tensile force acting on the wing section is defined as ¯ T = EA U . Then, considering the linear deflection terms, we have x     P EAu − EAe v x x y xx     ( )     M GJ + T k θ  x   x  = (105)         M EI w − EI v y yy xx yz xx     M − EAe u − EI w + EI v z y x yz xx zz xx ∂ ( EAu − EAe v ) x y xx = − f (106) x ∂ x [( ) ] ( ) ( ) ¯ ∂ GJ + T k θ x d M ∂ M z z ¯ ¯ ¯ ¯ ¯ = − m − − P V + ¯ m w − − P v − P V + m W x x x z x x x x x z x ∂ x dx ∂ x ( ) ( ) ¯ d M ∂ M y y ¯ ¯ ¯ ¯ ¯ + − P W − ¯ m v + − P w − P W − m P V (107) x x y x x x x x y x dx ∂ x ¯ ¯ ∂ ( EI w − EI v ) ∂ m ∂ ( P w + P W ) yy xx yz xx y x x x x = f + + (108) z ∂ x ∂ x ∂ x ¯ ¯ ∂ ( − EAe u − EI w + EI v ) ∂ m ∂ ( P v + P V ) y x yz xx zz xx z x x x x = f − + (109) y ∂ x ∂ x ∂ x where I = Ak .

xx Since ¯ V and ¯ V are generally small, so they can be neglected. Then, ¯ P = T so that the equations above become x xx x [( ) ] [ ] ¯ ∂ GJ + T k θ d ( − Te − EI W ) x y yz xx = − m − + ¯ m w x z x ∂ x dx [ ] [ ] ¯ ∂ ( − EAe u − EI w + EI v ) d ( EI W ) y x yz xx zz xx yy xx ¯ ¯ − − T v + m W + − T W − ¯ m v (110) x z x x y x ∂ x dx ¯ ∂ ( EI w − EI v ) ∂ m ∂ [ Tw + ( EAu − EAe v ) W ] yy xx yz xx y x x y xx x = f + + (111) z ∂ x ∂ x ∂ x 14 of 26 American Institute of Aeronautics and Astronautics ∂ ( − EAe u − EI w + EI v ) ∂ m ∂ ( T v ) y x yz xx zz xx z x = f − + (112) y ∂ x ∂ x ∂ x Finally, if the vertical bending deflection is not too large so that the nonlinear effect is not dominant, then in the simplest form these equations are expressed as [( ) ] ∂ GJ + T k θ x = − m (113) x ∂ x ∂ ( EI w − EI v ) ∂ m ∂ ( Tw ) yy xx yz xx y x = f + + (114) z ∂ x ∂ x ∂ x ∂ ( − EAe u − EI w + EI v ) ∂ m ∂ ( T v ) y x yz xx zz xx z x = f − + (115) y ∂ x ∂ x ∂ x The aeroelastic equations of motion for axially loaded structures with no deflection-dependent terms are given by ( ) ( ) ∂ ( EAu − EAe v ) π x y xx 2 e e = m + ρ c cos Λ u − δ ( x − x ) F w + F v (116) ∞ tt e x x z y ∂ x 4 [( ) ] 2 ( ) ∂ GJ + T k θ π π x 2 2 2 2 2 4 = c q c cos Λ cos Γ + mr + ρ c e cos Λ + ρ c cos Λ θ m ∞ ∞ ∞ tt k m ∂ x 4 128 ( ) π − me + ρ c e cos Λ w − δ ( x − x ) [( F y + F x ) w − ( F z + F x ) v ] cg ∞ m tt e x e y e x x e z e x ( ) ( ) 1 1 1 2 f + δ ( x − x ) m r θ − m y w + m z v + δ ( x ) m y w + I θ − m z v (117) e e tt e e tt e e tt f f tt tt f f tt e xx 2 2 2 ( ) ( ) ∂ ( EI w − EI v ) ∂ π yy xx yz xx 2 2 = c q c cos Λ + c q c sin Λ cos Λ + me + ρ c e cos Λ θ L ∞ m ∞ cg ∞ m tt ∂ x ∂ x 4 ( ) ( ) π 2 e e − m + ρ c cos Λ w + δ ( x − x ) − F θ − F w ∞ tt e x y x { } ∂ + δ ( x − x ) [( F y + F x ) θ − ( F z − F y ) v ] + δ ( x − x ) ( m y θ − m w ) e x e y e y e z e x e e e tt e tt ∂ x ( ) [ ( )] 1 1 ∂ 1 ∂ ( Tw ) x f + δ ( x ) − m y θ − m w + δ ( x ) I w + (118) f f tt f tt xtt yy 2 2 ∂ x 2 ∂ x ( ) ( ) ∂ ( − EAe u − EI w + EI v ) ∂ π y x yz xx zz xx 2 2 2 2 = c q c cos Λ − c q c cos Λ sin Γ − m + ρ c cos Λ v a ∞ m ∞ ∞ tt ∂ x ∂ x 4 ( ) { } ∂ e e + δ ( x − x ) F θ − F v − δ ( x − x ) [( F z + F x ) θ − ( F z − F y ) w ] e x e x e z e y e z e x z x ∂ x ( ) [ ( )] 1 1 ∂ 1 ∂ ( T v ) x f + δ ( x − x ) ( − m z θ − m v ) + δ ( x ) m z θ − m v − δ ( x ) − I v + (119) e e e tt e tt f f tt f tt xtt zz 2 2 ∂ x 2 ∂ x 15 of 26 American Institute of Aeronautics and Astronautics B. Discretization The partial differential equations can be solved by weak-form solutions such as the finite-element method (FEM).

Using interpolation functions, the displacements can be expressed as n u ( x , t ) = N ( x ) u ( t ) (120)

∑ u i

i = 1 n θ ( x , t ) = N ( x ) θ ( t ) (121)

∑ θ i

i = 1 n w ( x , t ) = N ( x ) φ ( t ) (122) w i

∑

i = 1 n v ( x , t ) = N ( x ) ψ ( t ) (123) v i

∑

i = 1 [ ] [ ] > > where N ( x ) , N ( x ) , N ( x ) , and N ( x ) are the interpolation functions, φ = , ψ = , and u θ w v i w w i v v i x , i i x , i n is the number of elements.

The discretization results in a spring-mass-damper system in the form n n n M ¨ x + C ˙ x + K x = 0 (124) i i i i i i

∑ ∑ ∑

i = 1 i = 1 i = 1 where M is the elemental mass matrix, C is the elemental damping matrix, K is the elemental stiffness matrix, and i i i [ ] > x = is the nodal displacement vector.

i u θ w w v v i i i x , i i x , i The mass matrix due to the structural mass, aerodynamic mass, engine mass, and fuselage mass is given by f s a e M = M + M + M + M (125) i i i i i 16 of 26 American Institute of Aeronautics and Astronautics where   > N N 0 0 0 u u   ∫   > 2 > l i 0 N r N − N e N 0  cg w  θ s θ k θ M = m   dx (126) i  > >  0 − N e N N N 0 cg θ w w w   > 0 0 0 N N v v   > N N 0 0 0 u u ( )   ∫  > 2 c >  l i 0 N e + N − N e N 0 π θ m w  m  a 2 θ 32 θ M = ρ c cos Λ   dx ∞ i 4   0 > > 0 − N e N N N 0  m θ w  w w > 0 0 0 N N v v     0 0 0 0       ∫   ∂ c ∂ c ∂ c  > m > m > m l i 0 N c cos Γ N N c cos Γ N N c cos Γ N  θ w v  θ θ θ 2 ∂ θ ∂ w ∂ v tt tt tt + q c cos Λ   ∞ ∂ c ∂ c ∂ c   > 1 L > 1 L > 1 L   0 − N N − N N − N N θ w v  w w w  cos Λ ∂ θ cos Λ ∂ w cos Λ ∂ v  tt tt tt    > ∂ c > ∂ c > ∂ c a a a 0 − N N − N N − N N θ v w v v v ∂ θ ∂ v ∂ w tt tt tt     0 0 0 0          0 0 0 0   + dx (127)   ′ ′ ′ > ∂ c > ∂ c > ∂ c  m m m   0 N c tan Λ N N c tan Λ N N c tan Λ N  w v w θ w w  ∂ θ ∂ w ∂ v   tt tt tt   ′ ′ ′  > ∂ c > ∂ c > ∂ c m m m 0 − N c sin Γ N − N c sin Γ N − N c sin Γ N θ w v v v v ∂ θ ∂ w ∂ v tt tt tt   0 0 0 0     > 2 > > 0 N r N − N y N N z N  θ e w e v  e θ e θ θ M = m (128)   e i  > >  0 − N y N N N 0 e θ w  w w  > > 0 N z N 0 N N e θ v v v x = x e   0 0 0 0    f  > > > ¨ 1 0 N I N Θ N m y N − N m z N  xx θ i f f w f f v  f θ θ θ M = (129)   i ′ ′ f 2  > > >  0 N m y N N m N + N I N 0 f f θ f w yy  w w w w  ′ ′ f > > > 0 − N m z N 0 N m N + N I N f f θ f v zz v v v v x = 0 The damping matrix due to structural damping and aerodynamic damping is given by s a C = C + C (130) i i i 17 of 26 American Institute of Aeronautics and Astronautics where s s − 1 C = 2 M Φ ζ ΩΦ (131) i i     0 0 0 0   ( ) ( )     ′ ′   ∂ c ∂ c ∂ c ∂ c ∂ c  > m > m m > m m ∫  0 N c cos Γ N N c cos Γ N + N N c cos Γ N + N l i  θ w v  θ θ w θ v ∂ θ ∂ w ∂ w ∂ v ∂ v t t xt t xt a 2  ( ) ( )  C = q c cos Λ dx ∞ i ′ ′   1 ∂ c 1 ∂ c ∂ c 1 ∂ c ∂ c > L > L L > L L  0 − N N − N N + N − N N + N    θ w v w w w w v  cos Λ ∂ θ cos Λ ∂ w ∂ w cos Λ ∂ v ∂ v t t xt t xt    ( ) ( )   ′ ′ > ∂ c > ∂ c ∂ c > ∂ c ∂ c  a a a a a 0 − N N − N N + N − N N + N θ w v v v w v v ∂ θ ∂ w ∂ w ∂ v ∂ v t t xt t xt      0 0 0 0           0 0 0 0   ( ) ( )  + dx (132) ′ ′ ′ ′ ′   > ∂ c > ∂ c ∂ c > ∂ c ∂ c m m m m m  0 N c tan Λ N N c tan Λ N + N N c tan Λ N + N  w v   w θ w w w v ∂ θ ∂ w ∂ w ∂ v ∂ v  t t xt t xt   ( ) ( )   ′ ′ ′ ′ ′  > ∂ c > ∂ c ∂ c > ∂ c ∂ c m m m m m  0 − N c sin Γ N − N c sin Γ N + N − N c sin Γ N + N θ w v v v w v v ∂ θ ∂ w ∂ w ∂ v ∂ v t t xt t xt where ζ = diag ( ζ , ζ , . . . , ζ ) is the diagonal viscous damping ratio matrix, Φ is the structural dynamic eigenvector 1 2 n matrix, and Ω = diag ( ω , ω , . . . , ω ) is a diagonal matrix of the structural frequencies.

1 2 n The stiffness matrix due to structural stiffness including the nonlinear tension-induced stiffness, aerodynamic stiffness, and engine thrust-induced stiffness is given by s a e K = K + K + K (133) i i i i where   ′ ′ ′ ′′ > > N EAN 0 0 − N EAe N y u u u v   ( ) ∫ ′ ′   > 2 l i 0 N GJ + T k N 0 0   s θ θ K =   dx (134) i ′′ ′′ ′ ′ ′′ ′′  > > >  0 0 N EI N + N T N − N EI N yy yz w w w w w v   ′′ ′ ′′ ′′ ′′ ′′ ′ ′ > > > > − N EAe N 0 − N EI N N EI N + N T N y yz zz v u v w v v v v     0 0 0 0       ′ ′ ∫ ∂ c ∂ c ∂ c   > m > m > m  l i 0 N c cos Γ N N c cos Γ N N c cos Γ N  θ  w v a 2 θ ∂ θ θ ∂ w θ ∂ v x x K = q c cos Λ   dx ∞ i ′ ′ ∂ c ∂ c ∂ c   > 1 > 1 > 1  L L L  0 − N N − N N − N N θ  w w w w v  cos Λ ∂ θ cos Λ ∂ w cos Λ ∂ v   x x   ′ ′ > ∂ c > ∂ c > ∂ c a a a 0 − N N − N N − N N v θ v w v v ∂ θ ∂ w ∂ v x x     0 0 0 0          0 0 0 0   + dx (135)   ′ ′ ′ ′ ′ > ∂ c > ∂ c > ∂ c  m m m   0 N c tan Λ N N c tan Λ N N c tan Λ N  θ  w w w w v   ∂ θ ∂ w ∂ v x x   ′ ′ ′ ′ ′  ∂ c ∂ c ∂ c > m > m > m 0 − N c sin Γ N − N c sin Γ N − N c sin Γ N θ v v w v v ∂ θ ∂ w ∂ v x x   ′ ′ > e > e 0 0 − N F N − N F N u z w u y v   ′ ′   > > 0 0 − N ( F y + F x ) N N ( F z + F x ) N  x e y e x e z e  e w v θ θ K = (136)   i ′ ′ ′ ′  > e > > e >  0 N F N + N ( F y + F x ) N N F N − N ( F z − F y ) N θ x e y e θ y e z e  w y w w x w w v  ′ ′ ′ ′ > e > > > e 0 − N F N − N ( F z + F x ) N N ( F z − F y ) N N F N θ x e z e θ y e z e v z v v w v x v The effect of the tensile force can be clearly seen in the structural stiffness as the additive stiffness terms in both the 18 of 26 American Institute of Aeronautics and Astronautics bending and torsion stiffness matrices. The presence of a tensile force causes a geometric nonlinear tension stiffening effect which increases the natural frequencies of the structural dynamic modes. Conversely, a compressive force causes the natural frequencies to decrease. Static instability due to compressive loading could result when the compressive force-induced stiffness is large enough to offset the torsion or bending stiffness.

VI. Applications A. Rotary Wings Rotary wings are axially loaded structures that experience tensile loading due to the centrifugal acceleration. The structural stiffness of a rotary wing is given by [ ]   ( ) ′ ′ ′ ′ ′ ′′ ′ ′ ′′ > > > N GJ + EB γ N − N EB γ N − N EB γ N 1 2 3 w v θ θ θ θ ∫   l i   s ′′ ′ ′ ′′ ′′ ′′ ′′ K =   dx > > > − N EB γ N N EI N − N EI N 2 yy yz  w w w w v  0 θ ′′ ′ ′ ′′ ′′ ′′ ′′ > > > − N EB γ N − N EI N N EI N 3 yz zz v v w v v θ   ′ ′ > 2 N T k N 0 0 ∫ θ θ l i   ′ ′  >  + Ω dx (137) 0 N T N 0  w w  ′ ′ > 0 0 N T N v v ∫ L where T ( x ) = mr Ω dr is the centrifugal force, r is the radius, and Ω is the angular speed.

x Therefore, the system matrix equations at the angular speed Ω and zero speed can be expressed as ( ) 2 − 1 2 − 1 − ω x + M K + Ω M F x = 0 (138) n and 2 − 1 − ω x + M K x = 0 (139) Thus, the natural frequencies of a rotary wing are dependent on the angular speed of rotation according to 2 2 2 ω = ω + ( f Ω ) (140) n 0 where ω is the natural frequency at the angular speed Ω and ω is the zero-speed natural frequency, and f is the n 0 centrifugal stiffening factor.

As an example, the natural frequencies of compressor rotor blades in the NASA Ames 11-Ft Unitary Plan 11-Foot By 11-Foot Transonic Wind Tunnel three-stage compressor are computed. The experimental natural frequencies were obtained from modal analysis measurements. In addition, wind-on strain gauge measurements were conducted at various compressor speeds to obtain the variation in the natural frequencies with the compressor speed. The natural frequencies were computed from the Fast Fourier Transform (FFT) of the time signatures of the strain gauges at various compressor speeds. This allows the centrifugal stiffening factor to be computed. Table 1 shows the computed natural frequencies of the three-stage compressor rotor blades by the finite-element method presented herein as compared to the experimental results. The agreement between the computed natural frequencies and the experimental results is excellent.

19 of 26 American Institute of Aeronautics and Astronautics Mode FEM ω (Hz) Bench Test ω (Hz) FEM ω at 700 rpm (Hz) FEM f Wind-On f 0 0 n 1B 54.43 54.18 61.67 2.488 2.375 2B 167.32 160.92 174.51 4.252 4.842 1T 302.37 305.56 309.25 5.560 4.661 3B 370.29 374.61 375.25 5.24 9.110 2T 549.98 554.51 561.91 9.87 8.764 Table 1 - Natural Frequencies of Rotor Blades in NASA Ames Unitary Plan 11-Foot By 11-Foot Transonic Wind Tunnel Three-Stage Compressor B. Truss-Braced Wing Aircraft A recent flutter wind tunnel test of a dynamically scaled TBW was conducted in NASA Langley Transonic Dynamic Wind Tunnel (TDT). The model is 16% scaled, semi-span, wall-mounted model as shown in Fig. 3. It was reported that the flutter speed is dependent on the angle of attack at ± 1 deg. While there may be many nonlinear factors that could contribute to this observation, one possible exploration could be the geometric nonlinear tension stiffening effect of the main struts at a positive angle of attack, and conversely the softening effect at a negative angle of attack, assuming the axial force is below the critical buckling load.

Fig. 3 - Truss-Braced Wing (TBW) Model in NASA Langley Transonic Dynamic Tunnel (TDT) A static aeroelastic analysis is performed to compute the static deflections and the tension carried by the struts.

The finite-element model (FEM) of the TBW utilizes an equivalent beam/frame configuration that uses a total of five components. Only the right side of the aircraft is modeled, and the wing configuration is assumed to be cantilevered to the side of the fuselage. The three main physical components of the wing configuration are the wing, the strut, and the jury strut. The FEM uses 30 elements for the wing, 28 for the strut, and 5 for the jury strut. Two rigid elements are used to connect the strut and the jury strut to the wing. Figure 4 shows the FEM of the TBW.

20 of 26 American Institute of Aeronautics and Astronautics Fig. 4 - Finite-Element Model of Truss-Braced Wing (TBW) The FEM uses lumped mass and inertia properties provided by Boeing Research & Technology. The lumped mass data is converted into the running mass. Extensional stiffness of the wing is not provided and instead is estimated.

The engine-pylon mass is not available and therefore is not modeled. This initial model is aimed at understanding the geometric nonlinear effect of the tension stiffening, and is not intended to provide an accurate representation of the TBW model. Therefore, the results are used only for the purpose of illustrating the geometric nonlinear tension stiffness, but are not to be compared to the established TBW data.

Aerodynamic loads are computed by an aerodynamic analysis code. Figures 5 and 6 are the plots of lift and pitching moment distributions for the wing and the strut at the design cruise C = 0 . 81 at Mach 0.7 and an altitude of L 42,000 ft, corresponding to the 1g flight load or a cruise gross weight of 170,250 lbs.

Lift Distribution at 1g Pitching Moment at 1g 1400 50 Wing Wing Strut Strut −50 −100 −150 −200 Lift Distribution, lb/ft −250 −300 Quarer−Chord Pitching Moment Distribution, lb 0 −350 0 20 40 60 80 100 0 20 40 60 80 100 BBL, ft BBL, ft Fig. 5 - Lift Distribution Fig. 6 - Pitching Moment Distribution A nonlinear static analysis is performed to include the geometric nonlinear effect of tension stiffening. The solution is initialized with zero initial tension. After the first iteration, the tension is computed and the stiffness matrix is updated. The solution is iterated until the wing tip solution converges. The wing tip deflection from the linear solution with no tension stiffening is 1.5434 ft which indicates that the combined wing-strut structure is quite tiff. The nonlinear solution predicts a wing tip deflection of 1.5248 ft, which only reduces by about 1%. So the tension stiffening effect does not appear to be significant at 1g load. For comparison, the nonlinear wing tip deflections at 2.5g and -1g are computed to be 3.0622 ft and -1.8885 ft, respectively. These are compared to the corresponding linear wing tip deflections of 3.1009 ft and -1.9194 ft. The computed wing static deflections are shown in Fig. 7.

21 of 26 American Institute of Aeronautics and Astronautics 1g without Tension Stiffening 1g with Tension Stiffening 2.5g with Tension Stiffening −1g with Tension Stiffening Wing Vetical Deflection, ft −1 −2 0 20 40 60 80 100 BBL, ft Fig. 7 - Wing Static Deflection The axial extension of the strut for 1g, 2.5g, and -1g flight loads is plotted in Fig. 8. The tension carried by the strut is then computed from the axial strain and is plotted in Fig. 9. The static aeroelastic analysis shows that the strut on the average carries about 153,000 lbs to support a take-off gross weight of 170,250 lbs.

x 10 0.1 4 1g 2.5g 0.08 −1g 3 0.06 0.04 0.02 0 1g 2.5g Strut Tension, lb −0.02 −1g Strut Axial Extension, ft −1 −0.04 −2 −0.06 −0.08 −3 0 10 20 30 40 50 0 5 10 15 20 25 30 35 40 Strut Axial Position, ft Strut Axial Position, ft Fig. 8 - Axial Extension of Strut Fig. 9 - Tension in Strut An eigenvalue analysis is performed to compute the no-load structural dynamic natural frequencies of the TBW and the structural dynamic natural frequencies with tension stiffening due to the loaded structure at 1g, 2.5g, and -1g loads. The natural frequencies of the first ten modes are shown in Table 2. It can be seen that the natural frequencies with tension stiffening at 1g and 2.5g are larger than the no-load natural frequencies. The increase in the natural 22 of 26 American Institute of Aeronautics and Astronautics frequencies from the no-load natural frequencies ranges from about 1% to 8% for the 1g flight load case and from about 1% to 15% for the 2.5g flight load case. The -1g flight load case exhibits a significant reduction in the natural frequencies ranging from about 1% to 24%, particularly at higher modes.

No-Load Frequency Frequency Frequency Mode Frequency (Hz) at 1g (Hz) at 2.5g (Hz) at-1g (Hz) 1 1.9413 2.0052 2.0718 1.8062 2 2.3662 2.4541 2.5483 2.2090 3 4.0659 4.1832 4.3246 3.8891 4 6.6962 6.7525 6.8226 6.6082 5 8.5328 8.6784 8.8275 8.2343 6 10.0313 10.0461 10.0611 9.9881 7 11.4553 11.4804 11.5082 11.4106 8 17.1444 17.2365 17.3415 15.5022 9 18.2639 19.8070 21.3600 16.6709 10 22.1617 22.2387 22.5115 17.9237 Table 2 - Natural Frequencies of TBW Aircraft A flutter analysis is performed for the no-load case and the three load cases. A 2% structural damping is assumed.

The no-load case shows a flutter occurring at modes 4 and 5 as these two modes coalesce together, as shown in Fig.

10. The flutter speeds are 589.00 knots equivalent airspeed (KEAS) for mode 4 and 574.85 KEAS for mode 5, as shown in Fig. 11. Figures 12 to 17 are the plots of natural frequencies and damping of the TBW aircraft for the three flight loads at 1g, 2.5g, and -1g. The flutter speeds for the three flight load cases are shown in Table 3. As can be seen, the flutter speeds are dependent on the flight loads due to the geometric nonlinear effect of tension stiffening.

The difference in the flutter speed is as much as 35 KEAS between the 1g flight load and -1g flight load and 15 KEAS between the 1g flight load and 2.5g flight load.

No−Load Damping No−Load Frequency 0.2 8 0 7 −0.2 −0.4 −0.6 Damping (G) 4 Frequency (Hz) −0.8 3 −1 −1.2 −1.4 0 100 200 300 400 500 600 700 0 100 200 300 400 500 600 700 Airspeed (KEAS) Airspeed (KEAS) Fig. 10 - No-Load Frequency Fig. 11 - No-Load Damping 23 of 26 American Institute of Aeronautics and Astronautics Frequency at 1g Damping at 1g 9 0.2 8 0 7 −0.2 6 −0.4 5 −0.6 Damping (G) 4 Frequency (Hz) −0.8 3 −1 2 −1.2 1 −1.4 0 100 200 300 400 500 600 700 0 100 200 300 400 500 600 700 Airspeed (KEAS) Airspeed (KEAS) Fig. 142 - Frequency at 1g Fig. 13 - Damping at 1g Frequency at 2.5g Damping at 2.5g 9 0.2 8 0 7 −0.2 6 −0.4 5 −0.6 Damping (G) 4 Frequency (Hz) −0.8 3 −1 2 −1.2 1 −1.4 0 100 200 300 400 500 600 700 0 100 200 300 400 500 600 700 Airspeed (KEAS) Airspeed (KEAS) Fig. 14 - Frequency at 2.5g Fig. 15 - Damping at 2.5g 24 of 26 American Institute of Aeronautics and Astronautics Frequency at −1g Damping at −1g 9 0.2 8 0 7 −0.2 6 −0.4 5 −0.6 Damping (G) 4 Frequency (Hz) −0.8 3 −1 2 −1.2 1 −1.4 0 100 200 300 400 500 600 700 0 100 200 300 400 500 600 700 Airspeed (KEAS) Airspeed (KEAS) Fig. 16 - Frequency at -1g Fig. 17 - Damping at -1g No-Load Flutter Flutter Speed Flutter Speed Flutter Speed Mode Speed (KEAS) at 1g (KEAS) at 2.5g (KEAS) at-1g (KEAS) 4 589.00 601.15 616.54 566.48 5 574.85 585.62 596.83 555.19 Table 3 - Flutter Speed of TBW Aircraft The flutter Mach numbers at an altitude of 42,000 ft are shown in Table 4. The Boeing Report indicates a dive speed of Mach 0.82. A flutter margin of 15% results in a minimum clearance flutter Mach 0.9430. Thus, the flutter Mach numbers are well above the clearance flutter Mach. It should be noted that this analysis does not include the engine-pylon mass an well as the correct fuel mass. So the computed flutter speeds in this study are optimistic.

No-Load Flutter Flutter Mach Flutter Mach Flutter Mach Mode Mach at 1g at 2.5g at-1g 4 1.2867 1.3133 1.3469 1.2375 5 1.2558 1.2793 1.3038 1.2129 Table 4 - Flutter Mach Number of TBW Aircraft at Altitude of 42,000 ft In general, increased wing loading results in an increase in the flutter speed. This is consistent with the observation from the wind tunnel test of the TBW in NASA Langley’s Transonic Dynamic Tunnel (TDT) which shows that the flutter speed at a negative angle of attack is consistently lower than that at a positive angle of attack at the same 3, 10 dynamic pressure.

VII. Conclusion This paper presents a study of the geometric nonlinear effect of axial loading on the aeroelasticity of a truss-braced wing aircraft structure. The main struts which restrain the high-aspect ratio wings generally carry a tensile force due 25 of 26 American Institute of Aeronautics and Astronautics to the lift force generated by the wings. The tension causes a geometric nonlinear stiffening effect which effectively raises the stiffness of the structure. This geometric nonlinear effect therefore causes the natural frequencies to increase, hence changes to flutter speed solutions. A nonlinear finite-element model of the trussed-braced wing aircraft has been developed. The nonlinear static deflection analysis is performed and yields a slightly lower wing deflection than the linear analysis. Eigenvalue analysis is performed and shows that the natural frequencies at the 1g and 2.5g flight loads are higher than the no-load natural frequencies. Conversely, the natural frequencies at the -1g flight load is lower than the no-load natural frequencies. A flutter analysis is performed and shows that two flutter modes exist. Modes 4 and 5 flutter at different airspeeds depending on the three flight loads. This is caused by the geometric nonlinear effect of tension stiffening due to the strut which affects the flutter frequencies. In general, increased wing loading results in an increase in the flutter speed. The study illustrates the importance of accounting for the geometric nonlinear tension stiffening effect in analyzing the truss-braced wing aircraft.

Acknowledgment The authors would like to thank the Fixed Wing / Advanced Air Transport Technology Project under the Funda- mental Aeronautics Program of NASA Aeronautics Research Mission Directorate (ARMD) for funding support of this work. The authors also would like to acknowledge Boeing Research and Technology for providing the Truss-Braced Wing aircraft models.

References Bradley, M. K. and Droney, C. K., “Subsonic Ultra Green Aircraft Research: Truss Braced Wing Design Exploration,” Contractor Report, The Boeing Company, June 2014.

Bradley, M. K., Droney, C. K., and Allen, T. J., “Subsonic Ultra Green Aircraft Research: Truss Braced Wing Aeroelastic Test Report,” Contractor Report, The Boeing Company, June 2014.

Bartels, E. R., Scott, R. C., Funk, C., J., Allen, T., J., and Sexton, B. W., “Comparisons of Computed and Experimental Aeroelastic Stability for the Boeing Truss-Braced Wing Wind Tunnel Model,” 32nd AIAA Applied Aerodynamics Conference, Atlanta, June 2014.

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

Nguyen, N., “Integrated Flight Dynamics Modeling of Flexible Aircraft with Inertial Force-Propulsion - Aeroelastic Couplings,” 46th AIAA Aerospace Sciences Meeting and Exhibit, AIAA-2008-0194, 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.

Nguyen, N., Ting, E., Nguyen, D., Trinh, K., “Flutter Analysis of Mission-Adaptive Wing with Variable Camber Continuous Trailing Edge Flap,” AIAA Science and Technology Forum, AIAA-2014-0839, January 2014.

Theodorsen, T., “General Theory of Aerodynamic Instability and the Mechanism of Flutter,” NACA Report 496, 1949.

Nguyen, N., Guist, R., and Muzzio, D., “Experimental Investigation of the Rotor Blade Vibration in the Three-Stage Compressor of the 11-By 11-Foot Transonic Wind Tunnel,” 31st AIAA/ASME/SAE/ASEE Joint Propulsion Conference, AIAA-1995-3139, July 1995.

Bradley, M. K., Droney, C. K., and Allen, T. J., “Subsonic Ultra Green Aircraft Research Phase II - Truss Braced Wing Design Exploration” Contractor Report, The Boeing Company, June 2014.

26 of 26 American Institute of Aeronautics and Astronautics