NASA Technical Reports Server (NTRS) 19970006701: Application of a Third Order Upwind Scheme to Viscous Flow over Clean and Iced Wings

10 Pages·0.45 MB·English

by  NASA Technical Reports Server (NTRS)

/"
/
I
_,j/' ....
NASA-CR-202616
AIAA 94-0485
Application of a Third Order Upwind Scheme
to Viscous Flow over Clean and Iced Wings
A. Bangalore, N. Phaengsook and L. N. Sankar
School of Aerospace Engineering
Georgia Institute of Technology
Atlanta, GA
32nd Aerospace Sciences
Meeting & Exhibit
January 10-13, 1994 / Reno, NV
For permission to copy or republish, contact the Amedcsn Institute of Aeronautics and Astronautics
370 LT.nfant Promenade, S.W., Washington, D.C. 20024
APPLICATION OF ATHIRD ORDERUPWIND SCHEME TO
VISCOUS FLOW OVERCLEAN AND ICED WINGS
A. Bangalore 1,N.Phaengsook I and L N. Sinker2
School of Aerospace Engineering
Georgia Institute of Technology, Atlanta, GA 30332
effects of turbulence. Rnaily, for the sake of computational
ABSTRACT
stability and efficiency, these methods used an implicit
time-marching algorithm that permit the use of very large
time steps. Reliable convergence to steady state in less
A 3-D compressible Navier-Stokes solver has been
than 1CPU hour ofcomputer time on a Cray Y/MP classof
developed and applied to 3-D viscous flowover clean and
computer, on • 100,000 point gdd. Reasonable correlation
icedwings.This method uses athird order accurate finite was observed between the calculations and the measured
volume scheme with flux difference splitting to model the data obtained byBragg and his=)workers [Ref. 4].
invtscld fluxes, and second order accurate symmetric
differences to model the viscous terms. The effects of
As these analyses were applied to a number of
turbulenceare modeled using ak-¢model. Inthe vicinityof icedwingcalculations, two drawbacks ofthese approaches
the solid walls, the k and = values are modeled using became apparent to the present investigators. Firstof all,
Gorski'salgebraic model. Sample resultsare presented for
the turbulence model used in these analyses is a simple
eurface pressure distributions, for untapered swept clean
algebraic eddy viscosity model, which uses the boundary
and iced wings made of NACA 0012 airfoil sections. The
layer thickness f> (or an equivalent quantity such as the
leading edge of these sections ere modified using a
distance from the wall, y) as the length scale of turbulent
simulated ice shape. Comparisons with the experimental
eddies. The algebraic model used the boundary layer edge
data obtained byProf. Bragg stthe University of Illinoisare
velocity U(x) (or an equivalent quantity y I¢olmax) as a
given.
measure of the velocity scale of turbulent eddies. Such a
INTRODUCTION simple representation of turbulence works well only for
simple attached or mildly separated flows. In the case of
iced wings, the flow field contains boundary layer overthe
The aerodynamic characteristics of liftingsurfaces
solid surfaces, reversed flow regions, and a free shear
such as wings and horizontal tails may be dramatically
layer coming off sharp corners and turns associated with
altered by the accumulation of ice over the leading edge.
the iceshape. The turbulent length and velocity scales for
Even l relatively small amount of ice accumulation may
free shear layers are clearly different from the wall
lead to dramatic rise in pressure drag, loss in lift, and in
bounded flows. Use of a simple eddy viscosity model to
tome cases, premature stall
such complex flows lead to inaccurate results inregions of
separation. For example, the length of the separation
Under the support of the NASA Lewis Research
bubble was often inaccurately predicted. Potapczuk has
Center, a research effort has been underway at Georgia
discussed the limitaUonsof the Baidwin-Lomax model, and
Tech on the effects of icing on the aerodynamic
has recommended phenomenological fixes to this model
pedormance ofliftingsurfaces. Both2-D and 3-D analyses
[Ref. 5].
have been developed, and calibrated [Ref. 1-3]. These
analyses used relatively simple algorithms fo sofve the
The second drawback of the previous analyses is
compressible Navler-Stokes equations. For example,
the irrecoverable filteringout and lossof information bythe
standard second order accurate central differences are
low-pass filters. These filters tend to smear out shock
usedtomodel the derivatives inthe governing equations. A
discontinuities, contact discontinuities, and shear layers.
fourth order low-pass filter is used at every time step for
The rate at which the filtering is applied is proportional to
removing any high frequency spatial oscillations. These
iul+a, where lulis the magnitude of localflow, and ais the
analyses also used simple algebraic turbulence models
speed ofsound.
such as the Baldwin-Lomax model to account for the
The present work Improves upon the existing
1Graduate Research Amiatant
numerical analysis, presented in Ref. 3, and attempts to
2 ProfessorofAerospace Engineering, Senior Member rectify the above two drawbacks. It uses aflux difference
AtAA.
echeme, which dynamically and adaptively reduces the
Copyright© 1994 by L N. Sinker, Publishedbythe smearing of contact discontinuities, and vorticity waves,
Amedcan InstituteofAeronsutics andAstronautics, Inc.
and leads to crisp resolution of free shear layer,; that arise
pem n.
in the flow [Fief. 6]. This method alsouses a sophisticated
k-¢ model to account for the effects of turbulence. Inthe where the operators 6! etc. represent standard
vicinity of solid wells, this model uses a set of
central diflerences suchIs
phenomenologically derived boundary conditionsfor kand
¢, Is described in Ref. 7. The improved analysis allows
transition to be detected automatically by the solver as 61_ . ,.!,= ,-I,2
m_lons ofvery lowturbulent Idnaticenergy.
(4)
The Improved eolver hIs been applied to clean
andloed, swept wing configurations tested by Prof. Bragg
The numerical flux vector _ differs from the
end hisoowod(em atthe University of Illinois,documented
InthePh.D. dissertation ofKhodadoust [Ref. 8]. physical fluxvector F in • manner specified by Roe [Ref.
6]. Let qL andqR be two estimates tothe flow field vector
The rest ofthispaper isorganized Is follows. First, q, to the left side end right side of the cell face (i+l/2,j,k).
the Roe scheme is briefly outlined. Next, the k-¢ Then,
formulation is described. Finally, several sample
applications of the improved solver to iced wing analysis
are presented I: ,.,,, . [;(q,) + ::(qR) +ITIAIT-'(q_ -q_)
2 2
MATHEMATICAL AND NUMERICAL FORMULATION
where,
The 3-D compressible Nsvier-Stokes equations
are giveninaCartesian coordinate system as q.., - q, _ q,.2 - q..,
3 6
...=.aa+0FaG aH 0R aS 0T
at 8x ay & _ _ &
(I)
Here q isa vector of conserved flowproperties; F, qL = q, + q,-1- q, + q, - q,-i
G and H are the inviscid flux vectors; R, S and T ere 3 6
viscous terms, which include the effects of turbulence {5)
_lrough the eddy viscosityconcept. Hero T is a matrix that contains the left
The calculations are carried out in • curvilinear eigenvectom of the matrix aF/aq, and A contains the
body-fitted coordinate system (F,,n,f;).In thissystem, the eigenvalues of this matrix. Those quantities must be
governing equations may bewritten as computedusingspecial *Roe" averages of qLandqR-
It may be shown that the above estimates of qL
and qR lead to s spatiallythird order aocurate formulation,
ot a_ an a_ a_ a_ o_
on uniformgridsandon mildlystretched grids.
(2)
Equation (3) represents a system of non-linear
where _ etc. ere related to their Cartesian Idgebralc equations forthe flow properties vectorq sttime
level 'rH-l'. Furthermore, each node (i,j+l,k) iscoupled to
counterparts bythe metrics oftransformation.
its nine neighbors. Thus, this system of equations are
The Roe Scheme: highly coupled. The classical Beam-Warming approximate
fectorization scheme was used solve these coupled
nonlinear equations, in• manner similar to that described
The objective of the lime-marching algorithm is to
inReference 1.
advance the flow field at cells (or nodes (i,j,k) from • lime
level 'n'tothe nexttime level 'n+l'. Forthis purpose, Inthe
K-t Turbulence Model
preeent work, equation (2) was discmtized Is follows:
The starting point is the Bou=mesq assumption for
the turbulent mtres_s.
+ "|" _.j.t + "n',,j.l + "_,"%l.k "
At
..% 2 ÷
(3) ¢3Xj aXk
2
The tud)ulent vls(:os_ Pl has to modeled. This is turbulence field, and isnot needed during subsequent time
steps.
achieved by modeling two physically conceivable
quantities, the turbulent kinetic energy k, and the
dissipation rate e. The basic high Reynolds number K-E Near solid walls, the classical high Reynolds
model integrates the followingtransportequations for kand number k-Emodel breaks down. In the present work, near
the solid waU, the k and c were assumed to have the
IL
followingvariation:
Dk __a (l_t_ ak +p_
k.Ay 3
Ok) ix
z= constant
Herer the constant A was evaluted using the
D,
computed values of k eswal nodes away from the solid
surface.
RESULTS AND DISCUSSION
The turbulent production P,isgiven by:
Nsvler-Stokes calculations were carded out using
the flowsolver developed inthe previous section, for swept
and unswept, untapered wing-alone geometries tested by
Bragget at. The airfoil sections are made of NACA 0012
aldoil sections. The wing tested had a chord length of 15
inches and a semispan of 37.25 inches. To study the
andthe turbulent viscosity _ is related to kand Eby"
effects of icing, the leading edge was modified with a
simulated ice shape. This ice shape corresponds to that
k2 measured atthe NASA Lewis Icing Research Tunnel. Iors
Ixt
NACA 0012 airfoil, st a freestream velocity of 130 mph,
£
angle of attack of 4 degrees, icing time of 4 minutes,
volume median dropletdiameter of20 microns, liquidwater
The basic k-zmodel constants am as follows.
content of2.1.grams per cubic meter and • temperature of
18 degrees. Bragg's measurements were done at zero
Cl_ -.09
sweep andata30 degree sweep.
C] - 1.43
The calculations used a 121 x24 x45 grid, with91
C 2 - 1.92 points on the airfoil at each span station, 14 spanwise
stations over the wing, and 45 points in the direction
ok - 1.0 normal to the wing surface. Roughly s quarter of these
pointsare inthe boundary layer region.Thus, the boundary
0z - 1.3
layer is only moderately resolved in these preliminary
calculations. Figure 1shows the internallygenerated body-
As inthe case of the Navier-Stokes equations, the
fitted gridused inthe present simulations. Figure 2 shows
k-t equations were cast on a curvilinear coordinate system
a typicalairfoilcross section, forthe icedwing case.
Inconservation form, and solved using an ADI procedure.
The source terms appearing inthe kand Eequations were Calculations have been carded for at least 16 of
lagged by one time step for simplicity. The turbulence
the possible 32 combinations: swept and unswept, clean
model equations were solved separately at the end of
and iced, ingle of attack of 4 degrees and 8 degrees, Roe
every time step, after the mean flow properties had been scheme and central difference scheme, Baldwin-Lomax
computed. This resulted inthe solution of a 2 x 2 system
model and k-_model. For the sake ofbrevity, only s small
f_t may be efficientlysolved. subsetof_ calculations carried out are discussed herein
detail.
The starting values for k and z were obtained
muming the flow to be tn equilibrium. That is, the
Rgure 3 shows the surface pressure distribution
producllon and dissipation terms inthe k-¢equations were
for the clean wing at 8 degrees angle of attack st 3
let to be equal. This results in • system of simultaneous
apanwtse locations, 42%, 56%, 72% end 89%. Inthese
akjebr_c equations for k and ¢, which may be solved to calculations, "the Roe solver was used, and a Baldwin-
initialize the turbulence field. These algebraic equations
Lomax model was employed. Because of the increased
Involve the eddy viscosity, which was computed using the operations counts associated with the Roe scheme, the
Bak:lwin-Lomaxmodel. Note that the Baldwin-Lomax model
calculations were typically 30% more expen.Q;veper grid
was used onlyst the start ofthe calculationsto Initialize the
3
point, per time step, compared to the original central obtained for swept, clean and iced wing configurations.
difference scheme. The Roe so_ver was alsos'dff,requiring Them issome additional overhead associated withthe Roe
smaller time steps than the central difference scheme. scheme and the k-= model. Work is needed on
Neverlhele_, celculaUons similartothose shown on figure Improvements to the algorithm to reduce the CPU time
8 could be obtained withthe Roe/Baldwin-Lomax solver in requirements. Additional systematic validation of the
about 1500 Iterations, requiring about 90 minutes of CPU Improved flow solver is now underway, prior to its
Wne on a Cray Y/MP system usinga single processor. adaptation for icedwing performance analysis.
Figure 4 shows the Roe/k-t solver results tor the ACKNOWLEDGMENT
lame Case. The results are comparable to the Baldwin-
Lomax region,except near the leading edge region, where This work was supported by the NASA Lewis
the Cp dletribution shows a small plateau, Indicative of a Research Center under Garnt NAG-3-768. Dr. Mark
small amount of local separation. The k-t model behaved Potapczuk isthe technicalmonitor.
well. The strong source terms cause the k-t equation set to
be stiff. These equations were therefore solved using a REFERENCES
smaller (local) time step than the basic flow equations.
Upper and lower limits were placed on the kand ¢values, 1. Wu, Jiunn-Chi, "A Study of Unsteady Turbulent
in orderto prevent the solver Irom diverging atthe start of Flow past Airfoils," Ph.D. Dissertation, Georgia Instituteof
the calculations. The flow solver had no other user Technology, Atlanta, GA 1988.
adiIj_table constants, and was run usingthe standard k-=
constants described in the previous section. The final 2. Kwon, O. J., and Sankar, L. N., "Numerical Study
results for the k and t fields were well withinthese upper of the Effects of Icingon Finite Wing Aerodynamics," AIAA
lind lowerlimits.Thus, these limitsdo not affect the quality Paper 90-0757.
of the final, converged solution. A typical k-¢ solution
required about 2 hours of CPU time, on a Cray Y/tAP 3. Kwon, O. J. and Sankar, L.N., "Numerical Studyof
system, using s single processor. Work is now underway the Effects of icing on Fixed and Rotary Wing
on reducingU_; time requirement. Aerodynamics," AIAA paper 91-0662.
Figure 5 shows the calculations and 4. Bragg, M. B. and Khodadoust, A., "Effects of
measurements for the swept icedwing at 4degree angle of Simulated Glaze Ice on a Rectangular Wing," AIAA paper
attack. In this case, the flow separates over the leading 89-0750.
edge ice shape, and the pressure distribution forms a
plateau. The Rce/k-t solverisseen to performsatisfactorily 5. Potapczuk, M., "Navier-stokes Analysis of Airfoils
at allthe stations. Itisseen that the plateau over the upper with Leading Edge Ice Accretions," Ph. D. Dissertation,
surface is well captured, while the plateau over the lower The University ofAkron, Akron, Ohio, 1989.
surface is notcorrectly captured. The solutiondownstream
of the plateau region is in excellent agreement with the 6. Roe," P. L., "Approximate Riemann Solvers,
mea_urmllents. Parameter Vectors and Difference Schemes," Journal of
Computational Physics, Vol. 43, 1981, pC)357-372.
Rgure 6 shows some preliminary results for the
7. Gorsfd,J. J., "ANew Near-Wall Formulation forthe
veloc_ field at selected chordwise locations, for the iced
|wept wing. Bo_ the k-_results and the original algebraic k-t Equations of Turbulence," AIAA Paper 86.0556, 1986.
model resultsare shown. The gridused inthese twocases
was identical. It appears that the k-¢ model does I better 8. Khodadoust, A., "An Experimental Study of the
Jobof predicting the velocity profile, both near the surface Row Field on a Semi-Span Wing with a Simulated Glaze
end near the boundary layer edge. The grid used inthese Ice Accretion," Ph.D. Dissertation, University of Illinoisat
two alrnulations was rather coarse inthe normal direction. Urbena-CharnpaJgn, 1993.
Additional calculations on finer grids are needed, before
one can conclusively say which turbulence model is
performing better. Nevertheless, the improved agreement
wnh_ k-cmodel is encouraging.
CONCLUDING REMARKS
An existing compressible Navler-stokse solver for
dean and k_d wing analysis has been upgraded through
the use ofstate of the art upwind difference schemes, end
a two equation k-t model. Encouraging _ have been
4
Fig : 1Computational Grid (151X42X43) for iced wing with 300 sweep
angle
NI,Z _ g:lnBl"l'lanlCS
o.,O I • Lmaim
..--. J striP(
- -- f_TED t4WI
Fig 2 : Simulated Glaze ice accretion
3
yb = 0.42
2.52
1.5
Cp 1
0.5
0
-0.5 0.4 0.6 0.8 1
-1 X/C
Z.5
2 y/b = 0.72
1.5
I
0.5
0
O.Z
-0.5
13
-I
Z.5
y/b = 0.89
2
1.5
1
0.5
0
-0.5
-1
Figure 3. Surface pressure Distribution over a Clean Swept Wing at 8 Degree Angle of Attack,
computed using Roe Scheme with Baldwin-Lomax Model, _ Theory, D Expt.
6
3
y/b =0.42
2.5 k-e Model
2
1.5
Cp 1
0.5
0
-0.5
-1
y/b =0.72
k-e model
2.5
2 y/'o= 0.89
k-e model
1.5
1
0.5
0
-0.5
-1
Figure 4..Surface Pressure Distribution over a clean Swept Wing using the Roe scheme and the k-c
turbulence model, ....__ Theory, I'1 Expt
7
2
1.5
1
CpO.5
0
0.5
-0.5
X/C
-1
2
y/b = 0.72
1.5
1
0.5
0
-0.5
-1
2
y/b = 0.89
1.51!
0.5
0
0.5
-0.5
-1
Figure 5. Surface Pressure Distribution over an Iced Swept Wing using the Roe scheme and the k-c
turbulence model, _ Theory, n Expt
8
0.18 ,,,- a
O.Z x/c = 0.04 B 0 16 -J- x/c= 0.04, k-e Model a
O.15 Baldwin-Lomax Model 0114 .J.
0.1Z ..J. n
y/c " o.1
oos 006
• ,_/1_,___,__jd_ D 0.04? __
-0.5
0.5
0 1 1.5 -O.S 0 0.5 1
uAJINF
0.12
0.2 x/c = 0.18, k-e Model
x/c = 0.18 0.1
Baldwin-Lomax Model
0.15 0.08
0.06
0.1
0.04
0.05 0.0Z
0
0
0 0.5 1
0 0.5 I 1.5
0.14
0.Z
x/c = 0.3, k-e model
0.1Z
X/C= 0.3
Baldwin-Lomax Model 0.1
0.15
0.08
0.06
0.1
0.04
0.05 O.OZ
0
! |
-O.OZ 0.5 1
0 0 0.5 1
1.5
Figure 6. Velocity Profiles at a number of
Stations for the Iced Swept Wing at 4 degrees
Angle of Attack, _ Theory, D Expt
9