Steady-State Computation of Constant Rotational Rate Dynamic Stability Derivatives

20 Pages·2000·1.1 MB·English

by  Park, Michael A.

AIAA 2000-4321
Steady-State Computation of
Constant Rotational Rate
Dynamic Stability Derivatives
Michael A. Park
George Washington University
Joint Institute for the Advancement of Flight Sciences (JIAFS)
Lawrence L. Green
NASA Langley Research Center, Hampton, Virginia
18th Applied Aerodynamics Conference
14-17 August 2000
Denver, Colorado

AIAA-2000-4321
STEADY-STATE COMPUTATION OF CONSTANT ROTATIONAL RATE
DYNAMIC STABILITY DERIVATIVES
Michael A. Park*
George Washington University Joint Institute for the Advancement of Flight Sciences (JIAFS)
[e-mail]
Lawrence L. Green :
NASA Langley Research Center, Hampton, Virginia
l.l.green @larc.nasa.gov
Abstract Introduction
Dynamic stability derivatives are essential to Dynamic derivatives quantify the aerodynamic
predicting the open and closed loop performance, damping of aircraft motions and are used to predict the
stability, and controllability of aircraft. Computational longitudinal short period, lateral pure roll. and lateral
determination of constant-rate dynamic stability Dutch roll behavior of the configuration. Analytical,
derivatives (derivatives of aircraft forces and moments empirical, and vortex lattice methods of estimating
with respect to constant rotational rates) is currently these derivative values are not suited to unconventional
performed indirectly with finite differencing of configurations or high-speed, compressible flows
multiple time-accurate computational t]uid dynamics dominated by viscous effects. Evaluating
solutions. Typical time-accurate solutions require unconventional configurations is of growing interest
excessive amounts of computational time to complete. due to the design and analysis of next generation
Formulating Navier-Stokes (N-S) equations in a attack, transport, and reusable launch vehicles.
rotating, noninertial reference frame and applying an Examples of these new, unconventional designs are the
automatic differentiation tool to the modified code has blended wing body and the X-33 configurations. A
the potential tbr directly computing these derivatives methodology of using high fidelity, noninertial Euler
with a single, much faster steady-state calculation. The and Navier-Stokes (N-S) calculations gives improved
ability to rapidly determine static and dynamic stability capability in predicting these dynamic stability
derivatives by computational methods can benefit derivative values in compressible flow on conventional
multidisciplinary design methodologies and reduce or unconventional designs.
dependency on wind tunnel measurements. The Due to cost and time limitations, it is impractical to
CFL3D thin-layer N-S computational fluid dynamics construct and test numerous wind tunnel models during
code was modified tot this study to allow calculations initial prototyping. Therefore, measurement of the
on complex three-dimensional configurations with effects of aircraft dynamics on preliminary
constant rotation rate components in all three axes. configuration aerodynamic forces and moments is
These CFL3D modifications also have direct limited. The application of automatic differentiation to
application to rotorcraft and turbomachinery analyses. a noninertial reference frame Euler and N-S code has
The modified CFL3D steady-state calculation is a new potential for providing designers with insight, gained
capability that showed excellent agreement with results from higher fidelity codes, into aircraft dynamics at the
calculated by a similar formulation. The application of preliminary design stage. This design stage is when
automatic differentiation to CFL3D allows the static control surface size and preliminary control laws are
stability and body-axis rate derivatives to be calculated being evaluated. Computational determination of these
quickly and exactly. derivatives is cheaper and faster than performing wind
tunnel measurements and will aid rapid prototyping
*Graduate Student. NASA Langley Research Center, and multidisciplinary design.
Multidisciplinar_()ptimizationBranch.MS 159,Hampton,Virginia, The modification of the CFL3D ] (Computational
MemberAIAA
Fluids Laboratory Three-Dimensional) computational
÷
Research Scientist, Muhidisciplinary Optimization Branch, MS fluid dynamics (CFD) code to perform calculations in
159.SeniorMemberAIAA
a noninertia[, rotating reference frame has the potential
Copyright©2000bytheAmerican InstituteofAeronautics. Inc.No to reduce the reliance on forced-motion wind tunnel
copyrightisassertedintheUnited StatesunderTitle 17,U.S.Code. and free-flight wind tunnel tests. Considerable
TheGovernment hasroyahy-free licensetoexerciseallfightsunder
previous work performed on turbomachinery has
the copyright claimed herein for government purposes. All other
fightsreservedbythecopyrightowner. demonstrated noninertial, rotating reference frame
1
American Institute of Aeronautics and Astronautics
fluid mechanicass a meansto greatlyreduce
study has been performed in an inviscid, Euler mode
computationtiamle(foranexamplseeeRef.2).Kandil
and a viscous mode with the N-S equations coupled to
andChuang
3"4have demonstrated noninertial reference the Spalart-Allmaras (S-A) turbulence model. I
frame calculations tor general motions on rolling
aircraft stability problems. Limache and Cliff 3devised CFL3D Noninertial Reference Frame Modifications
an efficient scheme for the special case of steady-rate There are two reference frames depicted in Fig. 2:
motion and applied this technique to stability and the inertial reference frame (denoted with upper-case
control work with a two-dimensional (2-D), symbols) and the noninertial frame (denoted with
unstructured grid code and the sensitivity equation lower-case symbols). The CFD grid (depicted as a
method. cube) is embedded in the noninertial reference frame.
The noninertial modifications to CFL3D were Positions relative to each of these two reference frames
initially validated in this study for a 2-D NACA0012 are quantified by three scalar quantities (X,)I, Z and x,
airfoil case with comparisons to previously published y, z) that describe location along three orthonormal unit
results by Limache and Cliff. 5The modified CFL3D vectors (1, J, K and i,j, k). Note that bold type face
was then applied to the lull three-dimensional (3-D) indicates vector symbols. The inertial frame is fixed in
Lockheed Martin Tactical Aircraft Systems-- space and the noninertial frame can translate and rotate
Innovative Control Effectors TM (ICE) _ configuration 6 with the rotation described with three scalar
(Fig. 1)with a turbulent Navier-Stokes calculation. components (60_,0_, _.) of the rotation vector o3. The
noninertiai frame and CFD grid tollow a curved path
Technical Approach (denoted as the curved, dashed arrow) as they
This study adopted the Limache and Cliff 5 simultaneously translate and rotate.
approach. There were two major aspects to this project. Each of these coordinate systems has advantages
The first was modifying CFL3D to perform and disadvantages. An advantage of the inertial frame
calculations in a rotating, noninertial reference frame. is that it is not moving (the existing stationary grid N-S
These CFL3D modifications included adding a source equations in CFL3D are only formulated for
term to the residual calculation and modifying the nonmoving frames). An advantage of the noninertial
boundary and initial conditions. The second aspect was reference frame is that it moves with the CFD grid;
the application of ADIFOR 7`s (Automatic therefore, the stationary grid tbrmulation of the
Differentiation in FORTRAN) to the latest parallel CFL3D N-S equations is already coded in this frame of
version of CFL3D. This code was used to compute reference with local (lower-case) variables. A
derivatives of aircraft forces and moments with respect disadvantage of the noninertial reference frame is the
to the flow angles and constant rotational rates in the current, stationary grid N-S equations are not
roll, pitch, and yaw axes. The application of ADIFOR formulated correctly because the CFD grid and its
to the unmodified version of CFL3D has been associated reference frame are rotating (e. g.,
performed successfully to calculate static stability accelerating) in this study to simulate aircraft
dcrtv•attve.s 0 (dertv...a.ttves ol atrcraft forces and moments constant-rate motion.
with respect to angle of attack and angle of sideslip). In order to correctly modify the stationary grid N-S
equations to calculate valid solutions with a translating
CFL3D Introduction and rotating CFD grid, the relation between the
The CFL3D code is a FORTRAN 77 (F77) descriptions of the same point in both reference frames
Reynolds-averaged thin-layer N-S flow solver tbr (b and B) must be sought. Note that all of these
structured-volume grids. CFL3D was written primarily derivations are performed at the instant in time when
at NASA Langley Research Center and is undergoing the unit normal vectors of both systems are parallel,
continuous development and improvement. The code which removes the necessity of a rotational coordinate
has the ability to compute inviscid Euler, laminar N-S, transform. Also, the noninertial reference frame is
and turbulent N-S calculations. The code employs translating and rotating at a constant rate: therefore
parailelization by decomposing the computational angular acceleration, 60, is zero. A nonzero value of 6_
domain into many separate blocks. These individual can be included to model more general motions, 3"4but
blocks are analyzed in separate processes that such was not the intention of this study.
communicate with each other by means of the Message There are two points in space of interest tbr this
Passing Interface (MPI) standard. Analysis for this derivation: The position of the noninertial frame
origin, C, expressed as a function of the inertial
Theuseof trademarks ornamesof manufacturers inthisreportis coordinates (X, Y, Z) and a fluid particle, B and b,
for accurate reporting and does not constitute an official expressed in the inertial and noninertial coordinates (X,
endorsement, either expressed or implied, of such products or
Y, Z and x, y, z), respectively. At any instant in time, it
manufacturers by the National Aeronautics and Space
Administration. is very easy to express the relation of the position of a
2
American Institute of Aeronautics and Astronautics
pointin bothcoordinatesystemsby additionof the _term is zero because the rotation rate is assumed
vectors.
to be constant in this steady-state tbrmulation.
B=C+b (I) dB dd all; j(_o×b)
- _----+ (5)
The next step is to find the relationship between the dt dt dt dt
instantaneous velocity of a point expressed in both
coordinate systems. The velocity is found by b =d +_;+_x/;+mx/,+exb +ex(mx/,) (6)
differentiating the expression for the vector relation of
position, Eq. (I), with respect to time. Note that there
will be an added complexity to computing the 11_=_;+b+ 2rex/; +e0x(mxb) (7)
derivative of any vector quantity expressed in the
Now the acceleration of the origin of the
noninertial frame, e. g., b, because the unit normal
vectors (i,j, k) are changing as a function of time, due noninertial, grid-fixed reference frame, C, must be
to rotation. To find the derivative of b, the product rule sought. In this formulation, the noninertial reference
is used on the multiplication of the scalar components frame (Fig. 2) is following a curved path (the dashed
(x, y, :) and unit vectors (i,j, k). The relation between arrow) through inertial space as it simultaneously
the instantaneous derivatives of the unit vectors translates (_') and rotates (co). The origin of the
and mxb is found by taking the limit on the derivative noninertial reference frame must accelerate to follow
as dt goes to zero. I0 this curved path. The expression for the acceleration of
a grid that is moving in a curved path with constant
dB dC db
- f-- (2) rotation rate is
dt dt dt
d =_xC', _ =r..ox-u (8)
B=C'+b+mxb (3)
Note that this reference frame origin acceleration is
Now that the velocity relationship, 'Eq. (3), has zero when _' is parallel to m.
been derived, the far-field boundary conditions of the Now that the _;term is known, the expressions for
noninertial CFD problem can be discussed. For free
the difference between the accelerations computed in
stream boundary conditions, the fluid particles are at
the inertial frame and the noninertial frame (CFD grid)
rest in the inertial frame: therefore their velocity, B , can be completed. This difference in acceleration is
computed by subtracting the acceleration of a fluid
is zero. The expression for the fluid velocity at the
particle in the noninertial frame, b, from the
CFD grid free stream boundary conditions is /_ .The
acceleration of the same particle expressed in the
aircraft, the nonincrtial reference frame, and the CFD
inertial frame, B. By accounting for this difference in
grid are all translating together at a velocity, C, which
acceleration (pseudo-acceleration, B-b), equations
is negative the free stream velocity, u_. The boundary describing the motion of particles measured in a
values for B and _'are substituted into Eq. (3), noninertial frame can correctly mimic the total
acceleration of these fluid particles in the inertial
which is rearranged to form Eq. (4).
frame.
b==u - _xo L (4)
B-b" = cox-u + 2elxb +mx(olxb) (9)
Therefore, the free stream boundary conditions can
The goal is to model constant-rate CFD grid
be described as the combination of a uniform flow
rotation and translation with a steady-state CFD
component, u_, and a rigid body rotation component, calculation. The difference in acceleration between an
¢.oxbh .The CFL3D boundary conditions at the near- inertial and a noninertial frame of reference is
field or solid surface boundary conditions are not employed to form a source term correction to the N-S
affected inthis noninertial formulation. equations in CFL3D. To illustrate the formulation of
The expression lbr acceleration is computed by the source term, the existing implementation of the
differentiating the velocity relation, Eq. (3). mThe time stationary grid N-S equations in CFL3D must be
examined. In the CFL3D code, the N-S equations are
derivatives of the b and /; terms are determined in the
expressed in a regular-spaced, Cartesian grid
same fashion as the derivative of the b term was
coordinate system. The generalized grid coordinate
derived in Eq. (3). This relation is valid ['or any vector
system that defines the problem (x, y, :) is internally
quantity expressed in a noninertial frame. Note that
mapped by CFL3D to this regular-spaced Cartesian
3
American Institute of Aeronautics and Astronautics
grid(_,
q, 4) with a coordinate transform. The N-S The N-S equations in CFL3D are written in
equations are written in the regular-spaced, Cartesian conservative fi_rm, so the vector source term for the
grid coordinate systcm as_ four momentum and energy equations must include the
volume of each computational cell (J-') and the
a0
+ I- -0 (10) density, p, of the fluid. The energy equation source
at 03_ Or/ 03_"
term is the dot product of the local flow velocity, /_,
where Q is a vector of the conserved variables. The with pseudo-acceleration,
conserved variables are a combination of density, p, 0
velocity components (u, v, w) and total energy per unit
volume, e. The vector {_ is the conserved variables (mx(mxb)+ 2mx/_-mxu ).,.
divided by J. s=P (mx(mxbl+2o_x[_-coxu_), (17)
J
e^=7Q=ffl [ppupvpwe] r (11) (¢0x(_xb)+ 2mxb-mxu ):
_;.(mx(mxb)+ 2¢oxb-mxu )
The Jacobean (J) of the coordinate transformation from
the Cartesian to thegeneralized coordinate system is
where b=[x 3' z] and /_=[u v w]. For
j _ (12) additional papers on the source term equations and
a(.,-,>,=,) associated physics of similar applications of this
theory, see Refs. 2-5. Refs. 4 and 5 also include
The inviscid flux terms are F. G, and H and the reference frame angular and translation acceleration
viscous flux terms are F,. G_.,and H,.. The /0 , _,, /._, terms.
. G,, and /'t, flux terms are created by dividing by
ADIFOR Automatic Differentiation
J in the same manner as Q was. For a nondeforming Automatic differentiation is a technique for
mesh (J is constant with respect to time), the solution
augmenting computer programs with statements for the
is advanced in time with the residual, R.
computation of derivatives. This technique relies on
laQ the fact that every function, no matter how
- R(Q ) (13) complicated, is executed on a computer as a
J 03t
(potentially very long) sequence of elementary
The residual is computed as operations such as additions, multiplications, and
elementary functions (e. g., sine and cosine). By
repeatedly applying the chain rule of differential
at(} -(}, ) + (14)
R(Q) =- 03( P' ) _ Or/ 03¢ calculus to the composition of those elementary
operations, derivative information can be computed
To permit noninertial calculations, a source term exactly and ina completely automated fashion.
(S) is added to the standard CFL3D residual The ADIFOR process is a technique that applies
calculation. the chain rule of differentiation to propagate, equation
by equation, derivatives of intermediate variables with
] 03Q
respect to the input variables. The ADIFOR tool has
- R(Q)+ S (15)
J 0t been developed jointly by the Center for Research on
Parallel Computation at Rice University and the
The source term is a vector with 4 nonzero
Mathematics and Computer Sciences Division at
components (the continuity equation is not affected). Argonne National Laboratory. In general, to apply
ADIFOR to a given F77 code, the user is only required
S =PEos , S, S=S] r (16) to specify those program variable names that
correspond to the independent and dependent variables
The momentum equation source terms (&, S,, and of the target differentiation. The ADIFOR tool then
S:) are the three components of the pseudo-acceleration determines the variables that require associated
derivative computations, lbrmulates the appropriate
(B -b ). The "work" done by these momentum sourcc
derivative expressions, and generates new F77 code for
terms must be included in the energy equation. The
the computation of both the original simulation and the
energy equation source term (&) is the dot product of
specified derivatives.
the local velocity with the pseudo-acceleration.
4
American Institute of Aeronautics and Astronautics
ThemodifiedversioonfCFL3Disabletocompute These 2-D NACA0012 cases shown in Figs. 3 and
theaircraftforcesandmomentassafunctionofthe 4 were chosen for initial validation of CFL3D.N1.AD.
threebody-axoisrthogoncaolmponen(pts,q, and r) of To improve convergence, a blend of half standard
the rotation vector. The application of ADIFOR to the CFL3D and half CFL3D Iow-Mach-number
modified version of CFL3D produced a code that preconditioning was applied for the 0.1 Mach (Figs. 3a
computed the forward mode derivatives of the aircraft and 4a) case. This preconditioning option was not
body-axis force and moment coefficients (CN, CA, CS, applied to the 0.5 Mach (Figs. 3b and 4b) or 0.8 Mach
CI, Cm, and Cn) with respect to the three body-axis (Figs. 3c and 4c) cases. Note that the CFL3D.NI.AD
rotation rates (p, q, and r) and the flow angles of attack derivative values are in excellent agreement with the
(c_)and sideslip (13). SFLOW values. For the O.I and 0.5 Mach cases, the
The latest beta version of CFL3D employs dynamic differences between CFL3D.NI.AD and SFLOW,
memory allocation and MPI libraries for ease of use although small, are most likely due to the formulation
and efficient, scalable parallelization. These differences between the flow solvers in CFL3D.NI.AD
implementations are not standard F77 features, and and SFLOW. The SFLOW code employs the
therefore previous releases of ADIFOR cannot handle hand-coded sensitivity equation technique and an
the code without manual preprocessing and unstructured grid discretization, whereas
postprocessing. The latest release of ADIFOR 3.0_uhas CFL3D.NI.AD is an automatically differentiated
reduced or eliminated much of the manual processing structured grid formulation.
associated with the MPI libraries: techniques for The 0.8 Mach case (Figs. 3c and 4c) shows poor
handling the dynamic memory allocation libraries are convergence properties. The poor convergence of
being developed. CFL3D.NI.AD at 0.8 Mach may be due to the
interaction of a shock, the flux limiter implemented in
CFL3D, and the automatic differentiation technique.
Examples and Results
The CFL3D smooth flux limiter Ituned to _;= I/3 was
The two examples in this study were a 2-D Euler
employed for the NACA0012 study. This poor
study of NACA0012 airfoil and a 3-D turbulent N-S
calculation on the ICE configuration 6 (Fig. 1). The convergence may be due to the automatic
differentiation technique attempting to formulate the
NACA0012 study will be detailed first, because it was
continuous derivative of a shock and flux limiter,
used for initial validation by comparisons to existing
methods. which does not have a continuous derivative. The 0.8
Mach case is also the worst comparison to SFLOW.
Only the final value for SFLOW is quoted in Ref. 5:
NACAO012 Dynamic Derivatives
The NACA0012 study focused on the effect of therefore the SFLOW 0.8 Mach case may or may not
pitch rate on the coefficients of normal force (Fig. 3) be fully converged. The convergence was not
and pitching moment (Fig. 4) at zero deg angle of improved by disabling multigrid calculations or
attack and zero pitch rate. The derivatives (CNq, Fig. 3, performing additional iteration cycles. At 0.8 Mach,
the final value of CFL3D.NI.AD and SFLOW differ in
and Cm,, Fig. 4) of these three and moment
normal force pitch rate derivative by 4.4_ (Fig. 3c)
coefficients with respect to nondimensiona[ pitch rate
were computed by the ADIFOR-generated, noninertial and in pitching moment pitch rate derivative by 8.9%
CFL3D code (CFL3D.NI.AD). The pitch rate (Fig. 4c).
derivatives are nondimensionalized by dividing by the
ICE Pathlines at Zero Rotational Rate
airfoil chord and multiplying by two times the free
Alter CFL3D.NI.AD was validated by comparison
stream velocity. The NACA0012 pitching moment
center is located at the leading edge of the airfoil. The with SFLOW, the ICE configuration (Fig. 1) flow
structure was examined with pathlines. These pathlines
convergence history of the derivative values is shown
in Figs. 3 and 4. The discontinuities in the derivative are shown to illustrate the changes in airflow structure
with increasing angles of attack. Pathlines for the
convergence history are due to mesh sequencing from
a coarser to a finer mesh every 500 iterations. A starboard half-span of the ICE configuration are shown
maximum of three levels of multigrid was employed in Fig. 5. The pathlines were seeded slightly ahead of
the sharp leading edge (just outside the boundary
on the finer meshes. The 2-D grid dimensions (49 ×
layer). The pathlines were computed from a full-span
13, 97 z 25, 193 x 49, and 385 x 97) are denoted for
N-S CFL3D solution on a grid with approximately 3
each mesh sequencing level. The derivative values are
million cells. These symmetric solutions in Fig. 5 were
compared to results computed by a similar method
calculated at 0.6 Mach, zero deg angle of sideslip, zero
published by Limache and Cliff (SFLOW), -_a panel
rotational rate, and various angles of attack. All ICE
method (QUADPAN), s2 and a vortex lattice method
CFL3D solutions in this study were computed with the
(VORLAX).I_
S-A turbulence model at a Reynolds number of
5
American Institute of Aeronautics and Astronautics
2,490,000per foot or 71,760,000per mean
(Fig. 6c) which resulted in static longitudinal
aerodynamcichord.Threelevelsof multigridwere
instability abovc 15 deg angle of attack. Note that
usedforthefinegridsolutionofthe0-15degangleof
CFL3D captures the radical change in Cm measured by
attackcasesandmultigridwasdisabledforthefine
WT (Fig. 6c).
gridsolutionofthe20-30degangleofattacckases.
Figure 7 shows the comparison among the lateral
Thestructuroef thesymmetrifclowdepictedin
angle of sideslip derivatives tbr three
Fig.5 aidsinterpretatioonf thesubsequefnigtures, central-finite-difference estimates from wind tunnel
whichdepicttbrce,momenat,ndstabilityderivative
data 6 (CD-WT) and an ADIFOR-generated CFL3D
informationat theseflight conditionsN. otethe
solution (CFL3D.AD). The CFL3D.AD derivatives are
attachefldowat5degangleofattack(Fig.5a).Weak dashed lines and the CD-WT derivatives are the
leadinegdgevorticafllowwaspresenatt10degangle
symbols with a central-finite-difference step of +2, +4,
ofattack(Fig5b).Thisinitialleadingedgevortex
and +6 deg angle of sideslip for the circle, square, and
structurgeainedstrengtaht15degangleofattack(Fig.
diamond, respectively. The +2 deg CD-WT data is
5c).Avortexburstdevelopendeatrhetrailingedgeat
connected with solid lines because the smallest
20degangleof attack(Fig.5d).Thisvortexburst
central-finite-difference step (+2) is presumed to be the
structuriesidentifiedbyanabrupsttreamwisinecrease
most accurate of the three finite difference step sizes
invortexdiameterT.heinitialvortexburststructure
for small sideslip disturbances. All three finite
intensitieadndmovedtbrwardat25and30degangles
difference step sizes are shown to give an indication of
ofattack(Fig.5eand5f).
the nonlinearities or measurement noise in the wind
tunnel data. The derivatives in Fig. 7 are presented in
ICE Forces, Moments, and Lateral Derivatives
the units of deg -_.
A summary of the figures depicting ICE forces,
The effects of the vortical flow structure (Fig. 5)
moments, and stability derivatives in this study is
can be seen clearly in the lateral force and moments
given in Table 1. Note that the elements assumed zero
angle-of-sideslip derivatives (Fig. 7). The initiation
would become significant at a nonzero angle of
and strengthening of vortical flow between 5 and 10
sideslip, roll rate, or yaw rate. The ICE angle-of-attack
deg angles of attack (Fig 5a and 5b) can be interpreted
derivatives are not presented in this study, but are to have sharply influenced the angle-of-attack trends of
presented in Ref. 6 tor 0-10 deg angles of attack. The
CS_ and Cnl_ (Fig, 7a and 7c) computed by CD-WT
orientation of the six lorces and moments, two flow
and CFL3D.AD. Then, the derivatives CS_ and Cnl3
angles, and three body-axis rotational rates is shown in
(Fig. 7a and 7c) dramatically reversed angle-of-attack
Fig. I. Figure 6 shows a comparison of longitudinal
trends above 10 deg angle of attack. The CFL3D.AD
forces and moment at 0.6 Mach, zero deg angle of
CI_ (Fig. 7b) dcrivative showed excellent agreement
sideslip, and zero rotational rate. The moment
with CD-WT for 0 to 15 deg angles of attack. The Cl_
reference center of the ICE is located longitudinally at
comparison deteriorated at higher (20-30 deg) angles
39% of the mean aerodynamic chord (MAC) and a
of attack.
distance of 16%, of the MAC below the body. Thc
As angle of sideslip varies, each wing experiences
comparison is the wind tunnel data (solid line, WT)
different effective leading-edge sweep angles. Due to
present in the ICE simulator database _ and CFL3D
the highly swept (65 deg) leading edge of the ICE
(dashed line). Coefficients of normal tbrce, axial force,
configuration the vortical flow field over the wing may
and pitching moment are shown in Fig. 6a, 6b and 6c,
be sensitive to changes in effective leading-edge sweep
respectively.
angle. Therefore, the calculation of a vortex burst
Figure 6 shows good agreement between WT and
CFL3D. The WT data has more detail because it was structure that tbrmed symmetrically at 20 deg angle of
attack (Fig. 5d) may be produced asymmetrically at
measured at approximately 1 deg increments, which
lower (10-15 deg) angles of attack. An asymmetric,
were smaller than the 5 deg increments of the CFL3D
calculations. There is no flow visualization information bursting vortex structure may have been responsible
for the dramatically reversed angle-of-attack trends in
available lor the WT data, but the CFL3D pathlines
will be used to infer the effects of flow structure on the the lateral derivatives (Fig. 7). The ICE configuration
does not have any vertical surfaces, so the magnitude
CFL3D calculations, which may also indicate the flow
structure effects on WT measurements. Note that the of CS_ and Cnl_ (Fig. 7a and 7c) was reduced as
compared to a configuration with vertical surfaces. The
initiation and strengthening of vortical flow between 5
and 15 deg angles of attack (Fig. 5a-5c) increased the small magnitude of CS_ and Cnl3 may have hindered
measurement accuracy and exacerbated comparison of
normal force (Fig. 6a). The increasing strength of the
CD-WT with CFL3D.AD.
vortex flow and the forward movement of the burst
location over the wing between 20 and 30 deg angles
of attack (Fig. 5d-5t), increased the pitching moment
6
American Institute of Aeronautics and Astronautics
ICE Dynamic Derivatives of attack and lower. These trends became less
Figures 8 and 9 show the dynamic derivatives consistent at 20, 25, and 30 deg angles of attack, which
computed by the DYNAMIC H code and corresponded with the initial indication of a symmetric
CFL3D.NI.AD. The DYNAMIC code utilized strip vortex burst structure in Fig. 5d.
theory and the results of the high-angle-of-attack The CFL3D.NI.AD differentiated flow solver had
stability and control prediction code HASC 15'l(' to convergence difficulties at 20, 25, and 30 deg angles of
calculate the dynamic derivatives. The HASC code attack. The 30 deg angle of attack case never reached a
employs VORLAX 13 and empirical corrections to steady-state value, so an average of the last 2 thousand
predict configuration lbrces and moments at various iterations is presented. These convergence difficulties
flow angles and rotational rates. The derivatives were may have been due to the presence of bursting vortex
computed at zero rotational rate, zero angle of sideslip, structures, with their inherent unsteadiness and
and various angles of attack. The CFL3D.NI.AD increased sensitivity to disturbances. These high
dynamic derivatives were computed assuming angle-of-attack conditions may be more suitable to a
rotations about the moment center of the configuration, time-accurate solution, but in the interest of
which is located slightly below the body. The minimizing computational resource requirements, that
longitudinal pitch rate derivatives CN,_ and Cm,_ are approach was not attempted in this study.
shown in Fig. 8a and 8d, respectively. The longitudinal
dynamic derivatives were nondimensionalized by ICE Pathlines at Nonzero Rotational Rates
dividing by the mean aerodynamic chord (345 in.) and Figure 10 shows the ICE configuration at 0.3
multiplying by two times the free stream velocity. Mach, 15 deg angle of attack and zero deg angle of
The rolling moment dynamic derivatives Clj, and sideslip, performing velocity vector rolls at various
Clr are shown in Fig. 8b and 8e, respectively. The rotational rates. In these velocity vector rolls, the
yawing moment dynamic derivatives Cn_,and Cnr are rotation vector was parallel to the free stream velocity
shown in Fig. 8c and 8f, respectively. The side force vector; this condition simulated a wind tunnel
dynamic derivatives CSi, and CSr are shown in Fig. 9a rotary-balance test. These solutions were computed by
and 9b, respectively. The lateral dynamic derivatives the noninertial, modified CFL3D (CFL3D.NI) code.
were nondimensionalized by dividing by the wingspan The rotational rate (f_) was nondimensionalized by
b (450 in.) and multiplying by two times the free multiplying by the wingspan b (450 in.) and dividing
stream velocity. There was no forced, oscillatory by two times the free stream velocity (u_), with a
motion wind tunnel data for comparison. positive rotational rate indicating the starboard wing
Both codes, CFL3D.NI.AD and DYNAMIC, was descending. The 0.2 and 0.4 rotational rate cases
showed fairly good comparison. Both of the (Fig. 10b and 10c) showed a much tighter vortex core
CFL3D.NI.AD pitch rate (q) derivatives (Fig. 8a and on the ascending, port wing than the descending,
8d) show a local maximum or a minimum near 5 deg starboard wing. The 0.4 rotational rate case (Fig. 10c)
angle of attack. Note that the CFL3D.NI.AD depicted a vortex burst on the descending, starboard
calculation of Cmq (Fig. 8d) was consistently more wing. From this point of view, the vortex wakes in Fig.
negative than the combined analytical and vortex 10b and 10c appear to be converging, but actually were
lattice method of DYNAMIC and VORLAX. This spiraling around the rotation vector.
trend agrees with those of both SFLOW and
CFL3D.NI.AD when compared to VORLAX for the ICE Rotary-Balance and Noninertial CFL3D
2-D NACA0012 case (Fig. 4). Comparison
The roll rate (p) derivatives (Figs. 8b+ 8c, and 9a) Figure 11 shows a comparison of wind tunnel
also showed a reversal of angle of attack trends at 5 rotary-balance data6 (ROT-BAL, solid line) and
deg angle of attack. The reversals of the q and p CFL3D.NI (dashed line) at 0.3 Mach, 15 deg angle of
derivative trends at 5 deg angle of attack corresponded attack, and zero deg angle of sideslip. Mach 0.3 was
to the indication of vortical flow at 10 deg angle of chosen to simulate the incompressible conditions of the
attack in Fig. 5b. These p derivatives also showed low-speed, rotary-balance tests. The ICE configuration
another local extreme at 15-20 deg angle of attack, was rotated about the moment center of the
which was slightly below the indication of vortex configuration, which is located slightly below the
bursting in the static pathlines (Fig. 5d). A roll rate body. The figure shows the change (A) in force or
creates differential angles of attack on each wing, moment coefficient between cases nonrotating and
which may induce asymmetric vortical burst structures rotating about the velocity vector. The rotation rate (_2)
at lower angles of attack than a zero-roll-rate about the velocity vector was nondimensionalized by
symmetric case. multiplying by the wingspan b (450 in.) and dividing
The yaw rate (r) derivatives (Figs. 8e, 8f, and 9b) by two times the free stream velocity (u=), with a
had consistent trends in angle of attack at 15 deg angle positive rotational rate indicating the starboard wing
7
American Institute of Aeronautics and Astronautics
wasdescendingN.otethatnonlineaerffectswith
RAM. By means of a batch queuing system, the
rotationaralteweremodeleidnCFL3D.NI.
16-processor O2K total wall time was achieved
Theincreasinenormaflorce(ACNF, ig.Ila)due through multiple 45 rain runs. The 16 processor O2K
torotationwasverysimilarbetweeRnOT-BALand had significant shutdown and restart overhead
CFL3DNI.Thechangienaxialforce(ACAF,ig.IIb) (approximately 10%), which adversely affects total
duetorotationwasverysimilarinmagnitudbeetween wall time for the CFL3D.NI.AD examples.
ROT-BALandCFL3D.N1b,utoppositiensign.The Note that CFL3D.NI required 0.5 hour (3.8c_)
changein pitchingmomen(tACreF,ig.llc) dueto more execution time than the original CFL3D
rotationwasassumetdo beanevenfunctionas steady-state execution wall time for the ICE
calculatedby CFL3D.NI,but ROT-BALshowed configuration with the S-A turbulence model. The
inconclusivtreendsT.hechangeinsideforce(ACS, corresponding wall time increase lot 2-D and 3-D
Fig.lld)duetorotationwasassumetdobeanodd Euler calculations due to noninertial modifications was
functionascalculatebdyCFL3D.NIb,utROT-BAL approximately 15%. The noninertial modifications had
showedinconclusivterendsT. hechangein rolling a larger penalty for Euler than turbulent N-S solutions
momen(tACI,Fig.lle)duetorotationwasthebest because N-S and S-A solutions required more
lateraclomparisoonfCFL3D.NwIithROT-BALT.he calculations per iteration than Euler solutions. The
changein yawingmomen(tACn,Fig.Ilf) dueto increased calculations per iteration of the turbulent N-S
rotationcalculatebdyCFL3D.NwIasmuchgreateirn solution masked the same number of noninertial
magnitudaendoppositiensignoftheROT-BALtrend. modification calculations per iteration of the turbulent
The nonlineareffectswith rotationalrateas N-S and Euler solutions.
calculatebdyCFL3D.NinI Fig.11canbecorrelatetod A time-accurate CFL3D solution that would
thecalculatioonf a vortexburststructuroeverthe emulate a CFL3D.NI solution was estimated to require
descendiwngingasillustrateindFig.10band10cT.he approximately 175 hours, or more than an order of
differencebetweenROT-BALand CFL3D.NI magnitude increase in wall time over a CFL3D.NI
(highlighteidnFigs.lib, llc,andIlf) iscurrently calculation. The central-finite-difference estimate wall
underinvestigationA. possibleexplanationof time was calculated by multiplying the CFL3D.NI time
ROT-BALasymmetrimesaybemodealsymmetrioers by II (one function plus ten perturbed solutions) to
installatiomnisalignmenTtsh.epoorcomparisoonfs yield 148.5 hours, which was scaled between the two
ROT-BALandCFL3D.NmI aybedueto rotation 02K computers assuming perfect, linear speedup with
aboutdifferentlocationstbr theexperimentaanl d a ratio of 3 worker processors to 13 worker processors.
computationcaalsesT.heCFL3D.NcIodesimulated In other words, 13.5 × 11= 148.5 and 148.5 × 3/ 13=
rotationaboutthereportedmomenctenterof the 34. The central-difference estimate required 9.7%
configurationw,hichis outsidethe model.The more wall time than CFL3D.NI.AD between 0 and 15
ROT-BALtestsmayor maynothaverotatedthe deg angles of attack. Compared to the 0-15 deg angle
modealboutthatmomenctentelrocation. of attack solutions, CFL3D.NI.AD required three to
four times the wall time at 20, 25, and 30 deg angle of
attack, due to differentiated flow solver convergence
Timing
difficulties. The vortex burst structures at the higher
Table 2 describes the processors, wall time, and
(20-30 deg) angles of attack (Fig. 5d-5f) may have
RAM required by the original CFL3D, CFL3D.NI, and
been responsible for the convergence difficulties.
CFL3D.NI.AD. The column labeled "Independents"
indicates whether function only (zero independents) or
Conclusions
function plus derivatives with respect to angle of
An initial application of ADIFOR to CFL3D with
attack, angle of sideslip, roll rate, pitch rate, and yaw
constant-rate noninertial modifications to compute
rate (five independents) were calculated. The column
constant-rate rotary stability derivatives was
labeled "'Processors" indicates the number of SGI
completed. This application was validated for a 2-D
Origin 2000 TM (O2K) processors employed lor the
NACA0012 Euler case by comparison to the SFLOW
calculations. All three parallel versions of the CFL3D
code, a similar formulation. ADIFOR-generated
code employed in this study use one of the processors
noninertial CFL3D derivatives of a 2-D NACA0012
lor administrative tasks, so the number of actual
airfoil showed good comparison with existing methods
computing processors is one less than the number
at 0.1 and 0.5 Mach. Symmetric vortical flow
quoted in the "Processors" column. The four-processor
structures for the ICE configuration were identified by
runs were performed on a NASA Langley
means of computational flow visualizations of
Multidisciplinary Optimization Branch four-processor
turbulent N-S calculations at 5-30 deg angles of attack.
O2K with 4 Gb RAM. The 14-processor runs were
The nature of these vortical flow structures was
performed on a HPCCP 16-processor O2K with 12Gb
correlated to the behavior of forces, moments,
8
American Institute of Aeronautics and Astronautics