Flutter Analysis of the Transonic Truss-Braced Wing Aircraft Using Transonic Correction

Flutter Analysis of the Transonic Truss-Braced Wing Aircraft

Using Transonic Correction

∗ Nhan Nguyen NASA Ames Research Center, Moffett Field, CA 94035 † Jason Fugate Stinger Ghaffarian Technologies, Inc., Moffett Field, CA 94035 ‡ Juntao Xiong Stinger Ghaffarian Technologies, Inc., Moffett Field, CA 94035 § Upender Kaul NASA Ames Research Center, Moffett Field, CA 94035 This paper describes a flutter analysis method for the Transonic Truss-Braced Wing aircraft using a vortex- lattice method coupled to an unsteady transonic correction method to account for unsteady aerodynamics in transonic flow. A steady-state vortex-lattice model of the Truss-Braced Wing aircraft is developed using vortex- lattice code VSPAERO. A transonic and viscous flow correction method is implemented in the VSPAERO model to account for steady-state transonic and viscous flow effects using transonic small disturbance code TSFOIL coupled to an in-house integral boundary layer code. In addition, a wing-strut interference correction method is developed to account for the transonic interference aerodynamics in the strut juncture region using high-fidelity CFD code FUN3D. A structural dynamic finite-element model of the Truss-Braced Wing aircraft is developed using BEAM3D in-house finite-element code and is coupled to the VSPAERO. The BEAM3D model includes a geometric nonlinearity due to the tension in the strut which causes a deflection-dependent nonlinear stiffness. An unsteady transonic correction method is developed to better capture the unsteady aerodynamics in transonic flow. The unsteady transonic correction method makes use of the Theodorsen’s theory to account for the amplitude and phase shift of the unsteady lift coefficient in transonic flow. A preliminary flutter analysis of the Truss-Braced Wing aircraft is conduct to illustrate the unsteady transonic correction approach.

I. Introduction The demand for green aviation is expected to increase with the need for reduced environmental impact of air travel.

Most large transports today operate within the best cruise L / D range of 18-20 using the conventional tube-and-wing design. This configuration has led to incremental improvements in aerodynamic efficiency over this past century. A big opportunity has been shown in recent years to significantly reduce structural weight or trim drag, hence improved energy efficiency, with the use of lightweight materials such as composites. The Boeing 787 transport is an example of a modern airframe design that employs lightweight structures. High aspect ratio wing design can provide another opportunity for further improvements in energy efficiency.

Research and development of high aspect ratio wing transport designs has placed a greater emphasis on the studies of aeroelasticity and flutter owing to the increase in the wing flexibility as the wing aspect ratio increases. These studies have sought to develop methods and tools for aeroelasticity by laying the foundation for modern high aspect ratio wing 1–3 aircraft such as the Transonic Truss-Braced Wing (TTBW). The concept of using truss structures to alleviate the wing root bending moment of an ultra-high aspect ratio wing has a long history in aviation with more than a decade of work focused on improving the understanding of the aeroelastic properties and structural weight penalties associated with the increase in the wing flexibility.

∗ Intelligent Systems Division, Technical Group Lead, Senior Research Scientist, Associate Fellow, nhan.t.nguyen@nasa.gov † Intelligent Systems Division, Aerospace Engineer, jason.fugate@nasa.gov ‡ Intelligent Systems Division, Research Engineer, juntao.xiong@nasa.gov § Computational Aerosciences Branch, NASA Advanced Supercomputing (NAS) Division, Senior Research Scientist, Associate Fellow, upen- der.kaul@nasa.gov 1 of 25 American Institute of Aeronautics and Astronautics The Subsonic Ultra Green Aircraft Research (SUGAR) Transonic Truss-Braced Wing (TTBW) aircraft concept is a Boeing-developed N+3 aircraft configuration funded by the NASA ARMD Advanced Air Transport Technologies 4, 5 (AATT) project. The TTBW aircraft concept is designed to be aerodynamically efficient by employing a wing aspect ratio of about 19.55 which is significantly greater than those of cantilever wing transport configurations. For example, the latest Boeing 777-X is reported to have a wing aspect ratio of 11. Without structural bracing, the increase in the wing root bending moment would require a significant structural reinforcement which would lead to an increase in the structural weight that would offset the aerodynamic benefit of the high aspect ratio wing. Thus, the design of a truss-braced structure is a multidisciplinary design optimization process that strives to achieve a delicate balance between aerodynamic efficiency and structural efficiency. In the SUGAR configuration, the trade between the aerodynamic performance and structural design results in a truss-braced configuration with the wings braced at approximately mid-span by two main struts. In addition, two jury struts, one on each wing, provide the additional structural reinforcement. Two versions of the TTBW configurations are currently being developed by Boeing; a Mach 0.745 version and a Mach 0.8 version. Figure 1 is an illustration of the Mach 0.745 TTBW aircraft.

Figure 1. Boeing SUGAR Transonic Truss-Braced Wing (TTBW) Aircraft Concept The development of the TTBW aircraft is funded by the NASA AATT project through a NRA (NASA Research Announcement) selection of Boeing Research and Technology as the lead organization. Multidisciplinary design analysis and optimization (MDAO) studies have been conducted at each stage to improve the wing aerodynamics, structural efficiency, and flight performance using advanced N+4 turbofan engines. These MDAO studies have refined the geometry of the wing and configuration layout and have involved trade studies involving minimizing induced drag, profile drag, and wave drag due to the addition of the main strut and jury struts.

Owing to the high aspect ratio flexible wings, significant aeroelastic interactions are expected for the TTBW aircraft. Aeroelastic interactions can result in adverse aerodynamics that can lead to drag penalty at off-design per- formance. The TTBW, therefore, could be a prime candidate for the performance adaptive aeroelastic wing (PAAW) 6–8 technologies such as the Variable Camber Continuous Trailing Edge Flap (VCCTEF) system currently developed under the PAAW research element of the AATT project. The PAAW technologies provide a mission-adaptive wing capability to improve aerodynamic efficiency and structural efficiency of a transport aircraft wing by active controls to minimize drag and to suppress structural responses. A MDAO framework is currently being developed for analyz- ing high aspect ratio flexible wing aircraft configurations such as the TTBW and for the applications of the PAAW technologies including the VCCTEF. This framework is intended to include a suite of aerodynamic tools of varying fidelity as well as finite-element models of these configurations in order to address aero-structural, aeroelasticity, and aeroservoelasticity aspects of the TTBW development. Using this framework, a maneuver load alleviation capability of the VCCTEF for the TTBW has been investigated with the flap layout as shown in Fig. 2. Thus, options may exist to further improve structural efficiency using the PAAW technologies.

While there are on-going efforts in developing high-fidelity Reynolds-averaged Navier-Stokes (RANS) compu- tational fluid dynamics (CFD) models of the TTBW for aerodynamic design and performance analysis, these high- fidelity CFD models are typically not coupled to the structure of the TTBW. There is a need to conduct aeroelastic 2 of 25 American Institute of Aeronautics and Astronautics analysis of the TTBW configuration to address flutter requirements and aero-structural coupling which can affect aerodynamics and vehicle stability and control. Thus, structural coupling is an important element of the TTBW devel- opment. To enable this analysis capability, one approach is to use a lower-order aerodynamic tool that can still capture sufficiently the transonic and viscous flow effects on the vehicle for coupling to a finite-element model of the TTBW structure. This approach can provide a rapid analysis capability that blends itself to vehicle optimization with aeroe- lastic constraints such as load and flutter. This is a practical approach considering that currently a vehicle optimization with flutter constraints cannot be handled by high-fidelity CFD tools due to the extremely high computational cost of running unsteady CFD simulations coupled with finite-element models. The industry standard for flutter analysis is the doublet lattice method available in NASTRAN, but this tool does not account for the transonic and viscous flow effects. Corrections to the section lift curve slope can be made in NASTRAN by using steady state RANS CFD solutions. During a design optimization, as the vehicle configuration changes, new steady state RANS CFD models would have to be created. Thus, even this approach can still present a large computational cost.

Figure 2. TTBW with Notional Layout of VCCTEF System The TTBW configuration provides a structural solution to high aspect ratio wing aircraft designs. The long slender wing employs 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 span load carried by the wing and are generally loaded in tension. Under a negative-g flight condition such as during a dive, a load reversal could occur that could put the struts in compression. The compressive loading would require a design consideration 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 presence of axial loading causes the bending and torsional stiffnesses to change. A tensile loading will result in an increase in bending and torsional stiffnesses. This tension-induced stiffening effect is a geometric nonlinearity which should be captured in aeroelastic analyses of the TTBW configuration. Thus, the TTBW structure is inherently a nonlinear structure. A recent study has investigated the effect of axial loading of the TTBW configuration on the flutter speed.

The study shows that the flutter speed has a nonlinear characteristic and is sensitive to the load factor or the angle of attack.

This paper presents a flutter analysis approach for the TTBW aircraft utilizing the potential flow solver VSPAERO coupled to an in-house nonlinear finite-element code BEAM3D. The VSPAERO models of the TTBW include both a low-fidelity vortex-lattice model and a mid-fidelity panel model for steady-state aerodynamics. Transonic and viscous flow corrections for the steady-state aerodynamics are implemented on the vortex-lattice model using a transonic small disturbance (TSD) code called TSFOIL coupled to an in-house integral boundary layer (IBL) code. In the region near the strut attachment to the wing, the flow involves a considerable degree of interactions between the wing and the strut. A high-fidelity CFD model of the TTBW is developed using FUN3D to investigate the wing-strut interference aerodynamics for the purpose of developing a wing-strut interference aerodynamic correction method to be applied to the VSPAERO model. Upon applying all the necessary corrections, the VSPAERO model shows an excellent agreement with wind tunnel test data for the Mach 0.745 version of the TTBW.

The BEAM3D finite-element model of the TTBW is based on a three-dimensional beam theory which includes all six degrees of freedom at a node. The BEAM3D model is coupled to the VSPAERO for static aeroelastic analysis.

The mass and stiffness properties of the TTBW are provided by Boeing. The BEAM3D model has been validated against NASTRAN static deflection analysis and frequency analysis. Within BEAM3D, unsteady aerodynamic mass, damping, and stiffness matrices are implemented using the Theodorsen’s theory for unsteady aerodynamics and the 3 of 25 American Institute of Aeronautics and Astronautics

steady-state lift curve slope c for steady-state aerodynamics computed by VSPAERO. An unsteady transonic correc- l α tion method to the Theodorsen’s theory is proposed in this paper to better capture the changes in the amplitude and phase shift in transonic flow. This modeling approach can provide a rapid aeroelastic analysis capability for flutter analysis without incurring the high computational overhead associated with a RANS CFD-based approach.

II. Aerodynamic Modeling of the Truss-Braced Wing

A. VSPAERO Models In order to develop a rapid aeroelastic capability for flutter analysis that blends itself in a vehicle MDAO process, a lower-fidelity aerodynamic model of the TTBW is necessary. A previous effort in modeling the TTBW using the vortex-lattice code VORLAX based on the aerodynamic superposition principle has shown that low-order aerody- namic tools can provide a reasonably accurate prediction of the aerodynamic performance of the TTBW. Due to the limitation of VORLAX in the previous effort, an improved vortex-lattice model as well as a panel model of the TTBW has been developed in VSPAERO. VSPAERO is a solver that includes both the vortex lattice method and the full panel method based on generalized vortex loops. The core VSPAERO solver is based on an agglomerated multi-pole approach, coupled with a preconditioned linear solver, to reduce solution times. Adaptive wakes, time-accurate un- steady analyses, and propeller modeling are all supported. The latest version of VSPAERO supports a loosely coupled integration of BEAM3D finite-element code to perform static aeroelastic analyses. VSPAERO is part of the OpenVSP design package and is freely available under the NASA open source license.

Figure 3. VSPAERO Model of TTBW Figure 4. Mesh of VSPAERO Vortex-Lattice Model of TTBW Two Mach 0.745 TTBW OML (Outer Mold Line) geometries are available: cruise or 1g shape geometry and jig shape geometry. The geometries for the Mach 0.745 TTBW are based on the CAD models of the 765-095-200_RJ version. The jig shape geometry is not a flight jig shape that accounts for the wing aeroelastic wash-out twist for the full-scale vehicle at the design flight condition, but is a wind tunnel model jig shape. Wind tunnel tests of the 1g shape geometry have been conducted in NASA Ames 11-Ft Transonic Wind Tunnel. Both geometries are modeled 4 of 25 American Institute of Aeronautics and Astronautics in VPSAERO. Figure 3 illustrates one of the TTBW geometries used in the VSPAERO models. The mesh for the vortex-lattice model of the TTBW in VSPAERO is shown in Fig. 4. The surface triangulation mesh of the panel model is shown in Fig. 5.

Figure 5. Surface Triangulation of VSPAERO Panel Model of TTBW ◦ Figure 6 shows the differential pressure coefficient contour at Mach 0.7 and the angle of attack of 4 for the VSPAERO vortex-lattice model. The shed wakes are shown in Fig. 6 to illustrate the time-stepping wake capability in VSPAERO. Figure 7 shows the surface pressure coefficient contour computed by the VSPAERO panel model.

Figure 6. Differential Pressure Coefficient Contour of VSPAERO Panel Model of TTBW Figure 7. Surface Pressure Coefficient Contour of VSPAERO Panel Model of TTBW B. Transonic and Viscous Flow Corrections Because of the missing transonic effect in the linear potential flow method, a method for transonic and viscous flow 12–14 corrections has recently been developed. In this method, a full-configuration aerodynamic model can be based 5 of 25 American Institute of Aeronautics and Astronautics on a potential flow method such as the vortex-lattice method or the panel method. The wing is discretized into several spanwise sections at which the section lift coefficients computed by the potential flow method are used to correct for 13, 14 the transonic and viscous flow effects. The transonic flow correction is handled by a transonic small disturbance (TSD) code called TSFOIL. This code is loosely coupled to an in-house integral boundary layer (IBL) code to correct for the viscous flow interaction with the transonic shock on an airfoil. The correction method is an iterative process to compute the incremental section lift coefficient due to transonic and viscous flow by a virtual re-twist of the individual wing sections to account for the accompanied change in the effective local angle of attack to produce the incremental section lift coefficient. Optionally, the correction method is also developed using the 2D Euler CFD code MSES with an integral boundary layer method developed by Mark Drela. The transonic and viscous flow correction method is shown in Fig. 8. An extensive validation of the transonic and viscous flow correction method has been performed to compare the method against the OVERFLOW results for a 2D airfoil, and when used with a vortex-lattice model against LAVA, a RANS CFD code developed at NASA Ames, and the Euler CFD code CART3D without the viscous flow correction for a full configuration aircraft. The agreement in the L / D result for the NASA Generic Transport Model (GTM) is less than 4% of the LAVA-computed value.

Figure 8. Transonic and Viscous Flow Correction Flow Chart Rigid GTM (M = 0.797) Rigid GTM (M = 0.797) 0.65 0.03 RANS (LAVA) RANS (LAVA) 0.028 VLM+TSD/IBL VLM+TSD/IBL 0.6 VLM+MSES VLM+MSES 0.026 0.55 0.024 0.5 0.022 L D 0.45 0.02 C C 0.018 0.4 0.016 0.35 0.014 0.3 0.012 0.25 0.01 0 0.5 1 1.5 2 2.5 3 0.2 0.25 0.3 0.35 0.4 0.45 0.5 0.55 , deg C L Figure 9. Lift Curve and Drag Polar for the Rigid Generic Transport Model Computed by TSFOIL/Integral Boundary Layer and MSES Coupled to VORLAX Compared with LAVA RANS Results 6 of 25 American Institute of Aeronautics and Astronautics Figure 9 shows the lift and drag coefficients of the GTM computed by VORLAX coupled to the transonic and viscous flow corrections and MSES. The major advantage is the computational efficiency of the method which is several order of magnitude faster than RANS CFD. This computational efficiency becomes highly important when the potential flow solver is coupled to a structural finite-element model for aero-structural modeling or flutter analysis which adds an additional computational overhead to the aerodynamic modeling.

C. Wing-Strut Interference Aerodynamic Correction The transonic and viscous flow corrections using the TSD/IBL method are generally valid for a single-element airfoil.

The TTBW configuration is a complex geometry that includes a strut juncture region where the effect of interference aerodynamics can influence the overall aerodynamic performance of the aircraft. As the strut approaches the wing from below, the transonic and viscous flow corrections using the TSD/IBL method are no longer valid due to the interactions between the wing airfoil and strut airfoil. To account for the interference aerodynamics between the wing and the strut, CFD models of the wing-strut configuration and the wing-alone configuration of the jig shape TTBW geometry are developed using FUN3D, as shown in Fig. 10. To isolate the interference aerodynamic effect for the wings and struts alone, the engines and the horizontal tails are removed from the models.

The Roe’s flux-difference splitting scheme and the Spalart-Allmaras turbulence model are used in the study. Figure 11 shows the grids for the full configuration with the engine, which is used for the grid independent study and for validation. The total number of nodes is 90 million. In this study, the mesh is comprised of tetrahedral elements and prism layer near the wall. The prism layer is used to resolve the turbulent boundary layer.

Figure 10. Jig Shape TTBW Wing-Strut Configuration and Wing-Alone Configuration in FUN3D Figure 11. TTBW CFD Mesh Surface pressure coefficients are computed for both configurations are various wing stations at Mach 0.3, 0.745, 0.8, and 0.85. Figure 12 shows the surface pressure coefficients at the wing station immediately inboard from the strut ◦ juncture at an angle of attack of 2 . The interference aerodynamic effect is seen to become increasingly significant as the Mach number increases and is most pronounced on the lower surface. A high speed flow region on the lower surface is indicated by a sharp decrease in the surface pressure coefficient especially at Mach 0.8 and Mach 0.85. This results in a decrease in the section lift coefficient.

A wing-strut interference correction method is developed to correct the VSPAERO model coupled to the TSD/IBL approach. The correction method is applied to the VSPAERO+TSD/IBL model to update the section lift, drag, and pitching moment coefficient as well as the lift curve slope c and the location of the aerodynamic center x / c of each l ac α wing section which are two important parameters for flutter analysis. Figure 13 shows the incremental values of c l α and x / c along the wing span at Mach 0.745 between the FUN3D results and the VSPAERO+TSD/IBL results. As ac 7 of 25 American Institute of Aeronautics and Astronautics can be seen, near the strut juncture region, there is a significant decrease in the value of c . As the flow becomes l α increasingly more transonic, the location of the aerodynamic center tends to gradually increase from the quarter-chord location ( x / c = 0 . 25 ) in subsonic flow to near the mid-chord location ( x / c = 0 . 5 ) when the flow becomes entirely ac ac supersonic. The incremental value of the location of the aerodynamic center is generally small and is less than 1%.

The wing-strut interference correction is generally small at Mach 0.745, but can become larger as the Mach number increases.

◦ Figure 12. Pressure Coefficients at Wing Station as Function of Mach Number at α = 2 6 6 TTBW Interference Wing c (M = 0.745, Re = 2.6x10 ) TTBW Interference Wing x /c (M = 0.745, Re = 2.6x10 ) l ac 0.01 0.2 0.008 0.1 0.006 0.004 -0.1 /c l ac c x 0.002 -0.2 -0.3 -0.4 -0.002 -0.5 -0.004 0 200 400 600 800 1000 1200 0 200 400 600 800 1000 1200 Wing Station, in Wing Station, in Figure 13. Incremental Lift Curve Slope and Aerodynamic Center Location at Mach 0.745 8 of 25 American Institute of Aeronautics and Astronautics D. Aerodynamic Analysis of Cruise Shape Geometry The VSPAERO vortex-lattice model is used for the aerodynamic analysis of the cruise shape TTBW geometry. The VSPAERO panel model is not used since it has not reached the same level of maturity as the VSPAERO vortex-lattice model with capabilities to couple to the BEAM3D finite-element m odel and the TSD/IBL method still not yet available. In addition, a CFD model of the full-configuration cruise shape TTBW geometry is developed in FUN3D. Wind tunnel test data of the cruise shape geometry in NASA Ames 11-Ft Transonic Wind Tunnel are available for validation of the VSPAERO and FUN3D models. Test data from Run 380 at Mach 0.745 and Reynolds number of 2.6 million based on the mean aerodynamic chord (MAC) with full wind tunnel model corrections are used for validation.

Figure 14 shows the unlabeled plots of the lift and drag coefficients computed by VSPAERO for Mach 0.745 and a Reynolds number of 2.6 million with and without all the corrections as well as FUN3D. These computed results are compared to Run 380 wind tunnel data. While the lift coefficient is somewhat under-predicted, the drag coefficient computed by VSPAERO alone is mostly induced drag and therefore is substantially under-predicted. With all the corrections applied to the VSPAERO model for transonic viscous flow and wing-strut interference aerodynamics, the lift and drag coefficients match remarkably well the wind tunnel data, although there is a small discrepancy in the drag polar at lower lift coefficients. For comparison, the FUN3D results are also plotted. The lift coefficient computed by FUN3D is in less agreement with the wind tunnel data and is higher than that computed by VSPAERO with all the corrections applied. The drag polar computed by FUN3D shows an excellent agreement with the wind tunnel data.

6 6 TTBW Lift Curve (M = 0.745, Re = 2.6x10 ) TTBW Drag Polar (M = 0.745, Re = 2.6x10 ) VSPAERO (1g) VSPAERO+TSD/IBL+Interference (1g) Boeing 11-Ft Wind Tunnel Data (Run 380) FUN3D L L C C VSPAERO (1g) VSPAERO+TSD/IBL+Interference (1g) Boeing 11-Ft Wind Tunnel Data (Run 380) FUN3D -3 -2 -1 0 1 2 3 4 5 C D , deg Figure 14. TTBW Lift Curve C vs. α and Drag Polar C vs. C at Mach 0.745 L L D Figure 15. TTBW Pressure Coefficient Contour Computed by FUN3D at Mach 0.745 and C = 0 . 73 L Figure 15 is the plot of the pressure coefficient contour computed by FUN3D at Mach 0.745 and the design lift coefficient C = 0 . 73. Figure 16 shows the plot of the computed L / D ratios as a function of the lift coefficient from L VSPAERO, FUN3D, and wind tunnel data. All the three curves of L / D ratios are in very good agreement.

The excellent agreement between the VSPAERO model and the wind tunnel data thus has validated the aerody- 9 of 25 American Institute of Aeronautics and Astronautics

namic modeling approach. The VSPAERO model is then used to extract the lift curve slope c for the wings and the l α struts for flutter analysis.

TTBW L/D (M = 0.745, Re = 2.6x10 ) L/D VSPAERO+TSD/IBL+Interference (1g) FUN3D Boeing 11-Ft Wind Tunnel Data (Run 380) C L Figure 16. TTBW L / D at Mach 0.745

III. Finite-Element Model of the Truss-Braced Wing

A. BEAM3D Model BEAM3D is an in-house finite-element code that is used to model a general wing structure using beam, rod, and rigid elements. The code contains a number of analysis modules: a static deflection analysis module, an eigenvalue analysis module, a flutter analysis module, and an aeroservoelastic flight dynamic simulation module. The static de- flection analysis module can be coupled to a flow solver such as VORLAX, VSPAERO, and CART3D to perform aero-structural and static aeroelastic deflection analyses. The eigenvalue analysis module computes the structural dynamic frequencies and the normal mode shapes of a wing structure. The flutter analysis module computes the fre- quencies and damping of the wing aeroelastic modes and performs a flutter speed search in the frequency domain using the Theodorsen’s theory for unsteady aerodynamics or in the time domain using the R. T. Jones approximation to gen- erate the aerodynamic mass, damping, and stiffness. The aeroservoelastic flight dynamic simulation module provides the capability to couple the dynamic aeroservoelastic state-space model of a wing structure to a nonlinear 6-degree-of- freedom flight dynamic model of an aircraft to perform time simulations of the aircraft response to atmospheric gust inputs and flight control surface inputs.

Figure 17 illustrates the TTBW finite-element model in BEAM3D. The left wing and main strut are modeled as the major structural elements. The jury strut is neglected in the model. The mass and structural properties of the finite- element model are made available to NASA by Boeing. The stiffness properties and the elastic definition of the wing and strut are included in the structural properties. BEAM3D models the mass and inertia properties as distributed or running mass and inertias. So, the point masses and inertias from the NASTRAN model provided by Boeing have to be converted into the equivalent distributed mass and inertias. Fuel weight is also modeled in the BEAM3D model.

The full fuel weight corresponds to a take-off gross weight (TOGW) of 138,000 lbs. The inertias of the fuselage and tails are estimated in OpenVSP.

BEAM3D includes a nonlinear capability to model the tension-induced stiffening effect of the struts which con- tributes to the overall stiffness matrix of the finite-element model. This is expressed as s a t K ( u ) = K + K + K ( u ) (1) i i i i s a where K is the total element stiffness matrix of the i -th element, K is the element structural stiffness matrix, K is i i i t the element aerodynamic stiffness matrix, and K is the element nonlinear stiffness matrix due to the strut tension i which is dependent on the displacement field u ( t ) of the structure. The aerodynamic stiffness is computed from the 10 of 25 American Institute of Aeronautics and Astronautics Theodorsen’s theory and depends on the dynamic pressure, wing section lift curve slope, and other pertinent unsteady aerodynamic parameters.

Figure 17. BEAM3D Finite Element Model of TTBW B. Static Aeroelastic Analysis For static aeroelastic analysis, an iterative solution method is implemented for the nonlinear tension-induced stiffness as shown in Fig. 18. At each iteration between the VSPAERO model and BEAM3D model, the iterative solution method iterates on the static deflection solution until it converges.

A static aeroelastic analysis is performed for the jig shape TTBW geometry. The VSPAERO model of the jig shape TTBW geometry with all the corrections is coupled to the BEAM3D model. Figure 19 shows the linear and nonlinear static aeroelastic deflections of the wings at Mach 0.745 and the design lift coefficient C = 0 . 73 computed L by BEAM3D. For comparison, the linear static aeroelastic deflections are also computed in the Boeing NASTRAN 3D model of the TTBW using the aerodynamic loads computed from the VSPAERO model upon convergence. Except for a minor discrepancy near the wing root, the linear static solutions of the wing bending deflection computed by NASTRAN and BEAM3D otherwise are in excellent agreement. The wing tip deflection at cruise is computed to be 23.3 inches by BEAM3D and 23.2 inches by NASTRAN. The nonlinear static solution of the wing tip deflection is computed to be 22.8 inches. Thus, the nonlinear tension-induced stiffening does not significantly alter the static aeroelastic deflections.

Figure 18. Iterative Solution for Nonlinear Finite-Element Model of Strut Tension Stiffening 11 of 25 American Institute of Aeronautics and Astronautics 6 6 TTBW Wing Static Aeroelastic Deflection (M = 0.745, Re = 11.9x10 , C = 0.73) TTBW Wing Static Aeroelastic Deflection (M = 0.745, Re = 11.9x10 , C = 0.73) L L 25 2.5 NASTRAN (Linear) NASTRAN (Linear) BEAM3D (Linear) BEAM3D (Linear) BEAM3D (Nonlinear) 20 BEAM3D (Nonlinear) 1.5 Bending Deflection, in 0.5 Pitch Axis Twist (+ nose down), deg -5 0 0 200 400 600 800 1000 1200 0 200 400 600 800 1000 1200 Wing Station, in Wing Station, in Figure 19. TTBW Wing Bending Deflection and Twist at Mach 0.745 and C = 0 . 73 L The nose-down wing twists about the aircraft pitch axis or the y -axis computed by BEAM3D and NASTRAN are in general agreement, but a discrepancy exists along the wing span with the wing twist computed by NASTRAN ◦ higher than that computed by BEAM3D by at most 0 . 3 outboard of the wing station y = 400 inches but lower than that ◦ computed by BEAM3D by at most 0 . 2 inboard of this wing station. This discrepancy could be due to the application of the distributed lift and pitching moment about the elastic axis in NASTRAN not being precisely on the elastic axis but rather along the nearest nodes. The overall nose-down wing twist at the tip is in better agreement between ◦ ◦ NASTRAN and BEAM3D with the NASTRAN computed value of 2 . 39 and BEAM3D computed values of 2 . 35 ◦ and 2 . 33 corresponding to the linear and nonlinear static solutions, respectively.

Figure 20 shows the strut axial displacement and tension as a function of the wing station at the design lift coeffi- cient C = 0 . 73. The axial displacement of the struts is almost 0.6 inches and the tension carried by the strut is about L 171,000 lbs.

TTBW Strut Axial Displacement (M = 0.745, Re = 11.9x10 , C = 0.73) TTBW Strut Tension (M = 0.745, Re = 11.9x10 , C = 0.73) L L 0.6 184000 182000 0.5 180000 178000 0.4 176000 0.3 174000 172000 Strut Tension, lb 0.2 170000 Strut Axial Displacement, in 168000 0.1 166000 0 164000 100 200 300 400 500 600 100 200 300 400 500 600 Wing Station, in Wing Station, in Figure 20. TTBW Strut Axial Displacement and Tension at Mach 0.745 and C = 0 . 73 L Figure 21 shows the unlabeled plots of the lift and drag coefficients computed by the VSPAERO model of the jig shape TTBW geometry coupled to the BEAM3D model. For comparison, the lift and drag coefficients computed by the VSPAERO model of the cruise shape TTBW geometry and the wind tunnel data are also plotted. The lift and drag coefficients computed by the VSPAERO+BEAM3D model of the jig shape TTBW geometry agree very well with the VSPAERO model of the cruise shape TTBW geometry and the wind tunnel data. The drag coefficient computed by the VSPAERO+BEAM3D model actually agrees slightly better than that computed by the VSPAERO model. The lift curve of the VPSAERO model matches the slope from the wind tunnel data very well, but has a slight offset.

12 of 25 American Institute of Aeronautics and Astronautics 6 6 TTBW Lift Curve (M = 0.745, Re = 2.6x10 ) TTBW Drag Polar (M = 0.745, Re = 2.6x10 ) L L C C VSPAERO+TSD/IBL+Interference (1g) VSPAERO+TSD/IBL+Interference (1g) VSPAERO+BEAM3D+TSD/IBL+Interference (Jig) VSPAERO+BEAM3D+TSD/IBL+Interference (Jig) Boeing 11-Ft Wind Tunnel Data (Run 380) Boeing 11-Ft Wind Tunnel Data (Run 380) -3 -2 -1 0 1 2 3 4 5 C D , deg Figure 21. TTBW Lift Curve C vs. α and Drag Polar C vs. C at Mach 0.745 L L D C. Frequency Analysis The first 12 structural dynamic modes for the TTBW computed by BEAM3D are presented for the 100% fuel condition corresponding to a TOGW of 138,000 lbs. The natural frequencies are computed with and without the tension-induced stiffening effect of the struts, referred to as 1g-loaded frequencies and no-load frequencies, respectively, as shown in Table 1. The BEAM3D code has been validated against NASTRAN frequency analysis. The lowest no-load natural frequency is 1.977 Hz corresponding to the first symmetric wing bending mode. For comparison, the lowest natural frequency computed by NASTRAN is 1.896 Hz for the same mode. The discrepancy could be due to the inconsistency in the mass used in both models. This will be verified later. The tension-induced stiffening effect increases the lowest no-load natural frequency by 2.6% to 2.028 Hz. Figure 22 shows the 12 natural mode shapes corresponding to the first 12 modes.

Mode No-Load Frequency (Hz) 1g-Loaded Frequency (Hz) st 1 Symmetric 1.9769 2.0282 st 1 Antisymmetric 2.0670 2.1066 nd 2 Symmetric 2.4592 2.5331 nd 2 Antisymmetric 2.4873 2.5468 rd 3 Symmetric 2.9034 2.9189 rd 3 Antisymmetric 2.8908 2.9242 th 4 Symmetric 4.4957 4.5194 th 4 Antisymmetric 5.3688 5.3865 th 5 Symmetric 6.7143 6.9530 th 5 Antisymmetric 5.5832 9.1957 th 6 Symmetric 7.0230 8.7344 th 6 Antisymmetric 9.2363 10.5985 Table 1. TTBW Structural Dynamic Natural Frequencies Computed by BEAM3D 13 of 25 American Institute of Aeronautics and Astronautics

Figure 22. TTBW Natural Mode Shapes

IV. Transonic Flutter Analysis

Transonic flutter is a topic of high interest for aircraft design. Doublet lattice methods are frequently used in flut- ter analysis which can predict low-speed flutter when the flow is entirely subsonic at low Mach number. Transonic flutter, on the other hand, is still an actively researched topic in transport aircraft design. Many methods for tran- sonic flutter have been developed ranging from the doublet-lattice method to high-fidelity CFD-based methods and 20, 21 any intermediate-fidelity methods in between. While the transonic flutter prediction can be more accurate with high-fidelity CFD methods than with unsteady potential flow methods, the computational cost is high. In many ap- plications that involve design optimization with flutter constraints, the computational cost associated with high-fidelity CFD presents a barrier. Therefore, computationally efficient methods for transonic flutter prediction continue to be of high interest to the aircraft design community.

Toward this end, a method for extending the classical Theodorsen’s theory of unsteady aerodynamics for incom- pressible flow to transonic flow is proposed in this study. The method proposes a modified Theodorsen’s function to correct for changes in the amplitude and phase shift of the circulatory lift in transonic flow. The amplitude and phase shift corrections are in general a function of the flow parameters such as the Mach number, the reduced frequency, and the instantaneous angle of attack as well the airfoil characteristics such as the thickness and the type of airfoil.

To derive the unsteady corrections, unsteady CFD simulations are performed using FUN3D for selected reduced fre- quencies and Mach numbers for the TTBW airfoil at the mean aerodynamic chord (MAC) wing station. The unsteady transonic correction method is then implemented in BEAM3D for transonic flutter analysis.

The unsteady transonic correction method is discussed further as follows: A. Unsteady Transonic Correction Consider the unsteady motion of an oscillatory pitching and plunging airfoil with an angular velocity ˙ θ and plunging velocity ˙ h where θ is the instantaneous pitch angle, positive nose-up, and h is the vertical displacement of the airfoil, 14 of 25 American Institute of Aeronautics and Astronautics positive downward. This is shown in Fig. 23. The oscillation of the airfoil causes a downwash on the airfoil. This downwash is equal to ˙ ˙ V = V θ + ( x − x ) θ + h (2) z ∞ e where x is the airfoil position along the chord from the leading edge and x is the location of the center of rotation.

e Figure 23. Downwash of Oscillating Airfoil The slope of the camber line is related to the downwash as [ ] ˙ ˙ dz V ( x − x ) θ h z e = − = − θ + + (3) dx V V V ∞ ∞ ∞ Based on thin-airfoil aerodynamic theory for incompressible inviscid flow, the circulatory lift coefficient due to the effective camber change is evaluated by [ ] ∫ ∫ π π ˙ ˙ C ( k ) c dz C ( k ) c ( x − x ) θ h l l e α α c = − ( 1 − cos ϕ ) d ϕ = θ + + ( 1 − cos ϕ ) d ϕ (4) l π dx π V V 0 0 ∞ ∞ This equation is valid for small-amplitude oscillations for which the Kutta condition at the trailing edge is generally assumed to be valid. For incompressible flow, c = 2 π . The function C ( k ) is the well-known complex-valued l α Theodorsen’s function, the exact expression of which is given by ( 2 ) H ( k ) C ( k ) == = F ( k ) + iG ( k ) (5) ( 2 ) ( 2 ) H ( k ) + iH ( k ) 1 0 ( 2 ) ( 2 ) ω c where H ( k ) and H ( k ) are the Hankel functions, and k = is called the reduced frequency.

0 1 2 V ∞ c Using the transformation x − x = ( cos ϕ − cos ϕ ) , this integral can be evaluated as e e [ ( ) ] ˙ ˙ 3 c θ h c = C ( k ) c θ + − x + (6) l l e α 4 V V ∞ ∞ Thus, the effective angle of attack of a pitching and plunging airfoil is equal to ( ) ˙ ˙ 3 c θ h α = θ + − x + (7) e e 4 V V ∞ ∞ The effective angle of attack can be extended to a swept wing section as follows: ∂ w e cos Λ ∂ θ 1 ∂ w c α = θ cos Λ − sin Λ + − (8) e ∂ ¯ y V ∂ t V ∂ t ∞ ∞ where w ( ¯ y , t ) is the wing vertical bending along the wing elastic axis, denoted by ¯ y , and Λ is the sweep angle.

The circulatory lift force acting on the airfoil in incompressible flow is then given by [ ( ) ] ˙ ˙ 3 c θ h l = c q c = C ( k ) c q c θ + − x + (9) c l ∞ l ∞ e α 4 V V ∞ ∞ 15 of 25 American Institute of Aeronautics and Astronautics As the airfoil undergoes a pitching and plunging motion, an opposing inertial force is created due the surrounding air acting on the airfoil. This unsteady lift force is called non-circulatory due to the apparent mass effect. The local angle of attack along the airfoil chord can be expressed as V z α ( x , t ) = = α ( t ) + α ( x , t ) (10) c 1 2 V ∞ ˙ ˙ ( x − x ) θ h e where α = θ + and α = .

1 2 V V ∞ ∞ According to Bisplinghoff et al., the velocity potential on the upper surface is defined as [ ] √ V c α ( x , t ) + α ( x , t ) ∞ 2 2 m φ ( x , t ) = α ( t ) + 1 − ξ (11) 2 2 2 ( x − x ) c m where ξ = and x = .

m c 2 The non-circulatory lift is evaluated as [ ] ∫ ∫ 1 1 √ ∂ ρ V c ∂ 1 1 ∞ ∞ l = ρ c φ d ξ = α + α + α ( x , t ) 1 − ξ d ξ nc ∞ 1 2 2 m ∂ t 2 ∂ t 2 2 − 1 − 1 2 2 [ ] ρ V π c π c ∞ ∞ ˙ ¨ ¨ ˙ = α ( x , t ) = ρ V θ + ( x − x ) θ + h (12) c m ∞ ∞ m e 4 4 Consider a volume of air cylinder with a diameter c that surrounds the airfoil. The non-circulatory lift force, π c therefore, is effectively equal to the apparent mass per unit length of the air volume ρ times the acceleration of the ∞ c c ˙ airfoil V α ( x , t ) acting at the mid-chord point due to the pitching and plunging motion. Let e = x − , e = − x , ∞ c m e m e 4 2 3 c and e = − x . Then, the total lift force is expressed as c e ( ) ˙ ˙ ( ) θ h π c ˙ ¨ ¨ l = l + l = C ( k ) c q c θ + e + + ρ V θ + e θ + h (13) c nc l ∞ c ∞ ∞ m α V V 4 ∞ ∞ In aeroelasticity, the Theodorsen’s airfoil convention is frequently used. The airfoil chord is c = 2 b and the location of the elastic center is at x = ( 1 + a ) b from the leading edge, as shown in Fig. 24. The parameter a is positive if the e elastic center is behind the mid-chord point and is negative if the elastic center is forward of the mid-chord point. The b aerodynamic center for incompressible flow is at the quarter-chord point x = . Thus, the total lift force according ac to Theodorsen is expressed equivalently as [ ( ) ] ˙ ˙ ( ) 1 b θ h ˙ ¨ ¨ l = 2 C ( k ) c q b θ + − a + + ρ π b V θ − ab θ + h (14) l ∞ ∞ ∞ α 2 V V ∞ ∞ Figure 24. Theodorsen’s Airfoil Convention V t ∞ Let τ = be a non-dimensional time variable. Consider an oscillating airfoil in transonic flow with a harmonic c angle of attack ¯ ¯ α = α + α sin ω t = α + α sin 2 k τ (15) 0 0 where ¯ α is the mean angle of attack.

16 of 25 American Institute of Aeronautics and Astronautics For transonic flow, it is hypothesized that the circulatory lift resembles that for incompressible flow but with a different amplitude and phase shift. To capture the changes in the amplitude and phase shift, a modified Theodorsen’s function is now proposed as ( ) ( ) ( ) t t t C k , M , = C k , M , F ( k ) + iC k , M , G ( k ) (16) T F G c c c ( ) t where C k , M , is a modified Theodoren’s function for transonic flow.

T c ( ) ( ) t t The functions C k , M , and C k , M , represent the amplitude and phase shift correction factors, respectively, F G c c and are assumed to be dependent on the Mach number, the reduced frequency, the airfoil thickness, and the type of airfoil. The circulatory lift coefficient in transonic flow is then modified according to [ ( ) ( )] e d α C G ( k ) d α e d α c G c ¯ c = ¯ c + c C F ( k ) α − α + + + l l l F c α c d τ 2 k d τ c d τ [ ( ) ( )] e e c c ¯ = c + c α + c α C F ( k ) sin 2 k τ + 2 k cos 2 k τ + C G ( k ) cos 2 k τ − 2 k sin 2 k τ (17) l l l 0 F G α α c c The quantity ¯ c is the mean unsteady lift coefficient due to the airfoil camber and the mean angle of attack ¯ α . It is l assumed that the inertial force acting on the airfoil is the same in transonic flow as in incompressible flow and is given by ( ) 2 ( ) π d α e d α e m m c = + = π k ∆ α cos 2 k τ − 2 k sin 2 k τ (18) l nc 2 d τ c d τ c The amplitude and phase shift of the circulatory lift coefficient are computed as √ [ ] [ ] 2 2 e e c c | c − ¯ c | = | ∆ c | = c α C F ( k ) − C G ( k ) 2 k + C G ( k ) + C F ( k ) 2 k (19) l l l l 0 F G G F c c α c c e c C G ( k ) + C F ( k ) 2 k G F − 1 c ] ∆ c = tan (20) l c e c C F ( k ) − C G ( k ) 2 k F G c where ∆ c = c − ¯ c − c with c and ¯ c computed by CFD.

l l l l l l c nc The phase shift is also given by ] ∆ c = 2 k τ (21) l lag c where τ is the lag time between the circulatory and quasi-steady state lift coefficients corresponding to the reduced lag frequency k = 0 which is expressed as ∆ ¯ c = c α sin 2 k τ (22) l l 0 α When the quasi-steady state lift coefficient reaches its maximum amplitude at some non-dimensional time τ = τ after the flow establishes a steady state, the circulatory lift coefficient reaches its maximum amplitude sometime later at τ = τ − τ . Thus, in general τ < 0 as the circulatory lift lags the quasi-steady state lift.

lag lag These two nonlinear equations can be solved for C and C for a given amplitude and phase shift of the circulatory F G lift which can determined from CFD simulations. The solutions are given by ( ) e c | ∆ c | cos ] ∆ c + 2 k sin ] ∆ c l l l c c c c [ ] C = (23) F ( ) e c c ∆ α F ( k ) 1 + 2 k l α c ( ) e c | ∆ c | sin ] ∆ c − 2 k cos ] ∆ c l l l c c c c C = [ ] (24) G ( ) e c c ∆ α G ( k ) 1 + 2 k l α c These functions can be computed for a given airfoil and can be implemented as a lookup table or in analytical forms using data regression.

17 of 25 American Institute of Aeronautics and Astronautics B. CFD Simulations of Oscillating Airfoil To conduct the unsteady transonic correction, the airfoil at the MAC wing station is used. OVERFLOW steady-state ◦ ◦ ◦ ◦ simulations are conducted for Mach 0.6, 0 . 7, 0 . 8 at the angles of attack of − 1 . 5 , − 1 , − 0 . 5 , and 0 to obtain the lift curve slope c for the unsteady transonic correction. The OVERFLOW results show an unexpected behavior at l α Mach 0.8 when the lift coefficient decreases precipitously from that at Mach 0.7. To better understand this behavior, FUN3D steady state simulations are conducted over a wider range of Mach number from 0.6 to 1 at an increment of ◦ ◦ ◦ 0.1 which also includes Mach 0.745 and a wider range of the angle of attack from − 2 to 1 at an increment of 0 . 5 .

The FUN3D results also show a similar behavior at Mach 0.8.

Figure 25 shows the unlabeled plots of the lift coefficient of the MAC airfoil computed by FUN3D and OVER- FLOW as a function of the angle of attack for various Mach numbers and the lift curve slope as a function of Mach number at various angles of attack. The lift curve at Mach 0.745 is nonlinear at the angle of attack greater than zero.

This is likely due to stall as the flow begins to separate at the trailing edge. The lift curve at Mach 0.8 is entirely nonlinear over the angle of attack range. The lift coefficient decreases precipitously above Mach 0.745. The lift curves at Mach 0.9 and 1 are generally linear. The OVERFLOW results are in general agreement with the FUN3D results, but there is some discrepancy between the Mach 0.6 and Mach 0.7 results. It should also be noted that the lift coefficients at zero angle of attack are quite large.

TTBW MAC Airfoil Lift Curve TTBW MAC Airfoil Lift Curve Slope o FUN3D = -2 o FUN3D = -1 o FUN3D = 0 o FUN3D = 1 o OVERFLOW = -1 o OVERFLOW = 0 FUN3D M = 0.6 o TSFOIL = -1 FUN3D M = 0.7 l l c c FUN3D M = 0.745 FUN3D M = 0.8 FUN3D M = 0.9 FUN3D M = 1 OVERFLOW M = 0.6 OVERFLOW M = 0.7 OVERFLOW M = 0.8 -2 -1.5 -1 -0.5 0 0.5 1 0.6 0.65 0.7 0.75 0.8 0.85 0.9 0.95 1 , deg M Figure 25. TTBW MAC Airfoil c vs. α and c vs. M l l ∞ α The lift curve slope is computed by cubic spline curvefitting of the FUN3D results. It can be seen that the lift curve slope is strongly dependent on the angle of attack. The lift curve slope has a maximum value at Mach 0.745 and the ◦ ◦ angle of attack of − 2 and decreases rapidly to a minimum value at Mach 0.745 and the angle of attack of 1 . For comparison, the OVERFLOW and TSFOIL results for some angles of attack are also presented on the plot of the lift curve slope.

It should be noted that the CFD simulations do not take into account the wing sweep which acts to reduce the effective Mach number by the factor of cos Λ where Λ is the sweep angle of the elastic axis at the MAC wing station.

So the Mach number in the simulations should be interpreted as the effective Mach number that the airfoil sees. The corresponding free-stream Mach number is higher than the indicated Mach number by the factor of 1 / cos Λ .

Figure 26 shows the velocity contour of the TTBW MAC airfoil which is illustrated by a notional airfoil at various ◦ angles of attack for Mach 0.8. At the angle of attack of − 2 , there is an accelerating supersonic flow region on the lower surface toward the leading edge while the flow on the upper surface reaches supersonic further downstream from the leading edge. There is a retarded flow region on the lower surface where the flow could be separated due to shock-boundary layer interactions. The supersonic flow on the lower surface is the likely cause of the decrease in lift at Mach 0.8. As the angle of attack increases, the accelerating flow region on the lower surface diminishes and the flow becomes entirely subsonic at zero angle of attack. As the angle of attack increases above zero, the shock becomes steeper and the shock-boundary layer interaction causes an increase in the retarded flow region aft of the shock due to the shock-boundary layer interactions. This could be a plausible explanation for the nonlinear behavior of the lift curve at Mach 0.8.

18 of 25 American Institute of Aeronautics and Astronautics Figure 26. Notional Airfoil Velocity Contour at Mach 0.8 Computed by FUN3D Unsteady simulations of the MAC airfoil oscillating about the quarter-chord location are conducted in FUN3D for Mach 0.6, 0.7, and 0.8 at four different reduced frequencies k = 0 . 02, 0.1, 0.2, and 0.3. The FUN3D mesh has about 30 thousand grid points. The mesh domain size is about 100 times the airfoil chord length. The Roe’s flux-difference splitting scheme and the Spalart-Allmaras turbulence model are used in the simulations. An optimized second-order backward finite-difference scheme is used in the time integration. To minimize nonlinear aerodynamic effects as seen ◦ ◦ in the steady-state results, the harmonic angle of attack is chosen with ¯ α = − 1 and α = 0 . 5 . So, the angle of attack ◦ ◦ oscillates between − 1 . 5 and − 0 . 5 . Within this angle of attack range, the steady-state lift coefficients at Mach 0.6 and 0.7 are linear, but the steady-state lift coefficient at Mach 0.8 is nonlinear. The unsteady transonic correction method is applied to the FUN3D-computed unsteady lift coefficients to compute the unsteady transonic correction functions C and C .

F G Figures 27, 28, and 29 show the unlabeled plots of the unsteady lift coefficient computed by FUN3D and the theoretical unsteady lift coefficient computed from the Theodorsen’s theory with and without the unsteady transonic correction method for the reduced frequency k = 0 . 1 at Mach 0.6, 0.7, and 0.8, respectively. The unsteady transonic ◦ correction method uses the steady-state lift curve slope at ¯ α = − 1 . It can be seen that the unsteady transonic correction method is able to correct for the changes in the amplitude and phase shift of the unsteady lift coefficient for the MAC airfoil. For reference, the theoretical unsteady lift coefficient computed from the Theodorsen’s theory is shown in black. The theoretical unsteady lift coefficient with the unsteady transonic correction, shown in red, demonstrates an excellent agreement with the unsteady lift coefficient computed by FUN3D for all the Mach number and the reduced frequencies.

As can be seen, an oscillating airfoil in transonic flow experiences a reduction in the amplitude and an increase in the phase shift in the unsteady lift coefficient relative to the theoretical unsteady lift coefficient. The unsteady lift coefficient response to the angle of attack generally traces out an elliptical trajectory. The elliptical trajectory of the unsteady lift coefficient computed by FUN3D exhibits a smaller inclination angle than that of the theoretical unsteady lift coefficient. This is due to the reduction in the amplitude of the unsteady lift response in transonic flow. The increase in the width of the elliptical trajectory represents the increase in the phase shift of the unsteady lift response in transonic flow. The elliptical trajectory of the theoretical unsteady lift coefficient with the unsteady transonic correction matches precisely that of the unsteady lift coefficient computed by FUN3D for Mach 0.6 and 0.7. A very small discrepancy 19 of 25 American Institute of Aeronautics and Astronautics between the two elliptical trajectories is noted for Mach 0.8.

Figure 30 shows the plots of the unsteady transonic amplitude correction function C and the unsteady transonic F phase correction function C . The unsteady transonic flow simulations exhibit a trend of decreasing the amplitude G and increasing the phase shift of the unsteady lift coefficient response. Thus, the amplitude correction function C is F expected to be smaller than unity while the phase correction function C is expected to be larger than unity. Figure G 30, indeed, confirms this trend. The amplitude correction function C indicates a limiting value of 1 as the reduced F frequency k tends to zero corresponding to steady-state solutions. The amplitude correction function for Mach 0.6 and 0.7 suggests an asymptotic behavior as the reduced frequency k increases, but this asymptotic behavior is missing in the Mach 0.8 case. The trend from the Mach 0.6 and 0.7 amplitude correction functions shows a monotonic decrease in the amplitude with increasing Mach number. The Mach 0.8 amplitude correction function, however, does not follow this trend. This could be due to the nonlinear aerodynamic behavior at Mach 0.8. Thus, increasing the reduced frequency k or Mach number causes the amplitude of the unsteady lift response to decrease. This behavior could be explained by a hypothesis that the accelerating flow as Mach number increases or the frequency of oscillation increases causes an increase in the response time for the circulation around the airfoil to re-adjust itself in order to re-establish the flow field in response to the pitch angle or the angle of attack input. This increase in the response time or lag causes the amplitude of the unsteady lift coefficient to decrease. The phase correction function C is generally greater than G unity and decreases monotonically as the reduced frequency k increases for Mach 0.6 and 0.7. The Mach 0.8 phase correction function does not follow this trend.

TTBW MAC Airfoil Unsteady Lift Transonic Correction (M = 0.6, k = 0.1) TTBW MAC Airfoil Unsteady Lift Transonic Correction (M = 0.6, k = 0.1) FUN3D FUN3D Theodorsen Transonic Correction Theodorsen Transonic Correction l l c c 0 20 40 60 80 100 -1.5 -1 -0.5 = V t/c Figure 27. TTBW MAC Airfoil FUN3D-Computed Unsteady Lift Coefficient with Unsteady Transonic Correction for M = 0 . 6 and k = 0 . 1 ∞ TTBW MAC Airfoil Unsteady Lift Transonic Correction (M = 0.7, k = 0.1) TTBW MAC Airfoil Unsteady Lift Transonic Correction (M = 0.7, k = 0.1) FUN3D Theodorsen Transonic Correction FUN3D Theodorsen Transonic Correction l l c c 0 20 40 60 80 100 120 140 160 -1.5 -1 -0.5 = V t/c , deg Figure 28. TTBW MAC Airfoil FUN3D-Computed Unsteady Lift Coefficient with Unsteady Transonic Correction for M = 0 . 7 and k = 0 . 1 ∞ 20 of 25 American Institute of Aeronautics and Astronautics TTBW MAC Airfoil Unsteady Lift Transonic Correction (M = 0.8, k = 0.1) TTBW MAC Airfoil Unsteady Lift Transonic Correction (M = 0.8, k = 0.1) FUN3D Theodorsen Transonic Correction FUN3D Theodorsen Transonic Correction l l c c 0 20 40 60 80 100 120 140 160 -1.5 -1 -0.5 = V t/c , deg Figure 29. TTBW MAC Airfoil FUN3D-Computed Unsteady Lift Coefficient with Unsteady Transonic Correction for M = 0 . 8 and k = 0 . 1 ∞ TTBW MAC Airfoil Unsteady Transonic Amplitude Correction Function C TTBW MAC Airfoil Unsteady Transonic Phase Correction Function C F G 1 2.8 FUN3D M = 0.6 FUN3D M = 0.6 0.95 FUN3D M = 0.7 2.6 FUN3D M = 0.7 FUN3D M = 0.8 FUN3D M = 0.8 0.9 2.4 0.85 2.2 0.8 F G C C 0.75 1.8 0.7 1.6 0.65 1.4 0.6 0.55 1.2 0 0.05 0.1 0.15 0.2 0.25 0.3 0 0.05 0.1 0.15 0.2 0.25 0.3 k k Figure 30. TTBW MAC Airfoil Unsteady Transonic Amplitude and Phase Correction Functions C and C as Function of Reduced F G Frequencies Computed from FUN3D Results C. Preliminary Transonic Flutter Analysis The purpose of the transonic flutter analysis for this paper is to demonstrate the use of the proposed unsteady transonic correction method in a flutter analysis. To perform a more detail flutter analysis using this method would require conducting more unsteady CFD simulations for other Mach numbers below 0.6 and above 0.8 up to Mach 1 and other reduced frequencies beyond 0.3 as well as for different airfoil sections along the wing. Instead, the currently available unsteady transonic amplitude and phase correction functions C and C computed from the FUN3D simulations are F G extended to lower Mach numbers by requiring that they tend to unity in the limit as the Mach number tends to zero corresponding to the incompressible flow for the Theodorsen’s theory. The steady-state lift curve slope also is extended by extrapolation to a lower Mach number where it matches the classical Prandtl-Glauert compressibility correction formula for subsonic flow in both the amplitude and slope. To deal with the effects of airfoil type and thickness, this would require additional unsteady CFD simulations. It may be reasonable to assume that the airfoil type used in the TTBW is the same for all the wing sections. So, the remaining effect to be accounted for is the airfoil thickness. The airfoil thickness for the five wing stations in the wing-strut interference aerodynamics investigation ranges from 11.6% near the root to 12.5% near the strut juncture. The MAC airfoil thickness almost falls within this range. Therefore, it may be reasonable to assume that the MAC airfoil is sufficiently representative of all the wing airfoil sections. This assumption will be verified later in the future work.

Figure 31 shows the unsteady transonic amplitude and phase correction functions C and C as a function of F G the Mach number and the reduced frequency k . The limited unsteady lift coefficient data computed by FUN3D still 21 of 25 American Institute of Aeronautics and Astronautics enable a reasonable representation of the unsteady transonic amplitude and phase corrections. As can be seen, these corrections vanish as the Mach number tends to zero.

TTBW MAC Airfoil Unsteady Transonic Amplitude Correction Function TTBW MAC Airfoil Unsteady Transonic Phase Correction Function 1 3 0.9 2.5 0.8 F G C C 0.7 1.5 0.6 0 0 0.5 1 0.1 0.1 0.8 0.8 0.6 0.6 0.2 0.2 0.4 0.4 0.2 0.2 k 0.3 k 0.3 0 0 M M Figure 31. TTBW MAC Airfoil Unsteady Transonic Amplitude and Phase Correction Functions C and C as Function of Mach Number F G and Reduced Frequency Computed from FUN3D Unsteady Simulations A preliminary flutter analysis is conducted using the unsteady transonic correction method. At sea level, the first flutter critical mode is associated with the second symmetric mode. No flutter is found for the anti-symmetric modes up to Mach 0.8. For the purpose of illustration, a flutter analysis is conducted at an altitude of 20,000 ft. The structural damping is assumed to be 2%. The tension-induced stiffening effect is included in the flutter analysis. Figure 32 shows the plots of the frequency and damping as a function of the Mach number at an altitude of 20,000 ft with the unsteady transonic correction method applied. The flutter speed corresponding to the second symmetric mode occurs at Mach 0.723.

TTBW Symmetric Mode Frequency w/ Unsteady Transonic Correction @ 20,000 ft TTBW Symmetric Mode Damping w/ Unsteady Transonic Correction @ 20,000 ft 10 0.1 -0.1 -0.2 -0.3 Damping G Frequency, Hz -0.4 -0.5 -0.6 1 -0.7 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 M M Figure 32. TTBW Frequency and Damping of Symmetric Modes as Function of Mach Number at Altitude of 20,000 ft with Unsteady Transonic Correction Method For comparison, a flutter analysis is also conducted using the classical method with the Theodorsen’s theory for incompressible flow and the Prandtl-Glauert compressibility correction formula for the lift curve slope of a swept wing which is given by 2 π cos Λ √ c = (25) l α 2 cos Λ 2 2 1 − M cos Λ + ∞ A . R . ε where A . R . is the wing aspect ratio and ε is the Oswald’s span efficiency factor.

Figure 33 shows the plots of the frequency and damping as a function of the Mach number at an altitude of 20,000 ft with the classical method. The flutter speed corresponding to the second symmetric mode occurs at Mach 0.685.

22 of 25 American Institute of Aeronautics and Astronautics Thus, the unsteady transonic correction method predicts a higher flutter speed than the classical method. This is due to the increase in the amplitude reduction and the increase in the phase shift due to the unsteady transonic correction method which accounts for the increase in the response time of an airfoil in transonic flow.

TTBW Symmetric Mode Frequency w/ Classical Method @ 20,000 ft TTBW Symmetric Mode Damping w/ Classical Method @ 20,000 ft 9 0.1 8 0 7 -0.1 6 -0.2 5 -0.3 Frequency, Hz 4 -0.4 Structural Damping G 3 -0.5 2 -0.6 1 -0.7 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 M M Figure 33. TTBW Frequency and Damping of Symmetric Modes as Function of Mach Number at Altitude of 20,000 ft with Classical Method Figure 34 is the plot of the flutter boundary defined by the flutter Mach number and flutter dynamic pressure computed by the flutter analysis with the unsteady transonic correction method and the classical method. This flutter boundary is associated with the second symmetric mode. The unsteady transonic correction method produces a larger flutter boundary as indicated by the dynamic pressure at the flutter speed. The flutter boundary computed by the unsteady transonic correction is steeper than that computed by the classical method. This flutter boundary appears to eventually dip below that computed by the classical method as consistent with a ’transonic dip’ behavior in transonic flutter. The ’transonic dip’ trend ceases to exist at about Mach 0.750.

TTBW Symmetric Mode Flutter Boundary , psf flutter q Flutter Boundary w/ Unstready Transonic Correction Flutter Boundary w/ Classical Theory 0.5 0.55 0.6 0.65 0.7 0.75 0.8 M flutter Figure 34. TTBW Flutter Boundaries Computed by Unsteady Transonic Correction and Classical Methods Due to the nature of the TTBW MAC airfoil which experiences a sharp reduction in the lift curve slope after reaching the maximum value at about Mach 0.745, the flutter speed computed by the unsteady transonic correction method also follows this pattern. The reduction in the lift curve slope causes a kink in the structural damping of the second symmetric mode as seen in Figure 35 which shows the frequency and damping plots of the symmetric modes 23 of 25 American Institute of Aeronautics and Astronautics

as a function of the Mach number at an altitude of 25,000 ft. This sharp change causes a trend reversal in the structural damping at about Mach 0.750, as can be seen in the structural damping plot. As a result, the flutter tendency begins to subside as the Mach number increases above Mach 0.750. The frequency plot reveals an interesting behavior. As the flutter tendency increases toward Mach 0.750, the second and third symmetric modes begin to coalesce as typically would be expected. As the lift curve slope drops precipitously from Mach 0.745 to Mach 0.8 as seen in Fig. 25, the two modes begin to move away from each other and the structural damping begins to decrease, thereby reducing the flutter tendency.

TTBW Symmetric Mode Frequency w/ Unsteady Transonic Correction @ 25,000 ft TTBW Symmetric Mode Damping w/ Unsteady Transonic Correction @ 25,000 ft 9 0 -0.1 -0.2 5 -0.3 Frequency, Hz -0.4 Stuctural Damping G -0.5 1 -0.6 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 M M Figure 35. TTBW Frequency and Damping of Symmetric Modes as Function of Mach Number at Altitude of 25,000 ft with Unsteady Transonic Correction Method

V. Conclusions

This paper presents a flutter analysis study for the Mach 0.745 version of the Boeing SUGAR transonic Truss- Braced Wing aircraft. The flutter study uses a suite of aerodynamic and finite-element tools to address transonic aerodynamics in steady flow and unsteady flow. An aero-structural model of the Truss-Braced Wing aircraft is de- veloped using VSPAERO vortex-lattice code coupled to an in-house BEAM3D finite-element code. A transonic and viscous flow correction capability is provided by a transonic small disturbance code and an integral boundary layer code. The steady-state aero-structural model demonstrates an excellent agreement in the lift and drag coefficients with the wind tunnel data.

An unsteady transonic correction method is proposed for transonic flow. This method provides two unsteady transonic correction functions to account for the changes in the amplitude and phase shift of an oscillating airfoil in transonic flow. The correction functions are then applied to the Theodorsen’s function to predict the unsteady lift on an airfoil. FUN3D unsteady simulations of a harmonic pitching of the mean aerodynamic chord airfoil are conducted to provide the unsteady lift coefficient for validating the proposed method. The theoretical unsteady lift coefficient computed by the unsteady transonic correction demonstrates an excellent agreement with the FUN3D- computed unsteady lift coefficient, thus validating the proposed method.

Steady state CFD simulations of the mean aerodynamic chord airfoil are conducted in FUN3D and OVERFLOW to compute the lift curve slope for use in the unsteady transonic correction method. The simulation results show a nonlinear aerodynamic behavior as the Mach number increases above Mach 0.745. This nonlinear aerodynamic behavior appears to be due to a high speed flow region on the lower surface of the mean aerodynamic chord airfoil.

The unsteady transonic correction method is applied to a preliminary flutter analysis. The flutter critical mode is associated with the second symmetric mode. The flutter boundary computed by the flutter analysis with the unsteady transonic correction method is greater than that computed by the classical method using the Theodorsen’s theory with the Prandtl-Glauert compressibility correction formula. The nonlinear aerodynamic behavior in the lift curve slope above Mach 0.745 causes a reversal in the flutter tendency of the second symmetric mode above an altitude of 25,000 ft.

In future work, the proposed unsteady transonic correction will be extended to the time domain method to enable the unsteady transonic correction method to be used for modeling flight dynamics of a flexible wing transport aircraft.

24 of 25 American Institute of Aeronautics and Astronautics

Acknowledgment

The authors wish to acknowledge NASA Advanced Air Transport Technologies project for the funding support of this work. The authors also would like to thank Dr. Cetin Kiris of NASA Ames Computational Aerosciences Branch for providing personnel support of this work.

References

Bhatia, M., et. al.. “Structural and Aeroelastic Characteristics of Truss-Braced Wing: A Parametric Study,” Journal of Aircraft, Vol 49, No.

1, 2012.

Gundlach, J. F., Tetrault, P. A., Gern, F. H., Nagshineh-Pour, A. H., Ko, A., Schetz, J. A., et. al., “Conceptual Design Studies of a Strut-Braced Wing Transonic Transport,” Journal of Aircraft, Vol. 37, No. 6, 2000.

Gur, O., Bhatia, M., Schetz, J.A., Mason, W. H., Kapania, R. K., and Mavris, D. N., “Design Optimization of a Truss-Braced Wing Transonic Transport Aircraft,” Journal of Aircraft, Vol. 47, No. 6, 2010.

Bradley, M. K. and Droney, C. K., “Subsonic Ultra Green Aircraft Research: Phase I Final Report,” NASA Contractor Report NASA/CR- 2011-216847, Boeing Research and Technology, April 2011.

Bradley, M. K., Droney, C. K., and Allen, T. J., “Subsonic Ultra Green Aircraft Research Phase II: N+4 Advanced Concept Development,” NASA Contractor Report NASA/CR-2012-217556, Boeing Research and Technology, May 2012.

Nguyen, N., “Elastically Shaped Future Air Vehicle Concept,” NASA Innovation Fund Award 2010 Report, October 2010, Submitted to NASA Innovative Partnerships Program, http://ntrs.nasa.gov/archive/nasa/casi.ntrs.nasa.gov/20110023698.pdf Nguyen, N. and Urnes, J., “Aeroelastic Modeling of Elastically Shaped Aircraft Concept via Wing Shaping Control for Drag Reduction,” AIAA Atmospheric Flight Mechanics Conference, AIAA-2012-4642, August 2012.

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.

Lebofsky, S., Ting, E., Trinh, K. V., Nguyen, N. T., “Optimization for Load Alleviation of Truss-Braced Wing Aircraft With Variable Camber Continuous Trailing Edge Flap,” 33rd AIAA Applied Aerodynamics Confererence, AIAA-2015-2723, June 2015.

Nguyen, N., Ting, E., and Lebofsky, S., “Aeroelasticity of Axially Loaded Aerodynamic Structures for Truss-Braced Wing Aircraft,” 56th AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics, and Materials Conference, AIAA-2015-1840, January 2015. Miranda, L. R., Elliott, R. D., Baker, W. M., “A Generalized Vortex Lattice Method for Subsonic and Supersonic Flow Applications,” NASA Contractor Report 2865, 1977.

Ting, E., Reynolds, K., Nguyen, N., and Totah, J., “Aerodynamic Analysis of the Truss-Braced Wing Aircraft Using Vortex-Lattice Super- position Approach,” 32nd AIAA Applied Aerodynamics, AIAA 2014-2597, June 2014.

Fujiwara, G., Chaparro, D., and Nguyen, N., “An Integral Boundary Layer Direct Method Applied to 2D Transonic Small-Disturbance Equations,” 34th AIAA Applied Aerodynamics Conference, AIAA-2016-3160, June 2016.

Chaparro, D., Fujiwara, G., Ting, E., and Nguyen, N., “Aerodynamic Modeling of Transonic Aircraft Using Vortex Lattice Coupled with Transonic Small Disturbance for Conceptual Design,” 34th AIAA Applied Aerodynamics Conference, AIAA-2016-3012, June 2016.

Chaparro, D., Fujiwara, G. E., Ting, E., and Nguyen, N. T., “Transonic and Viscous Potential Flow Method Applied to Flexible Wing Transport Aircraft,” 35th AIAA Applied Aerodynamics Conference, AIAA-2017-4221, June 2017.

Stahara, S. S., “Operational Manual for Two-Dimensional Transonic Code TSFOIL,” NASA Contractor Report 3064, December, 1978.

Drela, M., “A User’s Guide to MSES 3.05,” MIT Department of Aeronautics and Astronautics, July 2007.

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.

Ting, E., Chaparro, D., and Nguyen, N., “Development of an Integrated Nonlinear Aeroservoelastic Flight Dynamic Model of the Truss- Braced Wing Aircraft,” 58th AIAA/ASME/ASCE/AHS/SC Structures, Structural Dynamics, and Materials Conference, AIAA-2017-1815, January 2017.

Albano, E. and Rodden, W., “A Doublet-Lattice Method for Calculating Lift Distributions on Oscillating Surfaces in Subsonic Flows,” AIAA Journal, Vol. 7, No. 2, February 1969, pp. 279-285.

Bartels, R. E. and Stanford, B. K., “Aeroelastic Optimization with an Economical Transonic Flutter Constraint Using Navier-Stokes Aero- dynamics," Journal of Aircraft, Vol. 55, No. 4, July 2018.

Malik, W., Schetz, J. A. and Kapania, R. K., “Rapid Transonic Flutter Analysis for Aircraft Conceptual Design Applications," AIAA Journal, Vol. 56, No. 6, June 2018.

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

Nguyen, N. and Ting, E., “Inertial Force Coupling to Nonlinear Aeroelasticity of Flexible Wing Aircraft,” AIAA Dynamics Specialists Conference, AIAA-2016-1094, January 2016.

Bisplinghoff, R., Ashley, H., and Halfman, R., Aeroelasticity, Dover Publications, Inc., Mineola, New York, 1996, pp. 261-262.

25 of 25 American Institute of Aeronautics and Astronautics