Document
A Perspective on Computational Aeroelasticity
1970s to Now
Guru P. Guruswamy Computational Physics Branch NASA Advanced Supercomputing (NAS) Division Exploration Technology Directorate Ames Research Center Symposium on Classical to Computational Aeroelasticity Indian Institute of Science Bengaluru, India Sept 10-13, 2018 http://www.nas.nasa.gov/~guru
Acknowledgements
Research Partners Byun Chansup (Sun) : HiMAP, FSI, Parallel Computing Mark Potsdam (US Army) : HiMAP, grids for UH-60A Pieter Buning (LaRC) : OVERFLOW Peter Goorjian (ARC) : TSP, GO3D Doug Boyd (LaRC) : Deforming grid in OVERFLOW, HART-II GRID Dennis Jespersen (ARC) : MPI on Super-cluster Sabine Goodwin, Pradeep Raj, Vin Sharma (Lockheed) : L1011, F-16 Computations Shigeru Obayashi (Tohoku U) : ENSAERO Upwind Flow Solvers Lloyd Eldred (LaRC) : FEM, NASTRAN Neal Chaderjian (ARC) : Zonal Grid Options in HiMAP Rakesh Kapania, Dale MacMurdy (VPI) : FEM Yehia Rizk (ARC) : Hypersonic Fred Striz (Okalohoma U) : FEM, Frequency Domain Approach Manoj Bhardwaj (Sandia, DoE) : FEM Dave Findlay (US NAVY) : F18 Steve Dobbs, Gerry Miller, Iroshi Ide (Boeing) : B-1, X31 Appa (Northrup) Argyris (Stuttgart) : Flight Dynamics for CFD David Yeh (Boeing-Military) : Active Controls, X31 Eugene Tu (ARC) : ENSAERO, Stability Derivatives, Active Controls Lakshmi Sankar: (SA Turbulence model, Parallel version) Joseph Garcia (ARC) : Nonlinear CSD, Controls, Hypersonic Mike Myers (Lockheed): Flutter Ferat Hatay(Sun/Oracle): Parallel Computing Mehrdad Farhannia (President & CEO, Roxwood Medical) : HiMAP Max Platzer, Kevin Jones (Naval Post Graduate School, F-18 Abrupt Wing Stall) Program Managers (NASA) : HPCC, HSCT, SRW, NAVY, Air Force Mentors : Drs. Yang (Purdue), Olsen (AFWAL), Ballhaus (NASA), Ashley(Stanford)
Background
Structural deformations impact vehicle performance - High speed civil transports - Launch vehicles - Rotorcraft - Flexible thermal protection system Strong fluid/structure interactions occur due to - Flow separations, moving shock waves - Large structural displacements - Aeroelastic instabilities such as flutter Current analysis tools for design compute aeroelasticity within linear aerodynamic limitations (NASTRAN) - Euler/Navier-Stokes (ENS) based methods are needed Many current procedures use ENS as corrections to linear aerodynamic - Not adequate when flow non-linearities occur - Cannot account for non-linear phase angles - Not adequate for transient cases - Not suitable for active-controls Recently more tendency towards using Reduced Order Method - Limited to Euler equations - Not robust for Navier-Stokes equations (Stanford) - Not shown better than uncoupled modal approach (ENSAERO)
Key References
“Aeroelastic Time Response Analysis of Thin Airfoils by Transonic Code LTRAN2,” Computers and Fluids, (1981) “Efficient Algorithm for Unsteady Transonic Aerodynamics of Low-Aspect-Ratio Wings,” J. of Aircraft, (1985), “Efficient Algorithm for Unsteady Transonic Aerodynamics of Low-Aspect-Ratio Wings,” J. of Aircraft, (1985) “An Integrated Approach for Active Coupling of Structures and Fluids,” AIAA J., (1989) “Unsteady Aerodynamic and Aeroelastic Calculations for Wings Using Euler Equations,” AIAA J., (1990), “Vortical Flow Computations on a Flexible Blended Wing-Body Configuration,” AIAA Jl., (1992) .
“Direct Coupling of the Euler Flow Equations with Plate Finite Element Structures," AIAA J., (1995).)
“Navier-Stokes Computations for Oscillating Control Surfaces,” J of Aircraft, (1994), “Convergence Acceleration of a Navier-Stokes Solver for Efficient Static Aeroelastic Computations,” AIAA J., (!995) ”CFD/CSD Interaction Methodology for Aircraft Wings ", AIAA J (1998) “Aerodynamic Influence Coefficient Computations Using ENS Equations on Parallel Computers,” AIAA J., (1999) “Development and Applications of a Large Scale Fluids/ Structures Simulation Process on Clusters,” C & Fluids (2007) “Time-Accurate Aeroelastic Computations of a Full Helicopter Model RANS” J. of Aerospace Innovations, (2013).
“Dual Level Parallel Computations for Database to Design Aerospace Vehicles,” NASA TM-2013-216602, (2013).
“Dynamic Stability Analysis of Hypersonic Transport during Reentry,” AIAA J (2017) “Time Accurate Coupling of 3-DOF Parachute System using Navier-Stokes Equations,” J Spacecraft and Rockets, (2016) “N–S Equations-Based Aeroelasticity of Supersonic Transport Including Short-Period Oscillations, “” AIAA J, (2018)
Objective
To present a summary of frequency domain and time domain procedures for aeroelasticity by using non-linear flow equations - Transonic small perturbation theory - Euler/Navier-Stokes Equations - Suitable for frameworks - Focus on efforts at Ames Research Center since 1975 - 100 person-year effort!
Levels of Fidelity
Fluids Structures Detailed 3D finite Navier- Stokes Elements Simple 3D Euler Finite Elements Full 2D FEM & Potential 2D modal Transonic 1D finite Small Interfacing Elements Disturbance Linear Classical Analytical Complexity in geometry Beams Methods Complexity in physics Shape Look-up Functions Tables Aircraft Rotorcraft Hypersonic Vehicles
APPROACH
Solver Approaches
Transonic Small Perturbation Equations - Limited to transonic flows - Robust since no grid movements - Super fast and still good for conceptual design Euler/Navier-Stokes equations - Reynolds averaged Navier-Stokes (RANS) equations - Baldwin-Lomax & Spalart-Allmaras turbulence models - Diagonal form of Beam-Warming central difference solver - Stream-wise upwind algorithm of Obayashi and Goorjian - Structured grids, patched & overset - Implemented in NASA codes HIMAP, OVERFLOW Grids are validated for space and temporal accuracies Lagrange’s structural equations - Finite element and modal form Trajectory Equations 3 DOF - Parachutes Phugoid Motion – Supersonic Transport
Coupled Aeroelastic Equations of Motion
Lagrangian equations of motion • • •
} { } ]{ [ } ]{ [ } ]{ [ f d k d g d m = + +
Where [m], [g], and [k] Are mass, damping, and stiffness matrices {d} and {f} are the nodal displacements and aerodynamic force [ M], [g] and [k] are computed using FEM {F} is computed solving ENS equations Solved using direct time integration
Modal equations of motion
From Rayleigh-Ritz analysis, the displacement vector {d} can be expressed as:
} ]{ [ } { q d y =
] [ y ] [ y where is the modal matrix and {q} is the generalized displacement vector.
The final modal form of Lagrange s equations of motion is: • • •
} { } ]{ [ } ]{ [ } ]{ [ F q K q G q M = + +
where [m], [g], and [k] are modal mass, damping, and stiffness matrices respectively.
{F} is the generalized aerodynamics force vector defined as 2 T
} ]{ [ ] [ C A U y r
P where C p is the pressure coefficient [A] is the diagonal area matrix of the aerodynamic control points, r the free-stream density, and U is velocity.
Coupled Procedures Using CFD/CSD
Can be grouped into two categories based on type of fluid/structure coupling HB-Hybrid TA-time Methods Accurate − F(t) C C { f } C C F S F S Q(t) Q D D D D Every step Ad-hoc € In both methods CFD is computed time accurately In TA, airloads {f(t)} is directly from CFD in same time frame as CSD •• •• [ m ]{ q } + [ c ]{ q } + [ k ]{ q } = { f } CFD − { f } In HB, is modified using non-CFD (look –up tables!)
•• • − [ m ]{ q } + [ C ]{ q } + [ k ]{ q } = { f } = g ({ q }) + { Δ f } CFD - G is based on linear theory, look-up tables etc.
€ - CFD and CSD are not in the same time frame
Uncoupled Aeroelastic Procedure
• Onset of instability starts as a small perturbation phenomenon • Unsteady aerodynamic forces for all frequencies are computed using fast indicial response approach • Primary stability depends on rigid body pitch and plunge motions • Fast generation of indicial responses is accomplished using dual- level parallel computations • Stability analysis is performed in frequency domain using pre- computed unsteady aerodynamic data. (uncoupled analysis)
Two Degrees-of-Freedom Model
Rigid Body Plunge-Pitch modes
• Distances are measured from mid-length +ve towards tail
Stability Analysis Equations
• Frequency domain Eigenvalue Equations ( + ( + δ δ " % μ k [ M ] − [ A ] ) , = λ [ K ] ) , b $ ' # & α α * - * - k b = reduced frequency, ωb/U, ω is circular frequency in rad/sec, U is speed in feet per sec, and b is semi-length. δ = h/2b and α are displacements corresponding to plunging and pitching motions, The mass-to-air density ratio μ = m/(πρb ) where ρ is the air density and m is the total mass.
2 $ !
! $
x " % ω
0.5 C C 0
α h l δ l α & K =
[ ] #
M = # &
[ ] A =
$ '
[ ]
2 2 ω 0 & " x γ r ω γ " % %
π
− C − 2 C
α α α α
# &
m δ m α where w h , w α and w r are plunging, pitching and reference oscillatory frequencies. g a radius of gyration. The eigenvalue λ is defined as μ(1+ ig)(bω r /U) where g is the artificial structural damping.
• Solved using classical U-g method
Indicial Response Method
• History Ø Introduced for computing unsteady airloads - Lomax (1950s) Ø Extended to flight stability analysis – Tobak (early 60) Ø TSP based unsteady aerodynamics - Ballhaus &Goorjian (mid 70s) Ø Airfoil flutter boundary - Guruswamy and Yang (late 70s) Ø NS based wing flutter boundary – Guruswamy and Tu (mid 80s) Ø Stability analysis of spacecraft – Guruswamy (2016) • Assumptions Ø Unsteadiness is small-linear perturbation about a non-linear steady-state solution.
• Advantages Ø Data for multiple frequencies is extracted from a single response • Derivatives Ø Pulse transfer technique - Edwards (mid 80s) Ø Rotorcraft – Leishman (late 80s) Ø Non-linear perturbation - Cummings et. al. (current)
Classical Indicial Method
• Assuming sinusoidal pitching motion
+,#
( )
! # = ! + ! *
0 1 !" ! = − ! ! ( ! ) cos ( ! ! ) ! ! !
!"
!
• Similar equations apply for the moment coefficient
.
• Values for plunge motion are computed using the
h
relation that induced angle
= a
i
U
Typical Finite Elements Used
In - house NASA 1D, 2D and 3D CSD models - BEMBLD : 10dof rotating beam element - SPARH: shear panel for spars and ribs - PLTSHL: 18 - dof skin/plate/shell fem - TET3D: 12 - dof tetrahedron solid element (future) - NASTRAN SPARH BEMBLD PLTSHL Bembld * Shake test TET3D First flapping frequency
CFD/CSD Interactions (FSI)
Consistent load approach Illustration for transverse DOF
P(x)
F
F
F = ò p(x) s (x) dx 0 < x < l
i i where i denotes the deg of freedom p(x) is the distributed load and S i (x) IS ith shape function Work is conserved between CFD and CSD Extensions using virtual surface method Guruswamy, G. P.,“ Coupled Finite-Difference/Finite-Element Approach for Wing-Body - Aeroelasticity," AIAA 92-4680, n, September 21-23, 1992, Cleveland, Ohio.
Guruswamy, G.P. ," A Review of Numerical Fluids/Structures Interface Methods for Computations using High Fidelity Equations, "Computers and Fluids, (80) 200 2 (Editor Prof Bathe, MIT)
Fluid/Structure Interface
Lumped load approach - Fast, needs fine grids, adequate for uncoupled method Consistent load approach (conserves loads) - Accurate for coupled methods, expensive -1 T T -1 Y Y Y D Y [ ] • Transformation matrix : + K T = A S S S T Q Q • Deflections at cfd grid : = A S T Z Z • Nodal forces at fe grid : T = A S T T T T • • • • Q Z Q Q Z Q Z Z = = S A S A Cfd Fe grid Intermediate Grid Grid Y D, Q Q Q Q K, F( , ) Y = A = S A S Guruswamy, G.P, and Byun, C. "Direct Coupling of the Euler Flow Equations with Plate Finite Element Structures,'' AIAA J, Vol. 33, No 2, Feb 1995 .
Newmark’s Time Integration
Assuming linear acceleration: • •
}) ]{ [ } ]{ [ } ]({ [ } { w K v G F D q - - =
t t 1 -
ö æ
ö æ D
t t
ö
æ D
÷ ç
÷ ç
] [ ] [ ] [ ] [ + + = K G M D
÷ ç
÷ ç
÷ ç
6 2
ø è
ø è
ø è
• • •
t
ö æ D
q q v + = } { } {
÷
ç
t t t t D - D -
ø è
• • • • •
t
ö æ D
q q t q t q w + D + D + = } { } ){ ( } ){ ( } {
÷ ç
t t t t t t t t D - D - D - D -
ø è
• • • • • •
t t
ö æ D ö æ D
q q q q } { } { } { } { + + =
÷ ç ÷ ç
t t t t t t D - D -
2 2
ø è ø è
2 2 • • • • •
ö æ D ö æ D
t t
÷ ç ÷ ç
q q q t q q } { } { } ){ ( } { } { + + D + =
t t t t t t t t D - D - D -
÷ ç ÷ ç
6 3
ø è ø è
Above can be made implicit, assuming constant acceleration. This is NOT loose coupling as referred by Chopra’s, group at UMD, ,Barun’s group at Georgia Tech etc
FSI Implicit (staggered) algorithm
Algorithm a) Compute deformations using loads at t b) Using a) advance CFD from time t to t + ∆t c) Compute deformations using loads at t + ∆t d) Repeat a) b) c) until converge (use Newton method) Increases book-keeping and dds significant computational expense No impact on CFD and CFD (linear and geometrically non-linear) computations - In-house research - NASA sponsored research grants (UMD, Stanford)
Some Related Developments
(in-house and/or sponsored) Sheared grids for swept tapered wings (TSP) Characteristics BCs for supersonic flows (TSP & ENS) Stream-wise upwind algorithm for ENS solvers Pipe-line Gauss-Seidel algorithm for parallel CFD Parallel direct & sub-structure solver for CSD (NASTRAN) MPIRUN, for parallel fluid/structure/control simulation Parallel version of SA turbulence model Virtual surface method for FSI Virtual surface method for sliding CFD grid zones (controls) Area co-ordinate method for FSI(wing-box) Tcl/Tk, C++ for CFD/CSD communications (GUI) Staggered FSI (rotorcraft)
Dual-Level Parallel Computing
• Efficient single-job environment for multiple cases, each running on multiple cores: - Reduces system start/end overhead - Makes sure that cores do not overlap - All cases are completed at the same time, enabling designer to plan his work accordingly rather than “baby-sitting” multiple jobs • Facilitates fast generation of indicial response data for use in stability analysis Guruswamy, G.P., “Dual Level Parallel Computations for Large Scale High- Fidelity Database to Design Aerospace V ehicles,” NASA/TM-2013-216602 , Sept. 2013.
Load Balancing Scheme in HiMAP
Start Select A Block Yes Is > Optimum Memory Size ?
Go To No Partition Next Block Are Blocks No Completed ?
Yes Assign A Block To A Node Next Block Is Node Full ?
No Next Block Yes Next Node No Are Blocks Completed ?
Yes Stop LOAD BALANCING Node filling algorithm for cfd zonal grids Mesh partitioning for csm Results Of Load Balancing Typical 34 Zone CFD10m Mesh .
Number nodes were reduced To 28 from 34
Three Level Parallel Computing
Description of Higfidelity Multidisciplinary Analysis Process (HiMAP) A 3-level parallel, meta modular, multidisciplinary analysis software that runs on single-image shared / distributed memory supercomputers using MPIAPI middleware Accomplishments Can handle disciplines based on time accurately coupled high fidelity methods - Structured/unstructured grid-based Euler/Navier-Stokes solvers (ARC3D, GO3D,USM3D) with parallel multi-block moving grid - Modal and parallel finite element structures (NASTRAN) - Time domain controls - Winner of NASA Sofware Release award Applications Demonstrated for full F-18, L1011, 777, HSCT, X-31 and UCAV configurations MPIAPI(C++) CONTROLS FLUIDS STRUCTURES
Multizonal/Multidiscipline Parallel Communication
in HIMAP
3-level Parallel Communication In HiMAP
Effort towards Framework
RUNEXE - C++ Based Super Modular Process All modules are treated separately as objects and executed using C++ system commands Communication among modules is accomplished using I/O, MPI and/or TCL/TK On-the-fly based graphics, xmgrace, opengl Manages data for visualization (field view )
C++
PLOT START CFD FTOS CSD ?
STOF GRID CURRENT FORTRAN 95 HAS SOME SIMILAR CAPABILITIES
Illustration of RUNEXE Process
ARC2D TRI2D
C++
XMGRACE OpenGL
2D Transonic Small Perturbation Theory
ATRAN2S (1977)
Guruswamy, G. P. and Yang, T. Y.," Aeroelastic Time Response Analysis of Thin Airfoils by Transonic Code LTRAN2," Computers and Fluids, Vol. 9, No. 4, Dec. 1981,
3D Transonic Small Perturbation Theory
ATRAN3S (1984)
Guruswamy, G. P. and Goorjian, P. M. “ Computations and Aeroelastic Applications of Unsteady Transonic Aerodynamics About Wings," Jl. of Aircraft, Vol. 21, No. 1, Jan. 1984, pp. 37-43.
Fluid/Structures Interaction
HiMAP with NASTRAN NASTRAN/ANS4 Cp at M = 0.85 Elements Demonstrated for HSCT model with 5 million grid points fluid and 20k DOF ELFINI FEM Eldred, L. , Guruswamy, G.P. and Byun, C, "Parallel Aeroelastic Analysis Using ENSAERO and NASTRAN,” HPCCP/CAS Workshop 98, NASA/CP-1999-208757, Jan 1999.
Parallel FSI Demonstration for HSCT
Nastran Based Parallel Sub-structure Solver Is Developed Nastran/Ans4 Elements Tca6 20k Dof Fem Demonstrated for a HSR Model With 5 M Pt Fluid and 20K DOF FEM
Validation Of Unsteady Pressures
NASA TND 344 WING, M = 0.90, k = 0.26 UNSTEADY Cp AT 50% SEMISPAN WING IN FIRST BENDING MODE MOTION TOTAL LIFT MOTION Guruswamy, G. P. ,"Unsteady Aerodynamic and Aeroelastic Calculations for Wings Using Euler Equations ," AIAA Jl., Vol. 28, No. 3, March 1990, pp 461 - 469. (also AIAA Paper 88 - 2281)
Rectangular Wing
Aeroelasic validation NASA TMX-79, AR =5, M = 0.715, re = 4.5 million Modes Dynamic Aeroelastic Responses The picture can't be displayed.
MEASURED DYNAMIC PRESSURE q =1.31psi PLATE FEM AND MODAL GIVE SIMILAR RESULTS Guruswamy, G.P.,"Computational-Fluid-Dynamics and Computational-Structural-Dynamics Based Time-Accurate Aeroelasticity of Helicopter Blades," Jl. of Aircraft, Vol. 47, No. 3, May-June 2010
Validation For Aeroelastic Research Wing
HIMAP, M = 0.80, a = 2.98 deg, Re = 1.3M, q = 0.72 psi
Surfa c e Pres s ure s a t 71% Span SURFACE PRESSURES AT 71% SPAN -1 . 5 -1 -0 . 5 Cp FLEXIBLE Compu tatio n (fle x ible ) Compu tatio n (rig id) 0 . 5 RIGID Ex p erime n t EXPERIMENT 0 0 . 2 0 . 4 0 . 6 0 . 8 1 FEM MODEL OF WING - 400 DOF X/C COMPUTED WIND - OFF (COMPUTED) o EXPERIMENT Bhardwaj , M., Kapania , R., Reichenbach and Guruswamy, G.P., “ CFD/CSD Computational Interaction Methodology for Aircraft Wings, ” 6 MODES FROM GVT GIVES SIMILAR RESULTS AIAA Jl Vol 36, No 12, Dec 1998.
Vortex Induced Aeroelastic Oscillations of
Supersonic Aircraft
Guruswamy, G. P., ,"Vortical Flow Computations on a Flexible Blended Wing-Body Configuration," AIAA Jl., Vol. 30, No. 10, October 1992. pp 2497-2503
Demonstration for Full Aircraft
(HIMAP, L1011 Wind Tunnel Model, M = 0.85) 9m grid pts, 38 fluid zones 5 structural modes, 700 nodes with 3dof per node Cp deformed configuration Cp Typical aeroelastic computation using 34 nodes of parallel computer requires - 15 minutes with memory efficient load balance scheme or 10 CPU hrs without load balance scheme ( Potsdam, M.A. and Guruswamy G.P , " A Parallel Multiblock Mesh Movement Scheme for Complex Aeroelastic Applications, " AIAA 2001-0716, Jan 2001. )
X-31- Pressure Distribution
HiMAP Applications
L1011, 10M PTS MODALSTRUCTURES HSCT -12 M GRID POINTS ELFINI FEM F18E/F, 17M PTS FEM STRUCTURES UCAV 14M GRID POINTS
Rotating Blades
Unsteady Validation- Caradonna-tung Blade, AR =16, Re = 3.93m, RPM =1500, θ = 0.0deg
Aeroelastic Validation for Rotating Blade
Advancing Blade, Japan/MIT Ω =100 RAD/SEC, Μ = 0.40, Flexible Blade, Θ c = 0 Deg BEMBLD PITCH CHANGE RIGID BEARING SHANK Guruswamy, G.P.,"Computational-Fluid-Dynamics and Computational-Structural-Dynamics Based Time-Accurate Aeroelasticity of Helicopter Blades," Jl. of Aircraft, Vol. 47, No. 3, May-June 2010
HART II Configuration
o Advance Ratio = 0.15, RPM = 1041, Shaft Angle = 4.5 Wind Tunnel Model - 27m points with 32 near body grid blocks - 23 m points outer grids Guruswamy, G. P, "Time-Accurate Aeroelastic Computations of a Full Helicopter Model using the Navier-Stokes Equations," International Jl. of Aerospace Innovations, Vol. 5, No 3+4, Dec 2013, pp. 73-82
Air loads for Hart II Configuration
Baseline, Advance Ratio = 0.15, RPM = 1041 TA responses Fourier Analysis - Isolated blade (Full Configuration ) - 2m grid points, 50DOF FEM th Converged response at 6 revolution - 7200 time steps per revolution - 250 cpu hours - 4 hours wall clock time with 64 processors
Demonstration for Launch Vehicles
EFFECT OF MODES ON TRANSONIC AIRLOADS 1,000 Transonic aeroelastic responses computed with 30 minutes of wall-clock time by using 1,000 nodes each with 4-openmp cores and MPIEXEC utility developed at NAS Guruswamy, G.P., "Large-Scale Computations for Stability Analysis of Launch Vehicles Using Cluster Computers," Jl of Spacecraft And Rockets, Aug 2011.
Demonstration for Launch Vehicles
Unsteady Motions
Snap shot of unsteady pressure contours ( V = 0.005, L, H = 0.005L, k = 0.5, M = 1.8,) when h and v are maximum Effect of M on lateral forces Effect of M on longitudinal forces Each 42-case (13M grid pts 5 cycles) job required a total 25 hrs of wall clock Guruswamy, G.P Navier-Stokes based Unsteady Aerodynamic Computations of Launch Vehicles undergoing Coupled Oscillations., AIAA Jan 2013.
Trajectory Motion of Parachute
Trajectory Equations of Motion
The following assumptions are made: a) m a depends only on the canopy.
b) The centers of m c and m a are coincident.
c) m c and m p are at a fixed distance apart of L.
d) The centers of forces on canopy and its mass are coincident) e) The aerodynamic force on payload smaller compared to canopy The equations of motion governing the system are written as: " ! # & 12
+ # − = 4 (1)
.)
" !$ '()*+, +- ) '(+-*+,) " ! 6
+ − + & + 1 ,85# − 1 59:# = 4 (2)
5 $ 7 2 " !$ " ! ;
+ + 1 59:# + 1 ,85# = 4 (3)
5 7 2 " !$
Parachute Trajectory Motion
M = 2.0
Effect of structural damping on Responses during descent responses.
from M ∞ = 2.0.
Parachute Cluster
• Positions of canopies are initialized by applying transformations
and rotations to undeflected canopy using Config.xml input file of OVERFLOW - Radius of rotation for all canopies 2D - Canopy_1 by 15 degrees in the X-Z plane, - Canopy_2 by -15 degrees in the X-Z plane then -15 degrees in the X-Y plane, - Canopy_3 by 15 degrees in the X-Y plane
Parachute Cluster Grid
• Xray and Hole Cutting tools of OVERFLOW are used to blend near body grids with off-body grid.
Parallel Computations
Steady state computations on isolated canopy 12000 iterations, 26.5 miilion grid points 100 cases using 4000 cores with 40 cores per case requires 4.61 hrs - 1.8% more time than to run a single cases Coupling with trajectory motions is in progress
Example for Uncoupled Computations
Start Computations From Flutter Boundary Initial Conditions Of A Typical Wing Case_1 Case_m Case_2 - - - 0.5 . . . .
. .
Mn M1 M2 Mn M1 M1 M2 Mn M2 0.4 0.3 Flutter Speed 0.2 FLUTTER SPEEDS FOR ALL CASES BY U-g METHOD HiMAP EXPERIMENT 0.1 Mn N Modes = Number Of Modes X Number Of Frequencies Selected 0.6 0.7 0.8 0.9 1 1.1 1.2 Mach Number Case_m Each Case May Have Different Flow Conditions Note : 100 Gflop Performance On 1024 Nodes On O2000 Suitable Low- Fidelity Computations To Fill Design Space
Where CFD/CSD made a Difference
Transonic flutter-dip of transport aircraft Lateral vortex motion coupled with bending motion of blended wing body configuration (B-1, HSCT, BWB) Control-reversal due to moving shock-waves Jump in phase angles near shock-wave Leading edge vortex induced vertical tail oscillations (F18) Nacelle oscillations of aircraft in transonic regime (L1011) Blade vortex interactions of rotorcraft Flexible thermal protection system (in progress)
Phugoid Motion Simulation of a Supersonic
Transport using Navier-Stokes Equations
Oscillation in altitude due to the exchange between potential energy and kinetic energy is called phugoid oscillation. Beginning at the bottom of the cycle, pitch angle (θ) increases as the aircraft gains altitude and losses forward speed (V). During phugoid motion, the angle of attack (α) remains constant so that a drop in forward speed amounts to a decrease in lift and flattening of the pitch attitude.
Phugoid Motion Equations
Assuming that the phugoid motion starts with level flight, the equations of motion are written as: & − ( # # ! # = (1) *+ # !"
$ $ ) , # , where u is the change in the velocity from the initial velocity u 0 , θ is the flight-path angle, and g is acceleration due to gravity.
X u and Z u are defined as: -.
& = − 01 + 41 (2) # 2, 24 /# , -.
+ = − 01 + 41 (3) # 5, 54 /# , Equation system (1) are combined into a single ordinary & #*# # differential equation with u as a variable by using $ = ( which results in: + ( # ̈ # − & # − # = ,. , (4) #̇ # , In this work, Eq. (4) is solved using the Newmark’s time integration method
Typical Supersonic Transport
Effect of Mach Number on Phugoid Responses
Effect of Mach numbers on Effect of Mach number on change in speed at α = 5 deg. oscillation period α = 5deg
Dynamic Stability Analysis
Hypersonic Transport During Reentry
• Strong need exists for development of faster civil-transport - Supersonic transports - Boeing program, Next generation Concorde - Hypersonic Civil Transport (HCT) - NASA programs, European SKYLON • Successful NASA Space Shuttle Transport (SST) - Limited passenger capability, 6 crew members - Lower length to width ratio - Stable atmospheric re-entry trajectory • Current HCT - Planned for larger passenger capability, ~50 passengers - Larger length to width ratio - Aeroelastic stability plays more important role
Dynamic Stability Analysis
Background
• Stability characteristics of high speed vehicles during re-entry - Centre of pressure (x CP ) shifts significantly and affects stability - Trimming is not practical due to rapid changes - Often expensive approaches such as moving the fuel location are needed for stable flights • Computational approaches - Current stability analysis are limited to use linear aerodynamics - Fast methods based on Navier-Stokes flow equations are needed Typical Reentry Scenarios Statically Unstable Statically Stable Pitch-Plunge Dynamically Unstable Motion Statically Stable Dynamically Damped Statically Stable Dynamically Over -damped
Hypersonic Civil Transport (HCT)
• Configuration
- Topology of Langley Glide Back Booster (LGBB) • Grids are generated using OVERGRID following accepted engineering procedures (NASA/TM-2013-216601) • Body-fitted structured overset 19 near-body grid blocks (20 million points) - A normal spacing of 0.000025 of wing chord - Surface stretching factor of 1.125 - Typical y+ value (one grid point away from the surface) is around 1 • Cartesian overset back ground grid blocks (25 million points) - Resolution of Level 1 grid is 0.5 % of length - Outer boundary location 12 lengths
Parameters for Test Cases
(Final 30 seconds of typical atmospheric reentry) • Mach number (M ∞ ) 5.5 to 0.5 • Angles of attack (α) 12.0 to 2.0 degrees • Oscillating frequencies (f) plunge : 2 Hz, Pitch = 8hz • Coefficients validated - Steady pressure C p , Unsteady force coefficients • Stability parameters computed - Flutter speed and frequency
Reentry – Steady- State Computations
Full Configuration, M ∞ = 5.5 to 0.5, α = 12 to 2 deg • Using ~4,000 cores - 100 cases in 2.5 hrs wall clock time
Validation of Indicial Computations with Experiment
Oscillating NACA64010 Airfoil M ∞ = 0.8, Amplitude of α = 1.0 degs, Re_c = 2.0 million 6000 Steps
Validation of Indicial Computations
with Time Integration
Oscillating, 50% Semi span section of HCT Wing M ∞ = 0.9, Amplitude of α = 0.5 degs, Re_c = 2.0 million • Indicial response (IN) converged in 3000 steps • Time-integration (TI) required 3 cycles with 3600 steps per cycle for each frequency • Flutter speed differed by 5% • Computational speed-up of IN (TI/IN) = (3x5x3600)/3000 = 18
Indicial Response Computations for Reentry
Full Vehicle, 100 cases with M ∞ = 5.5 to 0.5, α = 12.0 to 2.0 degs
Moment Lift • Using ~4,000 cores.
- 100 indicial responses in 17.5 hrs wall clock time
Flutter Boundary Computations during Reentry
Elastic axis is at 67% of length from nose
100 cases with M ∞ = 5.5 to 0.5, α = 12.0 to 2.0 degs
x
α M ∞ = 2.5
x
α
Flutter Boundary Computations during Reentry
Effect on Center of Pressure location
100 cases with M ∞ = 5.5 to 0.5, α = 12.0 to 2.0 degs
Mass Center • Based on Classical Theory - Flutter speed decrease with increase in lift force - Flutter speed increase as x CP moves closer to mass center - Phenomenon is similar to dip in flutter speed for wing of supersonic aircraft in the transonic regime
LaRc BACT Model
Unsteady Computations on Oscillating Flap M ∞ = 0.77, δ 0 = 0.0 deg, Re_c = 3.86 million Sheared Grid Approach α = 4.01 deg, k= 0.22, δ β = 3.86 deg
Wing Body with Oscillating Control Surfaces
Active Controls using TSP
Modes
Demonstration of HiMAP for
Wing-Body-Control Aeroelasticity
Processors Fluid - 8 Structures -1 FLAP DOWN
Controls -1
FLAP UP -5 500 PRESSURE TIP DEFLECTION ITERATION
Summary of Results from Load Balancing Scheme
Grid In Millions Blocks In 10s Efficiency Factor(e) 777 SST L1011 UCAV Efficiency factor ‘E’ is defined as Where zo is number zones before load balancing and zn is the number of regrouped zones after load balancing.
Some Validation Efforts Since 1978
CFD- based frequency domain uncoupled methods Classical methods such as Theodorsen’s theory Linear aerodynamics theories - Kernel function method - Doublet-Lattice method (NASTRAN) Wind tunnel - NASA TND344: unsteady aero, rectangular wing - NASA TMX79: flutter of rectangular wing - NASA ARW: aeroelasticity of wing body - B-1 aircraft: vortex induced oscillations - F-18 aircraft: vertical tail buffet - HART II rotorcraft: blade vortex interaction - F5-wing: oscillating control surfaces - AFW: aeroelastic flexible wing, active controls - L1011: Nacelle oscillations - NASA/RAE/NLR – Moving control surfaces Flight tests - B-1 aircraft - F-18A -aircraft - UH-60A rotorcraft
What is Happening Elsewhere
Rigorous Math models are being developed for FSI - Not yet proven better than engineering approaches Reduced order models for flows - Still as good as modal approach for Euler - Not robust for use with Navier-Stokes equations Using unstructured CFD - Computationally very slow - Still has issues in viscous zones Lattice Boltzmann equations based NS - Bookkeeping can be nightmare!
- MDAO frameworks - One too many but still lack robust high fidelity MDA tools - Often end-up with low fidelity method - OpenMDAO supported by NASA has potential
Conclusions
A summary of about 35-years, 100-person-year effort on development and applications of accurate coupled and uncoupled aeroelastic procedures using non-linear aerodynamic is presented.
Case-by-case validations with either experiment or linear theory are shown to establish aeroelastic computations From this research it is observed that - Transonic small perturbation theory still useful for conceptual design - Euler/Navier Stokes equations and time accurate couplings are required to predict phase angles - Modal approach is adequate for responses - Staggered time integration is not necessary - Indicial method is better suited with NS than Reduced order modeling Given present computational resources, there is no need to hybridize CFD with linear/empirical aerodynamics - Introduces significant errors when inertial loads are present, which is common for flexible configurations - Not capable of predicting transient physics associated with coupling of non-linear flows with structures Classical Indicial method is robust and more efficient than reduced order method for RANS
Future Efforts
Coupled procedures need higher fidelity CSD - 3-D FSI - Composites for all aerospace vehicles - Visco-thermoelastic for spacecrafts Faster and more accurate CFD - Better turbulence models - Robust moving grids - Larger time steps with no CPU time penalty - Robust codes with real-gas effects - Multi-Phase such as Cavitational flows for PoGO of Launch Vehicles) CFD/CSD method time accurately integrated with stability and maneuver dynamics equations is needed for full simulation Appa, K., Argyris J.H and Guruswamy G. P., "Aircraft Dynamics and Loads Computations Using CFD Methods," AIAA 96-1342 .
Future Efforts
Hypersonic Inflatable Aerodynamic Decelerator
Simplified 1-D FEM Model Used Elsewhere 3-D Shell FEM Model Preliminary Results