Numerical Simulations for Landing Gear Noise Generation and Radiation

Department of Aerospace Engineering

ThePennsylvania State University

FINAL REPORT "Numerical Simulations for Landing Gear Noise Generation and Radiation" NASA Grant 210NRA-98-LaRC-5 NASA Langley Research Center Submitted by Principal Investigators: Philip J. Morris Boeing/A.D. Welliver Professor of Aerospace Engineering Lyle N. Long Professor of Aerospace Engineering April 2002 CHAPTER 1: LANDING GEAR AERODYNAMIC NOISE PREDICTION USING UNSTRUCTURED GRIDS Frederic J. Souliez*, Lyle N. Long t, Philip J. Morris _ and Anupam Sharma ° Department of Aerospace Engineering The Pennsylvzmia State University University Park, PA 16802 Abstract Aerodynamic noise from a landing gear in a uniform flow is computed using the Ffowcs Williams-Hawkings (FW-H) equation. The time accurate flow data on the surface is obtained using a finite volume flow solver on an unstructured grid. The Ffowcs Williams-Hawkings equation is solved using surface integrals over the landing gear surface and over a permeable surface away from the landing gear. Two geometric configurations are tested in order to assess the impact of two lateral struts on the sound level and directivity in the far-field. Predictions from the F.fowcs Williams-Hawkings code are compared with direct calculations by the flow solver at several observer locations inside the computational domain. The permeable Ffowcs Williams-Hawkings surface predictions match those of the flow solver in the near-field. Far-field noise calculations coincide for both integration surfaces. The increase in drag observed between the two landing gear configurations is reflected in the sound pressure level and directivity mainly in the streamwise direction.

Nomenclature c speed of sound "Graduate Research Assistant, Student Member AIAA * Professor, Associate Fellow AIAA Boeing / A.D. Welliver Professor, Fellow AIAA -1- D Wheel diameter Heaviside function, H(/) = 0 for f< 0 and H(]') = 1 for f> 0 I-I(/) Li refer equation 2 M local Mach number vector of the source M

IMI

Mr unit normal vector to the surface compressive stress tensor P pressure freestream pressure P_ acoustic pressure (p - p= ) p" root mean square pressure perturbation P RMS r distance from source to observer unit normal vector from source to observer ret retarded (source) time SPL Sound Pressure Level (dB) Lighthill stress tensor r0.

t observer time U freestream velocity U_ refer equation 2 ith fluid velocity component ith surface velocity component of integration surface f = 0 Vi -2- observer location vector x

8(s) Dirac delta function, unity forf = O, zero otherwise

density P free-stream density Po retarded time T rn 2 wave operator 1/c 2 32/_t 2 - V z Introduction The Ffowcs Williams-Hawkings (FW-H) equation has recently been used with permeable surfaces in order to predict aerodynamic noise 1'2. It is an inexpensive method to include the quadrupole source terms inside the FW-H surface without performing any volume integrations. This can significantly improve the accuracy of the noise predictions at locations where nonlinear interactions in the flow cannot be ignored. This is notably the case with highly turbulent flows such as high Reynolds number jets and wakes. It is also only slightly more expensive to use than a moving Kirchhoff surface (see 0zy6riik and LongS), but without the limitations of Kirchhoff methods.

The motivation for predicting the far-field noise generated by a 4-wheel landing gear stems from the increasing contribution of airframe noise to the overall sound level of an aircraft in its landing approach. Early studies in the 1970's by Heller and Dobrzynski 4 showed that high-lift devices such as slats and flaps, as well as deployed gears, generated noise levels 10 dB higher than those of an aircraft in its "clean" cruise configuration.

Aerospatiale (now EADS Airbus) investigated the noise produced by several Airbus airplanes, which seems to indicate that noise from high-lift devices is likely to dominate -3-

for medium sizeaircraft,while landinggearnoiseseems more of a problem for existing

and future high capacity aircraft. The importance of investigating landing gear noise is reinforced by Airbus Industry plans to extend the Airbus family towards a high capacity aircraft.

Heller and Dobrzynski carried out a series of tests with both scale models 4 and full scale models 5, which underscored the lack of detailed geometric features with model- scale experiments and their effects on high frequency noise. These early experiments also showed that there is an increase in noise radiation from tandem axle configurations, which is the second test case in the present study. However, it was also found during the full scale experiment that struts, braces and other small features contribute significantly to the overall sound level. A more recent work by Dobrzynski et al 6 where the impact of various gear sizes and configurations is measured, illustrates the difficulty in using scale- model results for full-scale noise predictions. The actual simulation of the landing gear flow field is also of interest since it potentially affects the inflow of flaps located downstream. This was experimentally shown by Stoker et al 7 during a wind tunnel investigation of the airframe noise radiated by a model-scale Boeing 777, in which case a second high-frequency noise source from the flap system is only seen in the presence of the landing gear. This landing gear - flap interaction noise source was even shown to increase significantly by using a highly detailed gear geometry.

As already performed in a previous study by the same authors 8, the goal here is to combine the flexibility of unstructured grids with the FW-H equation. We use the Parallel Unstructured Maritime Aerodynamics (PUMA) code for generating the flow data. PUMA has been validated in several instances for simulating time-accurate flow -4-

data 9'_°.Theaimin thepresent caseis to evaluate theimpactonthenoisedirectivityand

intensityof two landinggeargeometries (LDG1andLDG2). It is expected to observe

largerpressure fluctuationsand a more complexthree-dimensional flow in the case

involvingtwo additional struts(LDG2).

The Computational Grids The grids used for the simulation of the flow over both landing gear configurations were generated using the commercial package Gridgen by Pointwise, Inc.

Figure 1 and Figure 2 show an overall view of the meshes on the landing gear surfaces with and without lateral struts respectively. The first mesh consists of about 80,000 surface triangles, for a total of about 880,000 tetrahedra in the volume mesh. The second mesh reused as much of the previous grid features as possible. With two additional struts, the number of triangles on the surface went up to 135,0001 with about 1.2 million tetrahedral cells. Specific attention was given to the cell clustering between the front and rear wheels, in order to capture as much of the wake from the upstream wheel impinging on the downstream wheel. Flow separation from the fore wheel and wake impingement on the aft are expected to generate large unsteady pressure fluctuations and therefore noise. With the second geometry, great care was given to the mesh refinement between the two lateral struts, with the aft strut in the wake of the fore strut. The smallest geometric features were not overly simplified, since they have been shown to generate high frequency noise as explained in a later section describing the Ffowcs Williams - Hawkings equation.

-5-

As shown in Figure 3 for the secondgeometry(LDG2), a porous FW-H

integration surfacewasusedin additionto the flow datacollectedon the landinggear

surfaces themselves. This will help determine the magnitude of the quadrupole source

termfor this low Mach numberflow. Permeable FW-H surfaces wereusedfor both

geometries, with about13,000 triangles in the first case, and15,500 trianglesin thecase

includingtwo lateralstruts. This coarsening meshawayfromthe solidsurfaceis dueto

computerlimitations. It may not be able to supportthe higher frequency pressure

fluctuations.Theadvantage of these porous surfaces is thattheycancapture quadrupole-

like termswithouthavingto performany volumeintegration.FW-H surfaces can be

used in regions dominated by nonlinear effects(unliketheK_irchhoff formulations).

The Gibbs-Poole-Stockmeyer algorithm tl was used to speed-up the

communication process between CPUs.As shown in Figure4, thisprocedure dividesthe

domain into slicesthatminimizethenumber of messages between each processor, sothat

eachCPU exchanges datawith at mosttwo neighboring CPUs. The time stepfor the

unsteady simulations is determined by the smallest cell characteristic length. At a CFL

number of 0.95, this yields a time step of 0.86E-08 second for the first grid, and 1.90E-08 second for the second grid (due mostly to some improved CAD work in the original geometry file). The numerical conditions were dictated by the CFL3D run perfomaed at Langley: the Reynolds number based on the wheel diameter is 1.25 million, for a free stream Mach number of 0.2. The actual wheel diameter of the landing gear model is 9.4 cm.

-6- Flow Solver PUMA (Parallel Unstructured Maritime Aerodynamics) is the computer program that was used to run the unsteady calculations. It is a finite volume, Runge-Kutta time- marching code that solves the compressible Navier-Stokes equations and uses unstructured grids. It uses the Message Passing Interface (MPI) library for parallel implementation. Its scaling performance for the two configurations is illustrated in figure 5. The flop performance is slightly higher for the second case for any given number of CPUs, since the ratio of computation over communication is greater for a larger grid.

The facility used to perform the computation is the latest our two Cost effective Computing Arrays (COCOA and COCOA2) x2. COCOA2 is a Beowulf cluster comprised of 20 nodes each having dual 800 MHz Pentium III and 1 GB RAM. The cluster has dual fast-Ethernet on each node and all the nodes are connected using two HP2524 switches with channel bonding for increased data communication. These machines run Redhat Linux (version 7.0) and the gcc compiler.

For these simulations the inherent artificial dissipation provided by Roe's flux integration scheme acts as a sub-grid scale turbulence model. A parallel investigation on separated flow around a cone 13 shows that the implementation of a Large Eddy Simulation (LES) method using a Smagorinsky sub-grid scale model 14 may improve PUMA's accuracy to simulate both mean and turbulent quantities in the wake of a cone base flow. LES has already been used extensively to compute sound sources 15'16, but some recent work related to two-time statistics of LES data, would indicate that LES fields are too coherent if the eddy viscosity model does not include any random -7-

backscatter Iv. Onewayto circumvent this maybe theuseof a dynamicLES,which is

morelikely to yield enough backscattering to decorrelate thefluid motionatlargescales.

An exampleof the use of a dynamicsubgridscalemodel combinedwith a Ffowcs

Williams- Hawkingssolveris givenby Morriset a118 in anattemptto simulatethejet

noiseforcircularnozzles.

Simulation Results

Simulations were carried out over two cycles based on the expected shedding frequency of the wheel diameter. Each simulation took about 90 days on 24 CPUs. It is worth mentioning that the existing amount of data for both gear configurations put a strain on the available capacity in terms of storage requirements: about 40 Gigabytes of data have been collected for the two calculations described in this study. Data were sampled only for the second cycle to minimize the effects of the starting conditions.

Local time stepping was initially used to accelerate the convergence from free-stream conditions to a realistic state. This is achieved by assigning to each cell the maximum allowable time step for a given CFL number (pseudo time marching). Global time stepping is then turned on before unsteady data is sampled. In order to evaluate the total drag, the momentum deficit method is used by evaluating the velocity deficit in the wake of the landing gear. More details can be found in Rae and Pope 19. Figure 6 shows the average velocity deficit right behind the second gear configuration. In the second gear case, there is a good qualitative agreement with experimental results published by Stoker 7 in the high-fidelity landing gear configuration. Figure 7 shows the drag forces computed by integrating the pressure on the gear surface (pressure drag) and using Pope's wake -8-

deficit approach (labeledtotal drag). Resultsare shownfor both gearconfigurations

duringapproximately 7 milliseconds of simulated flow time,givingenough time for the

fluid to coverthreetimesthe geartotal length. Theforcecoefficients (drag,lateraland

verticalforce)arethe computed forcesdividedby the landinggearsurface areaandthe

dynamic pressure.The increase in overalldragdueto the introductionof the lateral

support strutsis largesincethesecomponents arenot aerodynamically profiled andare

comparable to flat platesfacingtheincomingfluid flow. Figure8 illustrates the lateral

forces stemming fromthepresence of these twostruts.

Oneexpects the far-fieldsound pressure levelto reflectthis unsteady loadingin

both its intensityand directivity. Figure9, which is a displayof the instantaneous

distribution of pressure onthelandinggearsurface in its second configuration, showsthe

pressure on thesesupportstruts,as well as on the wheels. Figure 10 showsa 3D

representation of somevortex filamentsshedding off various gearcomponents, and

highlighttheimpactof the upstream elements' wakeontogearelements at downstream

locations.Theeffectof thesevorticesis not completely captured by the FW-H surface

whichliesonthelanding gearitself. However thepermeable FW-Hsurface does account

forall theeffects induced by thesefilaments untiltheycross itsboundaries.

Far-Field Noise Prediction Only recently has the FW-H equation been used on a permeable surface, di Francescantonio was able to show that simply integrating the surface source terms on a porous FW-H surface does account for the quadrupole sources enclosed within the -9-

surface. The FW-H equationis written in the standard differentialform includingall

quadrupole, dipoleandmonopole source termsas

a _ (1)

2p'(x,t): x_Oxj [T,,H(f)]:_---_[L, fi(f) J+_[ (poU.)fi(f) ] Where Li and U. are defined as (2) U. =U ih i U i= 1- v i+- Po The subscript n indicates the projection of a vector quantity in the surface normal direction. Using the solution to the above equation given in Brentner and Farassat 21 and neglecting the quadrupole terms, the pressure fluctuation at a given observer location x and time t is (equation 3 below)

4zrp'(x,t)= f F mr L (,- )

+ dS

I

c rO2M,) J L;O-M )

I=o L r2(1-Mr) 2 JJ_,, FW-H code validation In the absence (at the moment) of experimental acoustic data to compare to, an already-proven method was implemented: the use of the CFD results to validate the FW- H sound predictions. As was done by the authors in a previous test case s, the pressure fluctuations computed by the flow solver PUMA at an observer in the near field were -10-

compared with the predictions givenby theFW-Hpost-processing utility. Althoughthe

near-field pressurefluctuationsare large, and likely contain a great amount of

hydrodynamic oscillations, the derivation of the FW-Hequation is suchthatall pressure

perturbations (acoustic andhydrodynamic) shouldbe recovered.Examples in the near

field aregivenby Farassat andBrentner 22in the caseof high-speed impulsivenoiseat

rotorbladetip Machnumber closeto 0.9. It is assumed thatata Machnumber of 0.2,the

quadrupole termsdo not contribute significantly to thefar-fieldnoise. The solution po

to the quadrupole term of the FW-H equation is: 4_" p_(x,t) = 32 (4)

ax, ax--; } I

-_ f>0 r The volume integration, if performed, must be carried out over a large volume and represents a large computational task. The far field approximation of equation 4 reduces to:

(5)

4xp;(x,t)-1 2 2 i

co,2 I Trraad

-_ f>O r However, there is in the present case an interest in capturing quadrupole effects in the near field, so that an exact result to the FW-H is needed instead of the far field approximation. Farassat and Brentner 23 decomposed the quadrupole noise term into three components varying with l/r, 1/r 2 and 1/r 3 respectively: -11-

(6)

,. . 1_2 If 4,'rpQ (x")=c- _ J I T'_ dad'r -_f>o r 3 i 3T_ -T.

+ 3---; I. r 2 dad_ _f>O ' 3L -7;,,.

+ _f>O There is a possibility that the second and third terms may contribute in a significant way to the near field pressure variations. This implies that in order to validate the FW-H predictions against the CFD results one may have to account for some of these nonlinear effects in addition to loading noise in the near field since the observer is in a highly perturbed propagating medium. In the current derivation of the FW-H equation, the quadrupole term is not computed (to reduce computing time and to limit storage requirements). However the porous FW-H surface shown in a previous figure has the ability to recover all nonlinear effects occurring within its own boundaries. Figure 11 is an illustration of the instantaneous pressure distribution on the permeable FW-H integration surface.

Observers were placed just above the landing gear main leg (x -- 2.68 cm, y = 0 cm and z = 17 cm for observer whose pressure is depicted in Figure 12), where the porous FW-H mesh is more refined, so that the FW-H predictions can take place using both the solid and the porous FW-H surfaces. An example of comparisons with the PUMA results at one of these near-field observer locations is shown in Figure 12. This shows good agreement between the porous FW-H surface predictions and the solver computation, whereas the solid FW-H surface misses by more than 50% some of the -12-

pressure fluctuations. This demonstrates the ability of the secondFW-H surfaceto

predictthe entirepressure oscillations, either of acoustic or hydrodynamic nature. It

tends to suggest that quadrupole effects may represent a significant contribution to the overall near-field sound level even at moderate Mach numbers. Unless a volume integration is performed over the entire CFD domain, the entire pressure perturbation cannot be exactly reproduced where nonlinear effects are important and where vortices flow across the permeable FW-H surface.

As expected, the agreement between the predictions from the two FW-H surfaces improves in the far-field. Figure 13 shows the pressure time history at 40 radii from the landing gear at a 50 degree angle with respect to the downstream axis. The field produced by the coarser FW-H surface off the landing gear does not reflect the same high-frequency fluctuations given by the predictions coming from data collected on the gear itself. In view of experimental results described earlier, it was decided to use the solid FW-H surface to investigate the far-field noise directivity, where high-frequency signals are thought to be significant. Figure 14 below illustrates the decrease of the RMS pressure signal as one moves away from the landing gear along the downstream axis.

Calculations were made at 20 observers from 25 wheel diameters down to 50 wheel diameters in the wake of the landing gear. Both FW-H surface data were used and compared with a trend line assuming a signal decaying with 1/r. As observed previously, the agreement between the two surface predictions improves with increasing distance from the landing gear. As one moves further away from the landing gear, it is seen by the observer as an acoustic compact source, and the signal intensity should decrease with the inverse of the distance from the source.

-13- Sound directivity patterns The sound directivity in the medium- and far-field requires the use of a parallel version of the FW-H post-processing program 8'24. For each radii away from the landing gear, 72 observer locations are defined, so that a resolution of 5 degree angle is obtained.

A total of 648 observer points were defined, and are illustrated in Figure 15, which shows the relative scale with respect to the landing gear. For both gear configurations, all three orientation planes were studied and the results are reported in Figures 16a to 16f, Figures 17a to 17f and Figures 18a to 18f for the first and second gear configurations from a streamwise, spanwise and vertical perpective respectively. Polar directivity plots at radial locations of 10, 15 and 20 radii from the gear are plotted separately from the locations further away (30 to 50 radii from the gear) for scaling issues. Sound Pressure Level (SPL) contours with a reference pressure of 6x104 Pa are also presented for radial locations varying from 25 to 50 radii from the landing gear. The scale for equivalent configurations is unchanged in order to allow for qualitative comparisons with respect to both directivity and intensity of the sound pressure signal. In Figures 16 to 18, the landing gear is not to scale, and is meant to illustrate which orientation axis is shown.

The drag augmentation is reflected in the sound directivity patterns of both configurations. The intensity of the RMS pressure is greatly increased along the streamwise direction for the second gear case (LDG2). This is due to the two support struts on the aerodynamic profile of the landing gear. Regarding the lateral noise, the two struts seem to interfere with the build-up of sound, so that the signal in the spanwise direction is less than that in the clean configuration, where varying lateral forces on the - 14-

gearleg create pressure levelsin thefar-fieldcomparable to thosealongthe streamwise

direction. In both cases, the near-fieldpressure perturbations are dominatedby the

fluctuations in drag. Little noiseis generated in the verticaldirectionsincemostof the

lift anddragvariations aregenerated on the maingearleg. The overallpressure field

looksmuchmore disturbed in the second gearconfiguration, illustratingthe complex

three-dimensionality of thenoise-generating flow pattern.

Conclusion

The flow field aroundtwo landinggearconfigurations of increasing complexity

hasbeenassessed. Thewakedeficitobserved behindthe landinggearis very similarto

thatexperimentally measured on comparable configurations.A parallelversionof the

FfowcsWilliams- Hawkings equation hasbeenimplemented usinginexpensive Beowulf

clusters to extract near-andfar-fieldsound information.Bothsolidandpermeable FW-H

integration surfaces havebeenused. Excellent agreement hasbeenobtained in thenear

field between the porousFW-H surface predictions andthe CFD solverresultswhere

hydrodynamic fluctuations areexpected to dominate andareof greater magnitudes than

thosetypical of acousticsignals. More work is neededin order to show that the

discrepancy observed between the solid andporousFW-H surfaces in the near-fieldis

linkedto short-range quadrupole-like effectseventhough theproblemthatwasdealtwith

in the present caseis a relativelylow Machnumberflow for this kind of effectsto be

significant.

The comparison of acousticpredictions producedby the two FW-H surfaces

improves astheobserver locationis movedfurtherawayin thefar field. Theincrease in

-15-

dragstemming fromthelateralstrutsis reflected in thenoiselevelanddirectivity. There

is a significantincreasein soundintensityin the streamwise direction,whereasthe

disturbance caused by thesegearelements seems to interfere with thevortexshedding off

thegearlegandtheresulting lateralsound radiation.

References i Singer, B. A., Brentner, K. S., Lockard, D. P. and Lilley, G. M., "Simulation of Acoustic Scattering from a Trailing Edge", AIAA Paper 1999-0231, 37 th Aerospace Sciences Meeting and Exhibit", Reno, NV, 1999.

2 Singer, B. A., Lockard, D. P., Brentner, K. S., Khorrami, M. R., Berkman, M. E. and Choudhari, M., "Computational Aeroacoustic Analysis of Slat Trailing-Edge Flow", IAA Paper 1999-1802, 5th AIAA/CEAS Aeroacoustics Conference, Greater Seattle, WA, 1999.

30zy6rfik, Y. and Long, L. N., "A New Efficient Algorithm for Computational Aeroacoustics on Parallel Processors", Journal of Computational Physics, Vol. 125, pp.

135-149.

4 Heller, H. H. and Dobrzynski, W. M., "Sound Radiation From Aircraft Wheel- Well/Landing-Gear Configurations", J Aircrafi, vol. 14, no. 8, Aug. 1977, pp.768-774.

5 Dobrzynski, W. M. and Buchholz, H., "Full-Scale Noise Testing on Airbus Landing Gears in the German Dutch Wind Tunnel", AIAA Paper 97-1597.

6 Dobrzynski, W., Chow, L. C., Guion P. and Shiells, D., "A European Study On Landing Gear Airframe Noise Sources", AIAA Paper 2000-1971, 6 th AIAA/CEAS Aeroacoustics Conference and Exhibit, Lahaina, HI, 2000.

-16-

7 Stoker, R. W. andSen,R., "An Experimental Investigation of AirframeNoiseUsinga

Model-Scale Boing777",ALAAPaper2001-0987, 39thAerospace Sciences Meetingand

Exhibit",Reno, NV, 2001.

8Long,L. N., Souliez, F. andSharma, A., "Aerodynamic NoisePrediction UsingParallel

Methods on Unstructured Grids", AL4Ok Paper 2001-2196, 7 th AIAA/CEAS

Aeroacoustics Conference, Maastricht, The Netherlands, 2001.

9 Modi, A. and Long, L. N., "Unsteady Separated Flow Simulations using a Cluster of Workstations", AIAA Paper 2000-0272, 38 th Aerospace Sciences Meeting and Exhibit", Reno, NV, 2000.

10 Sharma, A. and Long, L. N., "Airwake Simulations on LPD17 Ship", AIAA Paper 2001-2589, 31 st AIAA Fluid Dynamics Conference and Exhibit, Anaheim, CA, 2001.

11 Duff, I. S., Erisman, A. M. and Reid, J. K., Direct Methods for Sparse Matrices, Oxford University Press, 1986.

12 Long, L. N. and Modi, A., "Turbulent Flow and Aeroacoustics Simulations using a Cluster of Workstations", Linux Revolution Conference, Champaign, IL, June 2001.

t3 Souliez, F. J., Parallel Methods for Computing Unsteady Separated Flows Around Complex Geometries, PhD thesis, Penn State University, 2002.

14 Smagorinsky, J., "General Circulation Experiments with the primitive equations", Monthly Weather Review No. 91, pp. 99-165, 1963.

15 Piomelli, U., "Large-Eddy Simulation: Achievements and Challenges", Progress in Aerospace Sciences 35 (1999), pp. 335-362.

16 Seror, C., Sagaut, P., Bailly, C. and Juve, D., "On the Radiated Noise Computed by Large-Eddy Simulation", Phys. Fluids 13 (2001), pp. 476-487.

-17-

17 He,G.,Rubinstein, R.andWang, L. P.,"Effectsof EddyViscosityon Time

Correlations in LargeEddySimulation", ICASEReport No. 2001-10.

_s Morris,P.J., Scheidegger, T. E.andLong.L. N., "JetNoiseSimulations for Circular

Nozzles",AIAA Paper 2000-2080, 6thAIAA/CEASAeroacoustics Conference and

Exhibit,Lahaina, HI, 2000.

t9Rae,W. H. andPope,A., "Low-Speed Wind TunnelTesting",second edition,Wiley-

Interscience, 1984.

2odi Francescantonio, P., "A NewBoundary IntegralFormulation for the Predictionof

Sound Radiation", Journal of Sound and Vibration (1997) 202(4), pp. 491-509.

21 Brentner, K. S., and Farassat, F., "An Analytical Comparison of the Acoustic Analogy and Kirchhoff Formulation for Moving Surfaces" AIAA Journal, Vol. 36, No .8, Aug.

1998, pp. 1379-1386.

22 Farassat, F. and Brentner, K. S., "Supersonic Quadmpole Noise Theory for High-Speed Helicopter Rotors", Journal of Sound and Vibration (1998) 218(3), pp. 481-500.

23 Farassat, F. and Brentner, K. S., "The Uses and Abuses of the Acoustic Analogy in Helicopter Rotor Noise Prediction", Journal of the American Helicopter Society (1988) 33, pp. 29-36.

24 Long, L. N. and Brentner, K. S., "Self-Scheduling Parallel Methods for Multiple Serial Codes with Application to WOPWOP", Paper 2000-0346, 38 th Aerospace Sciences Meeting and Exhibit", Reno, NV, 2000.

-18- Z Figure 1 Surface mesh of first landing gear conflguration(LDG 1) Z Figure 2 Surface mesh of second landing gear configuration (LDG2) Figure 3 Partial view of the porous FW-H surface around'LDG2 gear configuration Z : CPU # Figure 4 GPS partitioning on LDG1 landing gear surface _icross 16 processors 1750 /. _

J

= LDG1 1500 ± LDG2 _.,_" . LDG1 ideal _,f 1250 - " =E 75O 5O0 25O 5 10 15 20 25 30 35 40 No. of processors Figure 5 Parallel speed up on COCOA2 for both landing gear configurations Z -1.07 -0.78 -0.50 -0.21 (u-U)/U : Figure 6 Average velocity deficit in the wake of the LDG2 configuration 0.14 - Total Drag - LDG1 Total Drag - LDG2 0.12 ....... Pressure Drag - LDG1 .... Pressure Drag - LDG2 0.1 .08 O O o_ 0.06 m°_._-_°_. °_._o m,_°_'_'_'m° _°_ D 0.04 0.02 _i lililllJlIIlllltllllll]*=llllll_lll I1| 00 0.001 0.002 0.003 0.004 0.005 0.006 0.007 0.008 Simulated Time (sec) Time-history of drag force coefficients for landing gear configurations Figure 7 LDG 1 and LDG2 using pressure and wake methods 0.012 - 0.008 f,%, " * _ /

/ \ _ /_

t- 0.004 o 0 Fy- LDG1 .......... Fz - LDG1 ....... Fy- LDG2 u.O-0.004 Fz - LDG2 • • °_" ° -0.0081 °,_, ,/ -0.012 0.001 0.002 0.003 0.004 0.005 0.006 0.007 0.008 Simulated Time (see) Figure 8 Time-history of lateral (Fy) and vertical (Fz) pr_essure force coefficients for configurations LDG1 and LDG2 Z P/P 1.020 1.004 0.983 0.971 0.959 0.946 j 0.996 0.934 0.910 Figure 9 Instantaneous pressure distribution for LDG2 configuration .Q (sec") 25471 17300 11750 5O0 Figure 10 Instantaneous vorticity filaments for the LDG2 configuration P/P 1.006 0.999 0.987 0.980 0,974 0.968 0.962 _ 0.993 0.956 0.950 Figure 11 Instantaneous pressure distribution on the permeabie FW-H surface for the LDG2 configuration 50- 40 - FWH1 - FWH2 • 30_ _' " e_ • PUMA I

20 ; I !i, r.,_

'° _

.lO__o - "..o." _ t

-30 ilillllPlitlt iJl_liJiilt I I I I I j I I I Ij -400 0.001 0.002 0.003 0.004 0.005 0.006 O.u07 Simulated Time (sec) Comparison in the near-field of the FW-H predictions from integration Figure 12 surfaces 1 and 2 versus the PUMA solver results FWH1 0.6 FWH2 0.5 0.4 0.3 0.2 _'o.1 D.

8 0 a.

n' -0.1 -0.2 -0.3 I '_1 -0.4 -0.5 -0.6 I I I I I I 1 ' ' l i r I ' l • i,ltlili_ill 0.015 0.016 0.017 0.018 0.019 0.02 Time (see) Comparison of the FW-H predictions from integration surfaces 1 and 2 in Figure 13 the far-field 0.025 • FWH 1 O FWH2 0.02 o l/r trendline ....0.015 a.

G.

0.01 0.005 i t i t I I I t I ] i ] I I J i k _ I I t I I , I, _ _ _ I 020 -25 -30 -35 -40 -45 -50 X/D Figure 14 RMS pressure signal predictions from integrationsurfaces 1 and 2 and trend lines

I

20 o" "o °

E • ." • " ." .'-..-'.,,.'-_%. - ; ".• : •

_-.. •. • . f__... • . •

_°i "",'(( )).'" "

- ...:: .:...

-_o_"" • • " : k \v/.._ : : •. : •

F" ." ",- "....__./- : : . •

L eo °io _eee_,@ D° oe e & eg BI -40

, I

-50 I, l,,,,l,t -50 -40 -30 -20 -10 0 10 20 30 40 50 60 x/R Figure 15 Observer locations for sound directivitycalculations 0.15 015 ....... 10D ....... 10D .... 15D .... 15D -- 20D -- 20D 0.1 0.1 0.05 0.05 a.

i f '1

5 0

5 0

'_ ,_1 i o,.

",,.1 -0.05 -005 -0,1 -0.1 =llll'llllllllllll' '11' -0,15

.o15 g_ _'2 " o -0.1 01.1 -0.3 -0.2 -0.1 0 01

P.. (Pa) P ,(Pa) .......... 3OD .......... 30D . . . . 35D .... 35D 0.02 0.02 ....... 40D ....... 40D .,/ ........ _-,_ .... 450 .... 45D f _*"_ -- _ _ _ "_"_ -- -- - 50D -- 50D .._'i:..<_' . - _ ..._.._.-... _'_.,_, ...'-"2'_'..- 001 001 ,'" t.", ",_""----"_ ",'-._._ ""-...

._.."" _"_ • -:, .'T. ,, ""...

• "_5_?- _', i / / " " :/. ,/ _"_, ', A V! ,/ \,: g. o Zo _ \ ',\ ),':

? x\_ ,,\ /_ J ?

,,. _..,._ _,," :' -0.01 -0,01 -0.02 -0.02 , i .... I , , , , I , , , , I , , -0.02 -0,01 0 0.01 -0.02 -0.01 0 0.01 P_ (Pa) P_ (Pa) SPL (dB) $PL (dB) ' 500 5O 5O 48 8 48.8 " 47.5 475 46.2 46.2 45.0 45.0 50,0 43.7 437 42.5 425 41.2 25 ,112 40.0 40.0 38.7 -25 0 25 50 Y/D Figure 16 (top) Medium-field RMS pressure, (middle) far-field RMS pressure and (bottom) SPL contour from streamwise perspective for LDG 1 (left) and LDG2 (right) 0,15 ....... 10D ....... 10D 0.15 .... 15D .... 15D -- 20D _'_'_ s,/._, ." -- 20D 0,1 0.05 11 _" -- \" # % I ':_,. I" i lk D.

5 0

i q i',

5 0

a.

a, / kl / ,' t -0.05 41.05 \. / _. i',._ \. i/ %l -01 41.1 *%._ p..*" , I , , , _ I i i i i I i i i ; I .... i ' ' ,Ik,llllllllldl11111llll -015 -o.15' -0.3 -0.2 -O.1 0 01 _.3 _2 _M 0 0.1 P_ (Pa) P_. (Pa) .......... 30D .... 35D .... 35D 0.02 0.02 ....... 40D ....... 40D .... 45D .... 45D 50D i .''' ....... _.._ 50D O01 0.01 .,'"" / ,..'_,,,_ "'"_._-_..... _..

,. I//',_- ...--..._ - ,-..,,, ,..

:' I :,Y _61 i / l!]/Jl;li i fl _ "/ "l_', I A ! I" / " : A • . / [iil / _11_( g.o _o

- _,,t_',ili

'., t t:t _/'.',',

i

,,, ,, "_ \_. ' _. ---'--_.L; _..- "' -0.01 -0.01 -0.02 -0.02 , I , , _ i I _ _ i i l r i a r I i r , I , , , b I i i i _ l I i i l 1 1 I -0.02 -0,o_ o o.o_ •0.02 -0o_ o oo_ P_ (Pa) Pm (Pa) SPL (liB) SPL (dB) 5O 5O 48,6 47,5 •. :: . 4., 46_.

45.0 : :'-::7_;._£::-:-::. :4, 43.7 42.5 50.0 41,2 25 25 .... ..,,.?;_,_ _,_,, 40.0 .,:.,>::skit"" %,_,7;::, , 'ii;//?/" _ t'_, i 38.7 ;';'i;!;'lllilil W l!',i)'; ,, ((/ ' ,

:'"N W '

',' ,N, ' ',X"," , ! i r -25 -25 • - -2: ...... .:::-- y- .... I , , , I , _ , I , i , , I ,

-5O5o ' ' ' '._ .... ' .... "5.... 5"' %0 .2_ o _5 50

XID XID Figure 17 (top) Medium-field RMS pressure, (middle) far-field RMS pressure and (bottom) SPL contour from sideways perspective for LDGI (left) and LDG2 (right) 0.15 ....... 100 015 ....... 10D .... 150 .... 15D "_"_. -- 20D -- 200 0.1 01 X .\ I \ \\ 0.05 ,..\ /',: 0 05 \k '\ I D,.

5 0 5 0

la.

_'_ 1 _ _ _:_ [ . . o.

.0.05 I -0.05 h.#' I / 't '--,--,,, / -0.1 -0.1 / i,\ i " i "1 ./ i /i , i _ ,._, , i .... i .... i .... i -015

•0.15_,_.... g._ .... .o'1 .... ; .... o'_' '

-03 -0.z -01 o Ol P_ (Pa) P ,(Pa) .......... 300 .......... 30D 0.02 0.02 .... 45D .... 450 .._. 1 .... _. '\ -- 50D ._' - - - _ .._ \ _, ............... ::-L.-. ..........

0.01 0.01 // / \_'_ I i (t i _,i,/ , ,,i ..;->_._ _\. _ ; _/il/" / I .'" ,/ _l I / A / I ! _/ {It' / #.o #.o

? ?

', ', '-,',_ 2';/. / / // i ) • ,.. ,, -.,'. ,_,,'/ ,.

-0.01 -0.01 Iii/ :% \ -0.02 -0.02 , i .... i i i h L I i i L h I i i -0.02 -0,01 0 0.01 -0.02 -0.01 0 0.01 P ,(Pa) P_ (Pa) Figure 18 (top) Medium-field RMS pressure, (middle) far-field RMS pressure and (bottom) SPL contour seen from above for LDG1 (left) and LDG2 (right)

AERODYNAMIC NOISE PREDICTION

USING PARALLEL METHODS ON

UNSTRUCTURED GRIDS

Lyle N. Long: Frederic Souliez t and Anupam Sharma t Department of Aerospace Engineering The Pennsylvania State University, PA-16802 Aerodynamic noise from a cone in a uniform flow is computed using the Ffowcs Williams-Hawklngs (FW-H) equation. The time accurate flow data is obtained using a finite volume flow solver on an unstructured grid. The FW-H equation is solved for surface integrals over a permeable surface away from the cone. Predictions from the FW-H code are compared with direct calculations by the flow solver at a few observer locations inside the computational domain. A very good qualitative match is obtained. Sound directivity patterns in the azimuthal and in the longitudinal directions are presented. The FW-H code is also validated against a model problem of a monopole in a uniform mean flow.

Nomenclature U_ refer Eq. 2 U_g_ coefficient of pressure Un c_ sub-grid scale constant in Smagorinsky model Dn c, u,_ e sound speed in quiescent medium Ua components of local fluid velocity d base diameter of the cone ui averaged streamwise perturbation velocity Uavg L vortex shedding frequency (Hz) uidi H(]) Heaviside function, H(f) = 0 for f < 0 and un local normal velocity of the source surface H(f) -- l for f > O Vn Dirac delta function refer Eq. 2 5(f) Li 6(f) = 1 for f = 0, otherwise 5(f) = 0 LM L_M_ Kronecker delta function, L_i 5_j Lr 5ij = 1 for i = j, otherwise 5ij = 0 density of the fluid M local Mach number vector of the source P freestream density of the fluid M IM I p0 density perturbation, p - po Mid_ P' M, vorticity (s -1) Mo Uo/c a angular frequency of the monopole source Mr MiFi w O 2 wave operator, {212 _ (_b-& - _72) M,e_ D2 fi unit normal vector to the surface, ni compressive stress tensor Introduction with potSij subtracted ECENTLY, the Ffowcs Williams-Hawkings P pressure (FW-H) equation has been used with permeable freestream pressure po surfaces for predicting aerodynamic noise. The pl acoustic pressure, p - Po application of FW-H in this manner effectively allows P'rms root mean squared pressure perturbation for the inclusion of the quadrupole source terms inside ret retarted time the surface without performing volume integrations.

Lighthill stress tensor T,j This has significantly improved the accuracy of noise t observer time prediction for cases where the contribution from angular location of the observer nonlinear interactions in the flow cannot be ignored.

U averaged streamwise velocity This is typical of highly turbulent flows, for example, freestream velocity Uo, Uoo high Reynolds number jets and wakes.

The FW-H equation requires time accurate data on, "Professor, Assoc. Fellow, AIAA. lnl_psu.edu.

and in the volume inside the permeable surface. This tGraduate Research Assistant, Pennsylvania State Univer- sity data is usually obtained by solving the Euler/Navier- Copyright _) 2001 by Lyle N. Long, The Pennsylvania State Stokes equations accurately in time. Since the FW-H University. Published by the American Institute of Aeronautics and equation uses data from within the FW-H surface, the Astronautics, Inc. with permission.

1OF9 AMERICAN INSTITUTE OF AERONAUTICS AND ASTRONAUTICS PAPER 2001-2196 outer grid can be made coarse without much loss of accuracy. Unstructured grids provide great flexibility in distributing the grid in the domain, and hence can be used to cluster the cells inside the FW-H surface.

This feature can be exploited to significantly increase the computation speed while keeping almost the same accuracy in predicting aerodynamic noise. This will also permit the modeling of complex geometries such as helicopter fuselages, landing gear, and flaps.

The goal here is to test the combination of unstruc- tured grids with the FW-H equation in predicting the aerodynamic noise. The test case is chosen to be ,_2-,+ ......_-- ,_ ,.., ;._._.., the flow over a cone. A cone has sharp edges which fixes the separation point. This makes the flow fairly Reynold's number independent.

We use the Parallel Unstructured Maritime Aero- dynamics (PUMA) 1 code for generating the time- accurate flow data. PUMA has been validated for time-accurate computations. 2-4 The ultimate aim is to predict the airframe noise from complex geometries such as landing gear, slats, and flaps. This cone case Fig. 1 Overall view of the 280,000 cell mesh.

may be considered as a benchmark problem.

during the acoustic prediction procedure, one does not The Grid have to take into account any phenomenon occurring The grid used for the simulation of the flow over a outside the integration surface. The surface can also cone of vertex angle 60 ° was generated using Gridgen. cross regions dominated by nonlinear effects.

Figure 1 shows an overall view of the mesh consisting During a run, the faces (triangles in this case) would of approximately 280,000 tetrahedra. The clustering be identified and flagged on each CPU, so that face was done around the cone and in the wake region with data would be output at a prescribed sampling rate increasing cell size towards the outer boundaries of the (around 50 kHz in the present case): the sampling was done in such a'm_ner that one had at least 20 computational domain. The reason for using Gridgen comes from one interesting feature of this commer- data points per wavelength, the shortest wavelength cial software: arbitrary surfaces can be created around being 10 times that of the simulated shedding fre- the cone (one within the CFD domain boundaries and quency. To avoid any redundant data, faces shared the other being the CFD domain boundary) and are between two adjacent CPUs had to be identified at sources for the meshing algorithm. It is possible to the beginning of each run, so that the number of faces export separately any of these closed surfaces in a sep- whose data are output is identical to the number of arate file, providing a means to extract flow data on the triangles on the actual FW-H surface. The grid par- surface using a FW-H module that was added to the titioning being done dynamically each time a run is unstructured solver. The smallest cylinder was used as initialized, the global cell indexing changes from run to a porous FW-H surface. At the bounding faces of the run, making it necessary to run the above flagging pro- CFD domain, Riemann boundary conditions were as- cedure any time the program is restarted. This makes signed at each face center, hence minimizing reflections the routine independent of the number of CPUs being from the boundaries into the computational domain.

used. Figure 2 illustrates the regions on the surface The large cells in the far-field also help dissipate any shared between 8 processors using the Gibbs-Poole- reflections. A no-slip condition was used at the solid Stockmeyer reordering algorithm. 5 As expected, each surface, even though the boundary layer was not re- region is a neighbor to at most two other partitions, solved due to computer limitations.

minimizing the amount of inter-processor communica- tion.

By using a set of faces that are actually used by the flow solver during the computation, there is no addi- The time step needed for a time-accurate solution is tional work required to extract the data needed for determined by the smallest cell characteristic length.

This is estimated to be one third of the cell volume the far-field noise. This type of FW-H surface also re- flects the true mesh clustering present where the flow divided by the maximum face area. For the grid variables are locally being computed: there is no loss described above, this yields a time step of 9.45E-08 in accuracy due to the interpolation onto a surface second at a CFL number of 0.9. The shedding fre- whose refinement might not be that of the computa- quency found during the experimental investigation of tional grid. Since only the surface terms are evaluated the flow is 36 Hz, for a Strouhal number equal to 0.171.

20P9 AMERICAN INSTITUTE OF AERONAUTICS AND ASTRONAUTICS PAPER 2001-2196 powerful machines, such jobs may require days, or even months to give results. Parallel computing using Be- owulf clusters offers an inexpensive way to handle such time-consuming simulations in reasonable amount of time.

Three facilities offering parallel computational power at Penn State were used for the computations - COst effective COmputing Array (COCOA), 2 CO- COA2 and LionX. 6 COCOA is a Beowulf cluster com- prising of 25 machines each having dual 400 MHz Pentium II processor. This facility was assembled by the authors and their colleagues in the Department of Aerospace Engineering at Penn State. The machines are connected via fast-Ethernet network which can support up to 100 Mbps bandwidth. A single Baynet- works 24-port fast-Ethernet switch with a backplane bandwidth of 2.5 Gbps is used for the networking. All the processors are dedicated to run parallel jobs. The Fig. 2 Partitioning of the FW-H surface across 8 operating system is Red Hat Linux. Message Pass- processors.

ing Interface (MPI) is used for parallel programming and the Gnu C compiler is used for compiling PUMA.

The Strouhal number was defined based on the cone di- Details regarding setting up and benchmarking of CO- ameter as St = fsd/U_. The numerical simulation is COA may be obtained from Modi and Long 2 and performed at Mach 0.2 at standard atmospheric pres- COCOA's website. 7 sure and temperature conditions, with an increased COCOA was primarily set up to make parallel com- viscosity to match the experiment's Reynolds number puting facility readily available to the CFD group of (50,000). Scaling the Strouhal number to the simu- the Aerospace Engineering Department at Pennsylva- lation's Mach number yields a shedding frequency of nia State University. The total cost of the cluster was 230 Hz. The computation of a complete shedding cycle just $80,000 in the year 1998, when it was set up.

requires roughly 46,000 iterations.

Since then this facility has been intensively used for various CFD simulations. COCOA2 is a newly assem- The Flow Solver - PUMA bled Beowulf cluster at Penn State. It has 21 nodes PUMA is a computer program, written in C, for each having dual 800 MHz Pentium III processors and the analysis of internal and external non-reacting 1 GB RAM each. The cluster has dual fast-Ethernet compressible flows over arbitrary complex geometries.

per node and all the nodes are connected using two PUMA uses the Message Passing Interface (MPI) to HP2524 switches with channel bonding.

run the code in parallel. It can be run on arbitrary Figure 4 plots the parallel speedup for COCOA and number of processors with very good scaling perfor- COCOA2 (1 Mflop = one million floating point opera- mance. Several papers 2,4 detail the benchmarking of tions per second). Fairly good performance is obtained the performance, and validation of PUMA.

considering the small size of the problem. Figure 4 PUMA is based on finite volume methods and sup- shows tile reduction in the flop rate per processor as ports mixed topology unstructured grids composed of the grid points are distributed over a larger number tetrahedra, wedges, pyramids and hexahedra. The of processors. This trend is typical of Beowulf clus- code may be run to preserve time accuracy for un- ters as the ratio of computation over communication steady problems, or may be run using a pseudo- decreases.

unsteady formulation to enhance the convergence to LionX is also a Beowulf cluster with 32 machines the steady state. Primitive flow quantities are com- (each having dual 400 MHz Intel Xeon processors).

puted at the cell centers. The code can be restarted These machines are connected via Myricom Myrinet from any point of time at which the solution is avail- with wire speed 1.28 Gbps. LionX also uses Linux with able from previous computations. All flow variables MPI for parallel programming. Performance com- are stored with double precision, but may be optionally parison and benchmarking results for LionX can be stored as single precision to save memory and commu- obtained from its website. 6 nication time at the cost of reduced precision.

CFD Results Parallel Machines After initializing all variables to the freestream val- Computational Aeroacoustics (CAA) codes are usu- ues, local time stepping is used to accelerate the con- ally very computationally intensive. Even with very vergence towards a physically realistic flow. This is 3OF9 AMERICAN INSTITUTE OF AERONAUTICS AND ASTRONAUTICS PAPER 2001-2196 IOO0 // 9OO COCOA ,/" ....... COCOAidea! .,, 8OO COCO,_. / .... COCOA,?. ideal .,' " /// _ 601) i // / / j.-'" :_ 40o!

I 300 _ I _/ I I•I 200i % 5 10 15 20 25 No. of processors Fig. 3 Parallel speed up for COCOA and COCOA2 Fig. 5 40 _ _ COCOA

J

2O 15 A , L _ , _ , , L ......... 21 - NO, of processors Fig. 4 Fig. 6 Average streamlines over one shedding cy- cle • done by assigning to each cell the maximum allowable Figure 7 shows the averaged streamwise velocity time step for a given CFL number based on each cell profiles computed by the original flow solver and those characteristic cell length. Global time stepping is then computed by the same solver combined with an LES.

turned on for several cycles before data are sampled, to ensure that the data on the FW-H surface follow In all three c&_es the magnitude of the reverse flow velocity is under predicted when compared with ex- the equations of motion. Figure 5 illustrates the vor- perimental measurements. The predictions agree fairly ticity patterns in the wake of the cone, showing strong well with Calvert's data in terms of the length of the recirculation phenomena. The noise from this recir- recirculation zone. Past the stagnation point, the re- culation is predicted by the FW-H module. Figure 6 sults including LES modeling follow the experimental is the averaged streamline contour over one shedding curve more closely than those computed without any period, illustrating the axisymmetric bubble that was turbulence model• observed during Calvert's experimental study. 8 In order to validate the solution, multiple compar- Figure 8 illustrates the variation of the pressure co- isons were made between the simulation and the exper- efficient Cp along the wake centerline. In this case, imental measurements. A basic Smagorinksy sub-grid the LES having the largest sub-grid scale constant scale turbulence model 9 was added to the flow solver Cs greatly over-corrects the pressure drop in the nea_ wake of the cone. The LES using a Smagorinsky in order to improve the predictions, since a large-eddy constant of 0.10 matches the measured pressure data simulation should yield better turbulent quantities.

4OF9 AMERICAN INSTITUTE OF AERONAUTICS AND ASTRONAUTICS PAPER 2001-2196 0.9 16[ AA&AA&&& 0.8 15= ,,, ,, 14_ • A AA 0.7 numerical _..,_-_" A A - numerical-C==0.10 ._/_.-. ,A 0.6 numerical - C= = 0.25 _'.'_ • = ,&'& 0.5 .'°S' 0.4 / • \ ,'-., 0.3 / -i \i,., , ,- _', • experiment / _, , .

"; 9F ' !'., .'fiX,', 0.2 •-1 0.1 _._7r_ L/ ,, ,,-" /, "_'v_c_..,'_.._ _,'\ 1 _ 2 3 -0.1 .,,,. ,,. xld ; ,r "-, -0.2 4 _ f'* • experimental "'- -0.3 3 _'- / _ numerical _i ....... numerical - C= ffi0.10 -0.4 AA AA AAA A -0.5 /' ......... numerical - C= ffi 0.25 -0.6 0J- t i i i I i i r I I i I , , I l _ l i 0 1 2 3 xJd Fig. 7 Comparison of the averaged streamwise Fig. 9 Comparison of averaged streamwise pertur- velocity with experiments for the flow solver with bation velocity with experiments with and without and without LES. LES.

yields unsteady velocity values that are closest to the experimental data. The solution without LES was se- lected to try to predict the far-field noise. It also leads 0.1 to the conclusion that a more advanced turbulence 0 model (dynamic LES, Detached Eddy Simulation) is ..., ..... ., ....

needed to simulate such separated flows, as found in -0.1 Strelets. 14 .o° / -0.2 Far-Field Noise Prediction e_ .. ooo°° / L) -0.3 The two commonly Used methods for far-field aero- i_/// * experiment dynamic noise predictions use the Kirchhoff equation --0.4 or the Ffowcs Williams-Hawkings (FW-H) equation.

,A / ." _ numerical While the governing equation in the 'moving surface'

-. // .-..:=-- ..::::=,,- = 01, 0

-0.5 Kirchhoff formulation 15 is a convective wave equation, -0.6 the FW-H equation is an exact rearrangement of the continuity and the momentum equations into the form of an inhomogeneous wave equation. Therein lies the strength of the FW-H equation over the Kirchhoff for- mulation. The FW-H equation gives accurate results Fig. 8 Comparison of the Cp coefficient with ex- periment for the flow solver with and without LES.

even if the surface of integration lies in the nonlinear flow region. This is typically the case in jets and wakes very well until the stagnation point is reached. These when the nonlinear region extends to large distances results are consistent with those found in other re- downstream.

lated investigations, using either the k-e turbulence In the Kirchhoff formulation the source terms are model l° or the k-e-v 2 model, zl These simulations were assumed to be distributed over a fictitious surface in compared against a set of experiments 12 at a lower the flow. The nonlinear effects (nonlinear wave prop- Reynolds number (42,000). Madabhushi z3 also used agation and steepening; variations in the local sound an LES with as many as 850,000 mesh points, but speed; and noise generated by shocks, vorticity, and completely over-predicted the length of the recircula- turbulence in the flow field) happening within the tion zone.

Kirchhoff surface are captured by the surface integra- Figure 9 shows that the averaged streamwise per- tion terms, but the Kirchhoff formulation requires the turbation velocity is not well predicted using any of integration surface to be placed in a linear flow re- the sub-grid scale constants. With the grid coarsening gion (i.e. far away from the body). This is difficult in the far wake, the fluctuating velocities are damped to achieve as most computational grids are generated very rapidly as one goes away from the cone base. It with the concern of minimizing computations. Usually, is the flow solver without any turbulence model that a fine quality mesh is used near the body with increas- 5oPg AMER[CAN iNSTITUTE OF AERONAUTICS AND ASTRONAUTICS PAPER 2001-2196 le_37 FW-H IXedk:t_*_ * ing cell size towards the outer boundaries. Therefore, ilytic_l .....

the quality of the solution available in the linear flow region is generally bad. The FW-H equation, on the other hand, works fine even if the integration surface is in the nonlinear flow region. A detailed comparison of the Kirchhoff and FW-H formulations is provided in Brentner and Farassat. _6

. i!ii-

-2e-_8 The solution of the full FW-H equation requires the _= -4e-O8 _ ; evaluation of two surface integrals and one volume integral. The surface integrations correspond to the "thickness" noise (monopole) and the "loading" noise -8e-o8 (dipole). The volume integration corresponds to the -le-.O, quadrupole term which accounts for the nonlinearity 0.5 1 1.5 time 2(tn _1 2.5 3 35 in the flow. 15,17 Evaluating the volume integral can Fig. 10 Validation of the FW-H code against the be extremely computationally intensive and difficult analytical solution for a stationary monopole in a to implement. Fortunately, the quadrupole term can uniform mean flow.

be safely ignored for most subsonic flows as is the case in the present study.

The quadrupole term is ignored in the present formu- Only recently has the FW-H equation been used on lation. The integrations are performed on the FW-H a fictitious (i.e. not the same as the body) permeable surface at retarded time. Since the FW-H surface is integration surface Is - exactly like the Kirchhoff ap- fixed relative to the body (the cone) for this study, and proach, di Francescantonio is demonstrated that when the flow Mach number is constant, the fol!owing terms the FW-H approach is applied on a Kirchhoff-type in the above integrals are zero : U_ = Mr = 0. The surface, the quadrupole sources enclosed within the standard time binning technique discussed by OzySriik surface are accounted for by the surface sources. It and Long 19 is used for obtaining pressure at the ob- should be noted that the "thickness" noise and the server locations.

"loading" noise as obtained from solving FW-H equa- The FW-H code and its Validation tion do not have any physical significance if the surface The FW-H code is written in Fortran 90. The of integration is chosen to be permeable (fictitious).

code was tested for a model problem - a stationary However, when the integration surface coincides with the body, these terms provide a physical insight into monopole in a unifocm mean flow. The FW-H surface is chosen to be a box made up of rectangular panels.

the source of sound generation.

The FW-H equation is written in the standard dif- The analytical solution to the mode] problem is eval- ferential form as uated at the center of each panel to obtain the time history of the primitive variables on the FW-H sur- face. The prediction from the FW-H code (using the 0 2 analytical data on the surface as input) is then com- O_p'(x,t) - oz_ozj [T_jH(I)] (1) pared with the analytical pressure perturbation at a 0 [L_(/)] + 0 point outside the surface. Figure 10 compares at an Oz, N [(poU.),_(Y)] arbitrary point (300 m, 0, 0) the pressure perturba- where Li and Un are defined as tion predicted by the FW-H code and that obtained analytically for a stationary monopole source with an Un = Ui_i :: Ui = (1 - P)vi + pUipo amplitude of 0.01 Pascals and a frequency of 2.267 Hz placed in a uniform mean flow of 0.3 Mach number.

Li = Pijrfj + pui(u,_ - vn) (2) The analytical solution to this problem is : and Tq is the Lighthill stress tensor. The FW-H equa- tion can be solved using the formulation in Brentner and Farassat, 16 and the solution can be written in an e exp(iw_-.)

¢(x,t) integral form as 4rr [(x + Uo(r, - t)) 2 + y2 + z211/2 (4) M0(=+U0(_.-e)) 1 + [(=+Uo(r._t))2+_2+z2p/2 rel_ where T. is given by + Mox- [(x 2 + (1 - M02)(y 2 + z2)] -- LM 7-, =t+ (5) + c(1 -M02) 6OF9 AMERICAN INSTITUTE OF AERONAUTICS AND ASTRONAUTICS PAPER 2001-2196 2e-o?

FW-H _edictloll using unit_JctUfe_ surf=c_ _lFl_ * Point No.

Analytlca_ 1.5e-07 5!

_ - .

1e-07 i !i Table 1 Coordinates of the observer locations for comparing FW-H predictions against PUMA.

-5e_8 I -_e-O? !

ordinates of the points are tabulated in Table 1. The

!/ i

._.5e-O?

cone has a base diameter of 0.02 m and a vertex angle of 60 ° . The center of the base of the cone is at the ori- -2e.07 gin and the vertex points upstream (positive x). The 0.5 1 1.5 Ik'ne 21in s) 2.5 3 3.5 FW-H surface is a cylinder of radius 0.05 m and length Fig. 11 Comparison of the FW-H prediction us- 0.175 m, centered at the origin.

ing unstructured surface grid against the analytical Figure 12 compares the pressure fluctuations at the solution.

four points listed in Table 1. Note that the PUMA pressure predictions have been shifted up by 20 Pas- The unstructured grid over the cone is created such cals. This is relatively a very small amount, about that there is an unstructured cylindrical surface en- 0.02% of the mean pressure. We believe that this closed in the computational domain (Fig. 1). This under-prediction by PUMA may be due to the dis- surface is chosen to be the permeable FW-H surface.

The elements of the surface are faces of the tetrahedra, sipation caused by inadequate clustering of grid cells.

It may also be due to the small sample size, and we and therefore, triangles. Since these triangles are cho- plan to do ensemble averaging. Note that this error sen from the unstructured mesh, the area and normal is of the order of magnitude of pressure pertubations varies from element to element. This, however, is not predicted by the FW-H code at any point inside the a problem because the FW-H equation only requires FW-H surface, which should actually be zero. How- information on a closed surface; it does not depend on ever, the prediction by the FW-H code agrees very well the structure of the elements constituting the surface.

qualitatively with the PUMA solution.

Clustering of the surface elements is desired to increase the resolution of the sources. The FW-H surface used Sound Directivity ....

for the present computation is the inner cylinder in The directivity of the noise from the cone was ob- Fig. 1. This grid was used with the model problem of tained by calculating the root mean squared (r.m.s.)

stationary monopole in a uniform mean flow to test if pressure perturbation for one shedding cycle at differ- the unstructured grid poses any problems. A perfect ent observer locations in azimuthal and longitudinal match is observed between the FW-H prediction and directions. Since the calculation for one observer loca- the analytical solution (Fig. 11). The comparison is tion is completely independent of any other location, it made at an aribtrary point (300 m, 0, 0). This con- is a perfect problem to run in parallel. Long and Brent- firms that an unstructured-mesh surface can be used ner 2° suggested some self-scheduling parallel methods as a FW-H surface without any loss of accuracy. Note for multiple serial codes. However, no parallalization that the first few seconds where the FW-H prediction was done for the noise prediction results presented does not match the analytical solution is the time it here.

takes for the sound to reach the observer. This delay Figure 13 plots the directivity pattern in the az- is more in Fig. 11 than in Fig. 10 because the unstruc- imuthal direction on the plane x = -0.1 m, which is tured FW-H surface is very small and hence, farther right behind the base of the cone. The pattern in Fig.

away from the observer point than the structured sur- 13 is symmetric because of the symmetry of the cone face used for Fig. 10.

about its axis. Since the FW-H equation cannot pre- Results for the Cone dict the pressure fluctuation inside the FW-H surface, we can compute the noise only outside the FW-H sur- PUMA is used to obtain time accurate data (the face. Therefore, the directivity patterns are plotted in primitive flow variables) on the FW-H surface. One an annular region outside the FW-H surface.

complete shedding cycle of the simulation is used for Figure 14 plots the directivity pattern in the longi- far-field noise prediction. Pressure at a few points out- tudinal direction on the z = 0 plane. Since the noise side the FW-H surface (in the near field) is collected is caused by both turbulence and fluctuating surface to compare with the predictions of the FW-H code.

forces, the directivity shows several lobes.

Four points distributed in the azimuthal direction near A conventional polar directivity pattern in the lon- the base of the cone and very close to but outside the gitudinal direction (z = 0 plane) is plotted in Fig. 15 FW-H surface were chosen for comparison. The co- 7oF9 AMERICAN INSTITUTE OF AERONAUT]CS AND ASTRONAUTICS PAPER 2001-2196 (1) .,_ FWH prmdi'clion -- PUMA c_IculaliOns ......

-13D _ -135 //' " a.

g_ i -140 I -145 .150 0001 00015 0O02 00025 0 003 00035 00o4 ,_ (in secs) (2) _, FVVH pmdl'c_on -- PUMA ca{culaHonl ......

Fig. 13 Directivity of the noise in the azimuthal _ 415 direction behind the base of the cone (x = -0.1).

43O

{

0.001 0.0015 O.OO2 00O25 0003 00035 0 O04 t{ • (m ,_.,_]

(3)

-130 FWH pr_cl_=_ -- -135 PUMA_ -140 //'/" _ -145 _ -150 f .1e5 -170 Fig. 14 Directivity of the noise in the longitudinal -175 ',... ,'" direction in the plane z = 0.

O_OOl 00015 0002 00025 0003 00035 0 004 k{r,_ On ,ec, I for observers at 10 different radial locations (r = 0.15 (4) .,oo - 0.24 m). In Fig. 15, the cone is pointing to the right; FWH pnKli_n -- PUMA CI{OJIi_O_ ......

the radial distance from the origin is equal to the r.m.s.

pressure and the angle (theta) illustrates the location of the observer point in the domain.

-120 Conclusions -130 i Aerodynamic noise from a cone has been studied l -140 as a model problem to test the possibility of using unstructured grids for noise prediction from compli- cated bodies like landing gears, slats etc. A finite .160 volume flow solver, PUMA has been used to obtain time-accurate flow data on a permeable FW-H sur- i i 0_11 0_OO15 o_on2 0_0025 0 O03 O_OD_ 0 OO4 face. The FW-H code was validated against a model problem of a monopole in a uniform mean flow. The Comparison of pressure fluctuation, p-poo Fig. 12 predictions from the FW-H code have been compared as predicted by PUMA and FW-H code at various locations listed in Table 1. 8 OF AMERICAN INSTITUTE OF AERONAUTICS AND ASTRONAUTICS PAPER 2001-2196 501, ,. , ,, , q - . L References tBruner, C. W. S. and Waiters, R. W., "Parallelization of 30 i i . : / " / ',: the Euler Equations on Unstructured Grids," AIAA Paper 1997- 20 i i ! _" : " - 1894, 35th Aerospace Sciences Meeting, Jan. 1997.

2Modi, A. and Long, L. N., "Unsteady Separated Flow Sim- ulations using a Cluster of Workstations," Paper 2000-0272, 38 _'h Aerospace Sciences Meeting & Exhibit, Jan. 2000.

3Sharma, A. and Long, L. N., "Airwake Simulations on LPD "_[-'°I._, o ....... I':i ........ _i ........... I _'__:)/i ........

17 Ship," Paper 2001-2589, 31 st AIAA Fluid Dynamics Confer- -20 : i ence and Exhibit, Anaheim, California, 2001.

-30 : '. . -} .

4Modi, A., Unsteady Separated Flow Simulations using a -_o_ _\ .

Cluster of Workstations, M.S. dissertation, The Pennsylvania -50 ...... State University, Department of AerospaCe Engineering, May -300 -250 -200 -I_o -loo -_o o 50 1o 1999.

p' cose(Pa) 51. S. Duff, A. M. E. and Reid, J. K., Direct Methods /or Sparse Matrices, Oxford University Press, 1986.

Fig. 15 Polar plot of sound directivity in z = 0 Shttp: / /cac.psu.edu/beatnic /Cluster/Lionx /perf.

plane at a few radial locations.

7http://cocoa.ihpca.psu.edu.

SCalvert, J. Ft., "Experiments on the Low-Speed Flow Past at four observer locations in the near field with direct Cones," Journal of Fluid Mechanics, Vol. 27, 1967, pp. 73-289.

calculations from PUMA. Noise predictions are made °Smagorinsky, J., "General Circulation Experiments with the Primitive Equations," Monthly Weather Review, Vol. 91, for a period of one shedding cycle. The comparison 1963, pp. 99-165.

is fairly accurate with only a small D.C shift error.

l°s. H. Johansson, L. D. and Olsson, E., "Numerical Simu- The directivity patterns of the noise from the cone lation of Vortex Shedding Past Triangular Cylinders at High are plotted in azimuthal and longitudinal directions.

Reynolds Number Using a k-e Turbulence Model," Interna- tional Journal For Numerical Me_hods in Fluids, Vol. 16, 1993, The sound directivity pattern has been shown to be pp. 859-878.

fairly complicated due to the complex physics inside 11Durbin, P. A., "Separated Flow Computations with the k- the FW-H surface.

e-v 2 Model," AIAA Journal, Vol. 33, No. 4, 1995, pp. 659-670.

X2A. Sjunnesson, C. N. and Max, E., "LDA Measurements of Velocities and Turbulence in a Bluff Body Stabilized Flame," Laser Anemomatry, Vol. 3, 1991, pp. 83-90.

13R. K. Madabhushi, D. C. and Barber, T. J., "Unsteady Simulations of Turbulent Flow Behind a 2"?iangular Bluff Body," paper 97-3182, 33rd AIAA Joint Propulsion Conference and Ex- hibit, Seattle, WA, 1997.

14Strelets, M., "Detachecl Eddy Simulation of Massively Sep- arated Flows," AIAA Paper 2002-0879, AIAA Aerospace Sci- ences Meeting and Exhibit, R.eno, NV, 2001.

lSFarassat, F. and Myers, M. K., "An Anaysis of the Quadrupole Noise Source in High Speed Rotating Blades," Com- putational Acoustics - Scattering, Gaussian Beams, and Aeroa- coustics, Vol. 2, 1990, pp. 227-240.

16Brentner, K. S. and Farazsat, F., "An Analytical Com- parison of the Acoustic Analogy and Kirchoff Formulation for Moving Surfaces," AIAA Journal, Vol. 36, No. 8, Aug. 1998, pp. 1379-1386.

17Faxassat, F., "Quadrupole Source in Prediction of Noise of Rotating Blades-A New Source Description," AIAA Paper 1987-2675, 1987.

lSdi Francescantonio, P., "A New Boundary Integral Formu- lation for the Prediction of Sound Radiation," Journal of Sound and Vibration, Vol. 202, No. 4, 1997, pp. 491-509.

19OzySriik, Y. and Long, L. N., "A New Efficient Algorithm for Computational Aeroaeoustics on Parallel Procesors," Jour- nal of Computational Physics, Vol. 125, 1996, pp. 135-149.

2°Long, L. N. and Brentner, K. S., "Self-Scheduling Par- allel Methods for Multiple Serial Codes with Application to WOPWOP," AIAA Paper 2000-0346, 38 th Aerospace Sciences Meeting & Exhibit, Reno, NV, 2000.

9oF9 AMERICAN INSTITUTE OF AERONAUTICS AND ASTRONAUTICS PAPER 2001-2196 Numerical Simulation of Subsonic Inviscid Flow past a Cone Using High-Order Finite Difference Schemes Jae Wook Kim" and Philip J. Morris _ Department of Aerospace Engineering The Pennsylvania State University University Park, PA 16802 ,_ Post-Doctoral Research Associate; currently, Research Professor, Department of Aerospace Engineering, Korea Advanced Institute of Science and Technology, Taejon 305-701, Republic of Korea _/Boeing/A.D. Welliver Professor, Associate Fellow AIAA Abstract The wake-dominated unsteady subsonic flow of Mach number 0.2 past a cone of vertex an- gle 60 ° is calculated numerically using high-order finite difference schemes on structured grids.

The three-dimensional compressible Euler equations are solved to simulate an inviscid flow that exhibits large fluctuations of pressure and velocity due to the shedding of vortices behind the cone. An axisymmetric structured grid system is used, that is generated by rotating a two- dimensional grid plane around a centerline. The grid singularity at the centerline, where the Jacobian and some grid metrics approach infinity, is avoided by changing the form of the flux vectors in the Euler equations without any asymptotic assumption or simplification. Fourth- and sixth-order finite difference schemes are used for the evaluation of spatial derivatives, and a fourth-order Runge-Kutta scheme is used for marching the solution in time. The complex wake motions behind the cone are investigated by visualizing the vorticity field. The mean flow pattern and periodic phenomena are analyzed and compared with the experimental data. This demonstrates the accuracy of the present approach to further analyses of wake-dominated flows past axisymmetric bluff bodies.

I. Introduction Calvert [1] studied experimentally the low-speed flow patterns and the associated periodic phenomena with flow past a cone of various vertex angles. He noted that there was little pub- lished work on incompressible flow past an axisymmetric blunt-based body, except some ex- periments carried out between the 1930's and 1960's. These earlier studies were restricted to flow visualization at low Reynolds numbers. He also discussed the periodic phenomena associ- ated with such flows that are presumably related to some kind of regular vortex shedding: however, the actual wake pattern is unknown above a Reynolds number of a few hundred. Cal-

vert'smeasurements ofmean andfluctuating velocity, andmean pressure forvarious cone an-

glesshowthatthewakes areallsimilar. Thisproblem hasreceived little attention in eitherex-

perimental or computational studies duringthelastthreedecades. Recently, Long et aI. [2] per- formed a simulation of a cone wake flow using low-order finite volume methods on unstruc- tured grids. The wake pattern and periodic phenomena in subsonic flow past a blunt-based body are of importance in the design of airframe components such as landing gears, especially for the reduction of aerodynamic noise: but, the existing knowledge of their flow characteristics is still far from complete. This paper seeks to provide some understanding of these complex flows through the simulation of the wake generated by an axisymmetric blunt-based cone.

The present paper focuses on the application of high-order finite difference schemes and structured grids to the computation of a three-dimensional, unsteady, inviscid, compressible flow past a cone at subsonic speeds. The scope of the present study is limited to the large scales of the wake turbulence and the shedding vortices, rather than the fine scales of turbulence whose simulation would demand huge computational resources and time. Therefore, the com- pressible Euler equations are used in the present computation. Fourth- and sixth-order numeri- cal schemes, based on central finite differences, are used for the solution of the Euler equations.

These methods have been developed for high accuracy and high resolution in the direct compu- tation of unsteady flows and acoustics [3, 4]. In addition, an artificial dissipation model and time-dependent boundary conditions are combined with the high-order schemes for a long- time stable solution [5-7].

The solutions are produced in a structured grid system generated by rotating a two- dimensional grid plane around the axis of symmetry. The axisymmetric structured grid has a grid singularity at the centerline, where the Jacobian and some grid metrics approach infinity.

When a conventional form of the Euler equations is considered, this singularity makes it hard to calculate the flux variables at the centerline and the flux derivatives near the centerline when wide differencing stencils are used across the centerline in a generalized coordinate system.

Thispaper proposes a waytoavoidthecenterline singularity by changing theformoftheflux

vectors in the Eulerequations withoutanyasymptotic assumption or simplification. In this

way,allthefluxvariables andderivatives maybeevaluated withoutlossof accuracy. Thesin-

gularitytreatment is a greathelpin thecalculation of theflux derivatives nearthecenterline: however, theEulerequations arestill indeterminate atthecenterline. In thepresent approach, thesolutions themselves atthecenterline areinterpolated fromtheneighboring values.

Theproblem conditions for thepresent computation areafree-stream Machnumber of0.2

andcone vertex angleof60 °.Thepresent workprovides a visualization ofthewakeflow field

with itsshedding vortices, andananalysis ofthemean flowpatterns andtheperiodic phenom-

ena.Theresults arecompared to Calvert's experimental data[1]for validation of theaccuracy

ofthesolution. A spiralpatternof thewakeflow behind theconeis shownclearly througha

visualization ofthewakevorficityfield.Theinitiallocation oftheregular vortexshedding and

itsStrouhal number aresimulated correctly. These successful comparisons indicate thatthepre- sent methodology iscapable ofdescribing three-dimensional, unsteady, bluff-body wakeflows.

In thenextsection, thegoverning equations andthetreatment ofthecenterline singularity

aredescribed. Thehigh-order finitedifference schemes arethenintroduced. Thegrid,andini-

tial andboundary conditions arethendiscussed. Bothqualitative andquantitative resultsfor

thewake flowfieldandacomparison withexperiment arethenpresented.

II. Governing Equations and Centerline Singularity The governing equations are the three-dimensional compressible Euler equations. The flux vector form of the Euler equations transformed to the computational domain may be expressed in generalized coordinates as,

where theflux vectors in generalized coordinates mayberepresented as

(2)

The conservative variables and the flux vectors in Cartesian coordinates are given by P pu pu 2 +p puv Q= pv E=[ pvu [,F= p , pw pw2+p

L.w J

(m, + [_(¢e, +p)q L(p , + p)w

where the total energy per unit mass is defined as e,=p/[(y-1)p]+(u2+v2+w2)/2 and T = G/G is the ratio of specific heats, y = 1.4 in the present computation. In the Eq. (2), the transformation Jacobian, J and the grid metrics, _: ,-.., (_. are given by a = x¢(y,z¢ - ycz, )+ x,, (ycz¢ - y_z¢ )+ x¢ (ycz, - y,z¢)' (3) with, (4) rl x fly rl, =JIy¢z¢-y_z¢ zcx¢-zcx¢ x;y¢-xcY;l.

_, _y _,] Iy,_z;-y¢z,7 z,x¢-z¢x, x,yc-xcy,] _ (r (: [.yCz,-y_z¢ zCx,-z,x¢ xCy,-x,y_J In an axisymmetric structured grid, the Jacobian approaches infinity at the centerline and the denominator in Eq. (3) becomes zero. Therefore, all the grid metrics expressed in Eq. (4) are indeterminate at the centerline. However, one may easily show that the ratios of the grid met- rics to the Jacobian remain finite or zero at the centerline, even though the Jacobian and some grid metrics themselves have infinite values. This suggests a way to remove the centerline sin- gularity. The key is to choose a different way of expressing the Jacobian and the grid metrics.

New variables that replace the conventional Jacobian and the grid metrics are defined by

(5)

J" -7 = x¢(yoz¢ - y,G )+ G (ycz_ - y_z¢ )+ xc(y_z _ - y_z_),

and

(6) The new variables with a superscript asterisk defined in Eqs. (5) and (6) have finite or zero val- ues in the entire computational domain, including at the centerline. Using the new variables, the expressions for the flux vectors given by Eq. (2) become (_=J'Q, ]_=_;E+_>'.F+_:'G, /r=qjE+_/;F+,*G, (_=(;E+(;F+_'.'G. (7) In this manner, the centerline singularity is removed and all the flux variables at the centerline may be evaluated. Especially, this treatment of the singularity benefits the calculation of the flux derivatives in Eq. (1) near the centerline when using wide differencing stencils across the center- line. However, the Euler equations still cannot be solved at the centerline since J" is zero there and the vector of the conservative variables (_ vanishes from Eq. (7). This means that Eq. (1) becomes indeterminate and it is impossible to integrate the solutions in time. However, the so- lutions (conservative variables) at the centerline may be interpolated from the neighboring val- ues.

III. High-Order Finite Difference Schemes High-order finite difference schemes with high-resolution characteristics are used in the present computation on a structured grid. The main scheme is a pentadiagonal type of central compact finite difference scheme [3, 4]. It is a generalization of the seven-point stencil Pad_ scheme used on the interior nodes. It may be expressed as Zf,'__+a f/_. + f/+af,:, + ' 1 , Pf,+_=; Z_.(f,+. -f,_.) (8) m=l where f is an objective function for the flux variables and f' is its spatial derivative on the i- th node. The grid spacing h is a constant independent of the index i in the computational do- main where all the grid points are equally spaced. Equation (8) may be solved by inverting a pentadiagonal matrix. The matrix must be completed at the boundaries. Therefore, non-central or one-sided formulations other than Eq. (8) are needed on the boundary and the near- boundary nodes to complete the matrix. These may be expressed as f0'+_z0.,f'+fl0.2f2' = _-_a0._,(f., -f0) for i=0 (boundary node), (9) em*0 eZ,.ofd+f,+_z,.2ff+fl,,ff: = 1 _-_a,.,,(f_-f) for/=1, (10) 1 $ fl2.ofo'+Cq.,f_'+Y[+a2,_f3"+fl2,,]_=_ _-_a2.,.(f.-f2) for i=2. " " (11) ra_0 .r_2 The coefficients in Eqs. (8)-(11) are listed in Table 1. They are optimized, as described in Refs. 3 and 4, to achieve maximum resolution characteristics with fourth-order accuracy (second-order in Eq. (9) for numerical stability). The optimized fourth-order compact finite difference schemes are used to evaluate the flux derivatives in the mean flow and radial directions. On the other hand, the flux derivatives in the azimuthal direction are calculated with a conventional central finite difference scheme for a simple handling of the periodic condition across the branch cut.

The standard central finite difference scheme is given by f,= 45(f,, - f_, )- 9(f+2 - f-2 )+ (f÷3 - f-3 ) (12) 60h Equation (12) has sixth-order accuracy, which is the highest order of truncation for the given standard seven-point stencil. Combined with thehigh-order finitedifference schemes in space, theclassical fourth-order four-stage Runge-Kutta scheme is used for marching thesolutions in time.

Asmentioned in Section II,thesolutions (conservative variables) areobtained atthecenter-

linenot by solvingthegoverning equations but by interpolation fromtheneighboring solu-

tions. In thepresent work,a fourth-order interpolation is used atthecenterline. It maybewrit-

tenas

f0 - 13(d + f_,)+s(A + f_2)-5(A +f-3) (13) where the negative indexes mean the values in the opposite direction across the centerline. Ac- tually, there are as many sets of the centerline solutions interpolated by Eq. (13) as half the number of grid points in the azimuthal direction. A unique value of the centerline solution is finally acquired by averaging these values.

High-order schemes in space and time resolve a wider range of wavenumber or frequency than low-order methods. However, even the present schemes do not resolve the high wavenumber or frequency range effectively, an adaptive nonlinear artificial dissipation model [6] is also used to remove the unwanted numerical oscillations that may develop from the unre- solved range. The artificial dissipation model is implemented only at the last (fourth) stage of the Runge-Kutta scheme in order to minimize computational costs. In addition to the stringent requirements on the high-order and high-resolution numerical schemes, an accurate and robust calculation depends heavily on the suppression of any waves that may result from unwanted reflections on computational boundaries. The boundary conditions for a time-dependent prob- lem should be physically correct and numerically well posed. Generalized characteristic bound- ary conditions [7] are used as the time-dependent boundary conditions in the present computa- tion. Non-reflecting inflow/outflow and the inviscid wall conditions are imposed at the boundaries ofthecomputational domain andthecone surface respectively.

Theaccuracy of thehigh-order compact finitedifference schemes, theadaptive nonlinear

artificialdissipation modelandthe generalized characteristic boundary conditions hasbeen

validated throughthe previouspublications of Refs. 3 to 7. Theywereappliedto multi-

dimensional steady/unsteady andinviscid/viscous computations in the variousbenchmark

problems: linearwaveconvection in Ref. 3,acoustic radiation fromaxisymmetric bafflein Ref.

4,noise generation fromairfoils inRef.5,shock-sound interaction in atransonic nozzle in Ref.6,

andwakeflowandacoustic radiation froma circular cylinder in Ref.7.It wasshownthatthey

produced veryaccurate numerical solutions in comparison with analytic solutions andexperi-

mental data. They canbeused effectively inthepresent computation.

IV. Procedure and Results of Numerical Simulation

In this section, the numerical simulation of an unsteady inviscid compressible flow past a cone of vertex angle 60 ° in a subsonic speed of Mach number 0.2 is described. The simulation uses high-order finite difference schemes and the modified form of the flux vectors in the Euler equations, in order to eliminate the centerline singularity.

A. Grid System and Computation Procedure For the present problem, the computational domain consists of two blocks and rotating a two-dimensional grid plane around the centerline as shown in Figs. 1 and 2 produces an axi- symmetric grid system. The number of grid planes in the azimuthal direction is 48. The number of grid points is 101x50x48 in block I and 161x51x48 in block 2. This gives a total of 636,528 grid points. The grid points are clustered along the centerline and near the edge of the cone base.

The minimum grid size at the center of the base is Ax/D = 0.007, Ay/D = 0.007 and Az/D = 27o<0.007/48 where D is the base diameter.

The timestepsizeisdetermined bytheCFL(Courant-Friedrichs-Lewy) condition [8]which

isgiven by, Chd' (14) • .2 .2 I ._ .2 ]' _--minlu']+lw']+lw I+ c(_/_; +_;2 +_; +_rz= + _ +.;2 +_/C2 +4-;2+C, J where the Courant number C = 1.0 is used in the present computation and c is the local speed of sound. In the evaluation of the time step, the operator 'rain' in Eq. (14) does not include the cen- terline as the denominator approaches infinity there and the Euler equations are not solved on the centerline as remarked in the previous sections.

A steady-state flow field is first acquired by solving the two-dimensional axisymmetric Euler equations. This is used as the initial condition for the present three-dimensional computa- tion. Figure 3(a) shows the pressure contours and Fig. 3(b) showsth e Mach number contours.

There is a recirculation bubble in cone base region. Calvert [1] deduced this characteristic form from his experimental measurements. The calculation, using the OpenMP [9] parallel libraries, is performed on an IBM SP2 parallel computer at the Pennsylvania State University Center for Academic Computing [10]. Occupying four CPUs, the actual computation time is 4.45 sec- onds/iteration, which gives 7 microseconds/iteration/number of grid points. The computation performs 500,000 iterations before a fully developed wake flow is generated.

B. Initial Triggering of Vortex Shedding In the early stages of the computation, it is difficult and time-consuming to initiate vortex shedding in the inviscid flow without some excitation. To start the first vortex shedding, the ve- locity normal to the base of the cone is excited in the form

(15)

where eistheamplitude ofthefluctuation, u, is the free-stream velocity, r is the radial distance, R is the base radius, 0is the azimuthal angle, and f_ is the frequency of the excitation. In Eq. (15), the distribution of the excitation velocity has four lobes with alternating signs on the cone base, and the amplitude decreases rapidly to half its initial value in one time period. In the present computation, the initial amplitude is chosen as e = 0.05 and the frequency is f_ = 0.171 u./D.

This corresponds to the regular vortex shedding frequency or the Strouhal number observed by Calvert [1]. Soon after the excitation vanishes, a non-axisymmetric flow field develops and the first vortex shedding occurs.

C. Investigation of Flow Field The unsteady wake flow may be visualized by the instantaneous pressure and Mach num- ber contours shown in Figs. 4 and 5 at the final time step (non-dimensional time is t ° = u.t/D = 64.04). Figure 4 shows the longitudinal distribution in a side view, and Fig. 5 shows the circum- ferential distribution in a rear view. The instantaneous three-dimensional unsteady flow field is very different from the two-dimensional axisymmetric steady flow field shown in Fig. 2. The pressure contours show many irregular vortices generated from the edge of the cone base and the Mach number contours show strong flow fluctuations behind the cone base. The flow pat- tern appears to be random or chaotic unlike the regular K_m_n vortex street [11] behind a two- dimensional bluff body. It can be seen that some vortices pair and merge mutually in the down- stream region. This results in a large-scale motion of the far wake. However, the vortices very near the outflow boundary are dissipated non-physically because of the lack of grid resolution.

In order to investigate the unsteady motion of the wake structure, time-traced snapshots of vorticity magnitude iso-surfaces are presented in Fig. 6. Figure 6 shows some vortex rings near the cone base that have non-axisymmetric distorted shapes, which seems to induce the chaotic flow fluctuations in the near wake and the spontaneous vortex shedding downstream. It is also shown that the vortex rings in the near wake region interact each other and change their shapes into the axial vortex tubes in some transitional region away from the cone base at a distance nearly equal to the wall diameter (x/D -_ 1). Further downstream, the vortex tubes shed from the transitional region possess a velocity component in the circumferential direction. This re- sults in a spiral motion of the far wake structure. Moreover, the vortex tubes sometimes twist, bend and stretch in the far wake region. These phenomena occur quite randomly and it is al- most impossible to pick up any two snapshots with the same instantaneous flow field during the entire computation. It is difficult to discern any periodicity from these time-traced observa° tions of the flow field. In the next subsection, a frequency analysis of the flow signals is pre- sented to quantify the periodicity.

D. Analysis of Flow Signals .......

The unsteady pressure and axial velocity along the centerline behind the cone are obtained as a function of time in the fully developed wake stage. By averaging the signals in time, the ax- ial variations of the mean pressure and axial velocity may be obtained. They are compared with Calvert's experimental data [1] in Figs. 7 and 8. Figure 7 shows that the mean pressure coeffi- cient from the present computation agrees with the experimental data well in the near wake re- gion up to the location of its minimum value. It then recovers faster to the free-stream value with a smaller overshoot than the experimental data in the far wake region. Figure 8 shows that the mean velocity curve of the present computation is a little more positive than the experimen- tal data and recovers slightly faster to the free-stream value too. The agreement between the predicted pressure coefficient and mean axial velocity is much better than that achieved in the low-order, unstructured grid calculations of Long et al. [2]. The mean stagnation point, where R = 0 and Cp = 0, of the present computation is slightly closer to the cone base than that of the experiment. It is likely that the small discrepancy between the present results and the experi- ment data is caused by the difference between the inviscid flow and the real viscous flow. The Euler equations have no physical viscosity that dissipates out the flow kinetic energy as the wakes move downstream, and the inviscid flow passes over a relatively smaller body in its ef- fective size than the real viscous flow since there is no boundary layer displacement effect.

Two variables are defined in Ref. 1 to characterize the fluctuating components of the flow.

They are denoted by, 8-_xlO0 and _-_xlO0 where u'= u- _ is the fluctuation of the axial velocity about its mean value and u, = lul is the magnitude of the axial velocity. Thus, 8 is a measure of the absolute level of the velocity fluctuations and # represents the local turbulence intensity. The axial variations of 8 and ¢ are shown and compared with the experimental data in Figs. 9 and 10, respectively. Figure 9 shows that the present results agree well overall with the experimental data except for a slight over- shoot in the near wake region and a small underprediction in the far wake region. Figure 10 shows the same kind of discrepancy between the present results and the experimental data in terms of the shortening of the axial development as is explained for the mean flow results in Figs. 7 and 8 above. The predicted overall amplitude of the fluctuations agrees very well with the measurements. It is much better than the predictions by Long et al. [2] particularly in the far wake region. However, the present grid is much finer in that region and the order of accuracy of the present numerical scheme is higher. A comparison of Figs. 8 and 10 shows that the region of the highest local turbulence intensities is in the vicinity of the stagnation point. This is consis- tent with Calvert's [1] observations. However, this high value is associated with the lower mean velocity atthislocation ratherthananincrease in theabsolute value of thevelocity fluctuations.

In fact, asseen in Fig.9,there is aslightdecrease in theabsolute valueofthefluctuations in this region.

Thepointof themean pressure minimumin Fig.7 corresponds to thetransitional region

(x/D ,_ 1) where the axial vortex tubes start to form as described in the previous subsection. It coincides with the point of the mean velocity minimum where the highest speed of reversed flow on the centerline occurs as shown in Fig. 8. This point has significance in estimating the starting position of the periodic vortex shedding. Calvert [1] remarked that a periodic wake mo- tion first appears in the region of the mean pressure minimum and the periodicity presumably arises from the instability of the free shear layer in an adverse pressure gradient. He also men- tioned that the periodicity is most prominent in the region of the mean pressure maximum. In order to quantify the periodicity, the frequency spectrum of the axial velocity is obtained at the point of the mean pressure maximum (x/D = 2.0). This is shown in Fig. 11. The highest peak in Fig. 11 shows that there is a strong periodic wake motion at a non-dimensional frequency of .fD/u_ = 0.171. This is exactly the Strouhal number found by Calvert [1]. The existence of this peak is taken as further evidence that the present computation succeeds in simulating the vor- tex shedding process accurately.

V. Conclusions The subsonic inviscid wake flow past a cone has been simulated using high-order finite dif- ference schemes on structured grids. The present computation, performed in an axisymmetric structured grid system is achieved by changing the form of the flux vectors in the Euler equa- tions to remove the centerline singularity in the generalized coordinates. This approach makes it possible to investigate the complex wake flow field and to obtain the accurate values of the flow properties. It is shown that the vortex rings in the near wake change their shapes into axial vortex tubes in the far wake through a transitional region where periodic vortex shedding be- gins. A spiral motion of the wake is found. The mean flow pattern agrees well with the experi- mental data when account is taken of the likely differences between the inviscid and the viscous cases. It is confirmed that the point of the mean pressure minimum is the starting position of the periodic vortex shedding. A spectral analysis of this periodic phenomenon shows that the Strouhal number is 0.171. This is in exact agreement with the experimental observation. On the basis of the very good agreement between the predictions and experiment, it is expected that the present methodology may be used for further analysis of wake-dominated flows past an ax- isymmetric blunt-based body.

Acknowledgment This research was supported by NASA Langley Research Center under Grant NRA 98- LaRC-5. The authors would like to acknowledge Center for Academic Computing at The Penn- sylvania State University for providing the supercomputer resources.

References 1. Calvert, J. R., "Experiments on the Low-Speed Flow past Cones," Journal of Fluid Mechanics, Vol. 27, part 2, 1967, pp. 273-289.

2. Long, L. N., Souliez, F., and Sharma, A., "Aerodynamic Noise Prediction Using Parallel Methods on Unstructured Grids," AIAA Paper 2001-2196, May 2001.

3. Kim, J. W., and Lee, D. J., "Optimized Compact Finite Difference Schemes with Maximum Resolution," A/AA Journal, Vol. 34, No. 5, 1996, pp. 887-893.

4. Kim, J. W., and Lee, D. J., "Implementation of Boundary Conditions for Optimized High- Order Compact Schemes," Journal o/Computational Acoustics, Vol. 5, No. 2, 1997, pp. 177-191.

5. Lockard, D. P., and Morris, P. J., "Radiated Noise from Airfoils in Realistic Mean Flows," AIAA Journal, Vol. 36, No. 6, 1998, pp. 907-914.

6. Kim, J. W., and Lee, D. J., "Adaptive Nonlinear Artificial Dissipation Model for Computa- tional Aeroacoustics," AIAA Journal, Vol. 39, No. 5, 2001, pp. 810-818.

7. Kim, J. W., and Lee, D. J., "Generalized Characteristic Boundary Conditions for Computa- tional Aeroacoustics," AIAA Journal, Vol. 38, No. 11, 2000, pp. 2040-2049.

8. Hirsch, C., "The Von Neumann Method for Stability Analysis," Numerical Computation of In- ternal and External Flows, 1st ed., Vol. 1, John Wiley & Sons, New York, 1992, pp. 283-341.

, http://www.openmp.org 10. http://natasha.cac.psu.edu/beatnic/SPinfo/SP_features.php 11. Schlichting, H., "Outline of Boundary-Layer Theory," Boundary Layer Theory, 7th Ed., McGraw-Hill, New York, 1979, pp. 24-46.

Table Caption Table 1. List of optimized coefficients for compact finite difference schemes Figure Captions Fig. 1.

Diagram of computational domain Fig. 2. Grid mesh system: (a) entire view, and (b) zoomed view Fig. 3. Zoomed side view of axisymmetric steady-state flow field for the initial condition: (a) pressure (p/p_), and (b) Mach number contours Zoomed side view of instantaneous flow field at the final time step: (a) pressure (p/p_), Fig. 4.

and (b) Mach number contours Fig. 5. Zoomed rear view of instantaneous flow field at the final time step (x/D = 1.0 from the cone base): (a) pressure (p/p_), and (b) Mach number contours Time-traced snapshots of vorticity magnitude iso-surfaces (17 levels of [Fo[D/u. from Fig.6.

0.0 to 0.002) in the interval of time t_,_ - t = 0.642 D/u. between two successive pictures

Fig.7. Axial variation of mean pressure along the centerline behind the cone

Fig. 8. Axial variation of mean axial velocity along the centerline behind the cone

Fig. 9. Axial variation of velocity fluctuations in terms of free-stream velocity Fig. 10.

Axial variation of velocity fluctuations in terms of the local axial velocity Fig. 11.

Spectrum of axial velocity fluctuations at the point of mean pressure maximum J. W. Kim and P. J. Morris Equation (8) Equation (9) Equation (10) Equation (11) al = 0.6511278808920836 ao,1 =-3.061503488555582 al,o = -0.5401943305881343 a2,o = -0.1327404414078232 a2 = 0.2487500014377899 ao_ = 5.917946021057852 a m = 0.8952361063034303 a2A = -0.6819452549637237 a3 = 0.006144796612699781 ao_ = 0.4176795271056629 al,3 = 0.2553815577627246 a2,3 = 0.7109139355526556 o_ = 0.5775233202590945 _,1 = 5.870156099940824 al,4 = 0.007549029394582539 az4 = 0.2459462758541114 It = 0.08953895334666784 _,2 = 3.157271034936285 Oa,o = 0.1663921564068434 a2,s = 0.003965415751510620 a_,2 = 0.7162501763222718 - f12,0 = 0.03447751898726934 flt_ = 0.08619830787164529 a'2,1 = 0.4406854601950040 a,2_ = 0.6055509079866320 fl2,4 = 0.08141498512587530 Table I J. W. Kim and P. J. Morris Y

D

Block 2

D Block I

5D

5D

Figure I

J.w. KimandP.J.Morris

Y

(a)

Figure 2-(a)

J.W.KimandP.J.Morris

Figure 2-(b)

2O J. W. Kim and P. J. Morris Figure3-(a)

J.W.KimandP.J.Morris

Figure 3-(b)

J.W.KimandP.J.Morris

Figure 4-(a)

J.w. KimandP.J.Morris

Figure 4-(b)

J.W.KimandP.J.Morris

Figure 5-(a)

J.W. Kim and P. J. Morris

Figure 5-(b)

J.W.KimandP.J.Morris

t =t_

I

l=t2 t=t7 _-1 1-.- x z-!L-- x l=13 ,,_ G _ :_i -_ _ - t_ tl4 t = t_ _ t=t9 j L__ x z_L.._x z,.-L.--x t=t5 t = t,_ __. ,

]__L_ "_ _

Figure 6

J.W.KimandP.J.Morris

0 1 i-i-i ...... i.... , -i-- i-- ;- J--i--i--i ..... _-_- i.-, --i_ 4 o-!-o- -O-i_ _ - ,- - i-- ;- -i ..... i--i--:-- _-_..-1._-_r._-_._._._-.;--_r._;._.-_._1_._.t_t_..;_t._.1._.t.-_._:..i_-_.;_._r_r-._w._1._-T_-;. r-_l-- -:_--.r --.r--.:---1---i-.- 1 I-..:---.L---_,---, ..._...;... ;....;.. J....:...;....,... J....;.._L.._....,:..._.._,...;..-.:....t..._...,_...;... J...,.Q.;..__.._...I O 01-,-:- :-;--:-':- :- ,--r :-:-_-_-:--:-,-)'_'_-_-_ _ + ..... I I--

_:..,._ ....... :.,,.._. ....... : ......... :..=.__..:._ .......... ,.._....... _ _ ;_

I _.2 ,, ,,, 7_:_:_, , ,

I

F"r-:---T"-;"-:":'-,----;'"r"-,"":'-,-'r7__T"T ''rT',''''',''r'';'::-':]

E ::'-'"=-:'-::_"': : -"'-*:"":::-'"-;:""_- """:-_ 1

t

....... ......... ,...... ........ .............

_-"- ' " ........ __- -.'---."--: '4-_ ' - ' ' ' " / -_,'e I__'I_4.- :- -,- ":- 2--;- - r -r - r - _ -/_ __ _ -_-- :1.:._-i_i::_]l-_ T- -, - "i" -i- -:- -| .L_.L--L.._,...J....;. , ' .. J...J.__.;,...L...L...;... j....L..I_...L...L._.[..i....;....;....;....__..;...L...; _. _..,',__.L...; ....

4,---_---_.--,,-- 4----_-- --'--q--,.b - _ - - -_- --,. -- ;---,_-- - -'-- 4 --- ,- - -,n- -_4-- - -_.---_. ---i -- -_ -o--,,.i--4--'*,_,_.,_.i- _..b- -, _ ---,- --, - -- 4 --- 1 -0 5 gNrt?$:_ : !::i::?:: • I---;.--'__-4-+-"th_O-_-O_Fd---:.---4 ----_ ...... ;---;----_ ..... ;---=- ; .--I -.i---_----_----',---W, i.--i---4---b-4--4---;---.;,.--i---i.-. ;.--..L. 4--, ;---;----i----!---4---;---_----i----i---;---;---_---4--.; ---;---;---,:-.- -06 .................................. _ _ 0.0 0.4 0.8 1.2 1.6 2.0 2.4 2.8 3.2 3.6 Figure 7

J.W.Kim andP.J.Morris

n a _-_---- .... :--_---_,,,---- ..... :,---;-,-;---_ .... ',----'-,-'----_-9-V--:'- _--i-- -_..___rese_:_._`_-._:_.:_._:_`._._-_::.._ ....

-..:...-:..._._.F..:...T...7..;...:... T...:...F..v.._......._...:....:...:._.I..7._... T.-.:_..i....

;, ---i-.4.--'.-.;.-.i.-.-i--...._---4.---.--i-..;..-4..- ,L--i---_---4..--i-._:::i ,- 0 2 -'---'--' ....... '--'--_-=--'- -' ..... :--'- ....... -'__--'--'--'- -:--',-- -.----:--- .-------,----:----,- -,--,--- .----:---, --.-- :----;---,---_--_-- F-I ....

.... ;-..,;...4--.._.-.1--.-;-...'-..._...'....L..'..._,..'..._-...'--._...'- __.._.-----.;.--: ....

I ........ g*-.._..- ;-..d....;-,.._..._'...-_...'... _...L.. _...I - -_....; ........ ;..._..._-...'...L.-a... _...'...;...a.-.-=...;..-;..._.-- g-.-.;...;..- I-..:_...,...O...:..._...,...'...-...'..._.. ___ _....'.._.',...:_.. _ i-__---_-._-._---i-.,.-__ ' '-_-',:---i---t --.-- , ; • _ : ; _ , + , _ . . ; , ; ', ', _ : ', : , _ _ ', _ ; _ . _ ', ; ', I-":"'_'"':"-'--:"-',-'":-'-0 ";_'"v-":'--_"_--_"": ""v":"'_"':' "'_" "':"":'":"',-"'r'",-":-'",'"=--"':'":'"','" v"-:--:-"1 I---:---£--'---_---:----:----;----_---,'_'--- _--O..)tC. ;...-...', ...a.-..'...."-.- "..._.--.'...." .-J...-a--- '-...:.-. :-.. a..-.:--.. ; ..._--- a,.-_--. :.--',-..

-0 6 I-ii;-TTiiIT 7T-I -;--iT-Ti-T_TTITTTiTTT:,--i "i--l-"':,-:,";" "0.0 0.4 0.8 1.2 1.6 2.0 2.4 2.8 3.2 3.6 7'" D Figure 8

J. W. Kim and P. J. Morris

:....;....:..,.J....;...._..A....J....L... L__ :.. _.,...._,....:.._&..._.._ _._..L..._... L.. J....;..._... J....:....;....;....L...;.... _..-._..._... ;....:....;....

181..._.-I-d_-t--;--:.L..L-_'--:-.:_.L-L_I-._-.J ..... ;.. _ L _ t _ .; _ .;_ .;_ _4_ _ L _ i.. " _ J -_: _ _:_ - :.- _ _. _ ,% _ : _ ..;_ _ -.-i--- G-4-.- J-.--i--- "_._._._._i_._._._-_t_._.._.i_..i_._.._&_i_._._._i..-_-_._._._...._ _--._.-..i-. -i--.

L. -:----F -;* - -- :'--- ,I-- --;--- -1- ..... 2:- ::::::::::::::::::::::: :-_--,--;- -.;'-r ' '-_'_ . -,-_- "!--!--:-r-f-_-ff-- ..............

......

• E:::.-: ....... ,'-.--' -: ....... _-.-L-: ........:...: ..,' ........ "--L-.-' ....... :.--: .-..

/ V ..;.. _..A.. L...'... L...L...;...-'-. __....', ... :..-.L._._...;...'....'_..._.. ,_.-..:....___.

_N_; .... ;-- -;-_-:-_--: ..... --'--:-_-_--:----_-:------ 8 "'_ .... _-. -i----;..-,.--4--..b--4---_,----_---b--_-.--;---;.--;.--+---_--.;----i--..i-.- - , , , , ;-'-:-'-';*"T'_-;-'-:-'-';----_'---;---7---:-'-_---.,----;----:---:---_--'-,,---y--T-'-;--, • ..-:.-..;...

E -\:!iiiil

2 _L "_":- :_:-I .--%--4----;---_---4----;,-.-;--- 4----; ,--:.--" ..- 4----;-.-4---_---_--- ;---_'---'---_.---'----'---"---.;---.:- ..-;-.-4-..4-.--'. *.'-,-._...,-.,:.-..;.-.._..q

I

0 --{---f---i- -I----i----k-i---t--_---i-i--¢---i-,-- +---t---i-f---i--i 'ffi-Tiji-ti-f,-Iii-i- 0.0 0.4 0.8 1.2 1.6 2.0 2.4 2.8 3.2 3.6

x."D

Figure 9 J. W. Kim and P. J. Morris 9O -.. :... J.. +.L .+. i.. _3_ _.. i ... ;_+. J___ .L.. + .'.+ +L_. _,-. + .L ...d... 3..+ .i ... :....L... ;...l.. + J_...;... _... j..+ .L... i+._ __...L+ + L+_ ;... 3-...l__. J....L..._.._ ,--i.--q+--4-.-' ....... _- ---_.+.:---i.--i ........ +. _.J t.._ _L . : i J t i ' 1 " ' ' • + , , .... ; . . , , , ; , . ; - .__: ....... ,_._-+.-:---,.--:.--.-..-!.._!._.+:_...,____...+.+._...

" ,'.- ,',-':-',--_-- F " :'- • "_ ...... ,--+'-t-': - _..-+- -+--,"- ,_- _ -', - -i,---+---,-.-+ -." -_a,r_;crr_v. - _- ,-- _- -_-- i ---_'_:_: ........_ +-+ .......+-+_::::_o1-,---_-+ ........_ i _ ........; _ _ .......... _ i -- V - "; - ":- -'- -:- - .:- _.'- - .'- K.- r-- a- q--i --:- -'--L-: +/- -_-;- --_: I. -_- -' - ..:- -: ..... ,-,-,.,_,__,__,_,_,_,.,__ __,_ ..........._ ........+_+__.:_,,___,,,=:; ........_;.- _. _, _: 5O '_- P-!i- Y * .......... Y_"i i 'A; ....... , - _+--.< ...... ;........... + :- ......._ ......... ; 4O _-: ,-._0-. + - • - ..'k--:- -:- -iY- "_,- -.:" -._ - _ - • N- ;- -: ..... _ - ; - + - -' - -_- -:-- _ - ; - ; - +- "- - •-.:.---_--.:-.-;-.-:----', . , . _.._::_."_\"i___!_;_._._!_+>_!_'_`_i_i._4_;_i_:_ --_ . ;.-._.--'-_ _...,._.._...+ ._..._..._ ..+... _+.+,. _ ........ ; .' -' --d. .......... ;...: ..... L ....... I , I , _ T "C'+T'";--" 2O

77;i;i;i17);7777111i;7i;)i!;iiiiii71

r--_.--_----_--+:---_-+--i ...... -'----i----i---i----L-+i--4----_----L.-i----:-.- _ - -;----i ....... -'- : _ -; +- i -: : -'- i r+-r.-.r...i-...r.._r++.i...i.-r.--;`.-r..t-..i...i.-.T.+T..i...T..T.T..T...r..i...r_..i.i...i...-t..`f.-.i...'...r..i...i...-V..

0.0 0.4 0.8 1.2 1.6 2.0 2.4 2.8 3.2 36 X./ F__) Figure 10 J. W. Kim and P. J. Morris = ................ ;- ":'-':";'r.-: .................. :-.----i--!-+++ ........... "...... i'"'"-'-'""-!-i ............. _ ..... :-- :---'- : - ...... ; :. _ ; _ : : ',q ............ , : : : : : ', : O .... _...... ,, --" ,,--',, "-: "_-':--;-r ........... _...... _"+" ":" "" _, "+;" ","" --:-'1 ........... ' ...... f ....... "_" "J, - -," "_" ,"-+ ........... ;" ..... J.-- ++,'-- - "+--:- -,_-_- _ .......... • ' ...... "...:-..;..;.:'-.:..:.L ! t : t ', : : , ', ! ', ', _ I ! , : : _ : : : , ; ,, : ; : ;, ......... _ ................. T'FT_ ....... ":...... ;'+"F"".'":","";'r', ........... t ..... ;+...:...:..+.:..:.:.: ..... _ ..... _-,,-_._,-:_.,_ ..... ,,---_ -__ ,_:.,_,_,_: .... _,........ _ ,,,_, ......... _ _ ;_; _ __:_.

_ ', ; : ,, ........ _ ....

.......... T ...... ,_"" _,'- + t'" ":" I+"," Y _.......... ] ...... :"+':'-+!'" _+T-,,-T _........ ', ...... ;--- -:--+_-._--:- -_-; _........... :- ..... ;-- --;-- -:-- _--;-_-:-

+

.......... ;...... ;+",'-'."; _,--:- _- _ .......... ;...... ,'---+--;--t- +-;--:-; ........ !...... i----i---_- .':--_ 4-i-: ........... ; ..... _----;--,:-._..'.:-: .......... ; ...... _--.;---_--;.;--:-'-_ ................ i--+_ ............ _............. _--_--'_-:.i ....... : ; : : : ', : : .... "_, ..... ._--_, _-- _ ..... _--- ,_-_-+--_ .... _ ........ 3- L -'-:" _ ........ L _ _' _:__'__ _ : : :;i : ; ! : . :}:l : : : : ', :::, : : :: .......... ; ...... ,'--'y";--:-y':-_-r .......... 3...... ;'"': ..... !. _..,-4--, ....... '...... _.--. -:-. -._-.._..;.._.. ;.., ........... ;. ..... .;+_..;...:.. _+ .;+_.;.

.......... : ...... ;----;.-.:...'.'.i.C-.:.! .......... : ...... '...._ : ._...-.:._'J .... .: . , ;. , _, ', : : '. : ; : '.

.w.._ ', : : i ; ', : I ' '. " ', ', , .......... . "" "i't-i ........... _,..... :"'-:--_,--".'-::-: O .......... ,_ ...... :'-"_"':---,_-';":":-_ ................. ÷..... '-_ ........ _---.'---.' _.._._._a ........... , .... j_. . :_ t. : .,, G .......... : ...... _....,'.+.;.. • _ ,, ........ , ; : .......... ; ...... ?"+_--"";+_-:'_'-_ ........ ;- -;'_ ....... ' .... ' _.-;- J ........ ;-- i + :- ; :-: :.

.......... ; ...... ,_- +" _,'" ",_"; "_ "';" _,"_..... _'- _ ........ _--,_ .......... ,_..... _,----;.. J-. J..L-:.

O _" ;-- -r - _ -:- _, _,-:'= -i- T , - ----'---'-'-''-'-'-'- ..........i......i-i--i--i-i--,_i- _ _, ........... i.....I--,-L-_---"--LL_- ........... + ..... -....... i..... _---i---!--.i--_4.;

-

01,,_

B

..... "'" ,++'I ......... ,_ -" ":,-++{+-:'- ,- ,_

<

i i i i i i;i i i i ilii[i __i;T_ 10 .3 10 .2 10 "1 0.171 100 101

,fD/,++

Figure 11

The Discontinuous Galerkin Finite Element Method

A Brief Review of the Method and its Applications in

Aeroacoustics

Preetham Rao* Department of Aerospace Engineering The Pennsylvania State University University Park, PA 16802, U.S.A.

Abstract The preliminary implementation and testing of the quadrature free Discontinuous Galerkin Finite Element Method in one and two dimensions are discussed. With an objective of understanding the method, algorithms to solve the linear scalar advection equation with periodic boundary conditions using the method have been written and executed successfully. A qualitative comparison of Lagrange polynomials and simple monomials as choices for the basis function set is made. Results from the algorithms using both Lagrange Polynomials and simple monomials as basis functions are also presented.

1 Introduction to the Discontinuous Galerkin

Method

Convection dominated problems arise in applications as diverse as aeroacoustics, gas dynamics, meteorology, oceanography, turbulent flows, viscoelastic flows, magneto hydrodynamics and electromagnetism, among many others. Devising robust, accurate, and efficient numerical methods for solving these problems is not a trivial task for two reasons. The first is that the exact solution of nonlinear problems develop discontinu- ities after finite time, and the second is that these numerical methods might display complicated and erroneous solutions at these discontinuities. Thus, while developing these numerical methods, care must be taken to guarantee that the discontinuities of the approximate solutions are physically relevant, and that the appearance of discon- tinuities in the approximate solution does not induce spurious oscillations [1].

The Discontinuous Galerkin Finite Element Method (DG method) is one of those methods that are being developed to successfully address the issues mentioned above.

The original DG method was developed in 1973, but it is only recently that enhance- ment and evolution of the method are taking place. The method has been proven to be well suited for high order accurate large time scale simulations. An important dis- tinction between the DG method and the conventional finite element methods is that "Graduate Research Assistant

theresulting formulation is localto theelement andthesolution is not reconstructed

by looking at the neighboring elements. Thuseach element in the DGmethodcan

bethought of asanindependent entity,merely requiring theboundary datafromthe

surrounding elements. Since the DGmethod incorporates numerical fluxesanddis-

continuous elements, it canbeconsidered asa generalization offinitevolume methods.

Owingto its unique finiteelement nature, theDGmethod hasnumerous advantages

overclassical finitevolume andfiniteelement methods suchas[1,2]: • TheDGmethod canbeused to obtainuniformly highorderaccurate solutions.

• Themethod ishighlyparallelizable since theelements arediscontinuous andthe

mass matrixis blockdiagonal andreadilyinvertible.

• It is suitedfor complex geometries sinceit canbe used on unstructured grids.

Themethod hasalsobeen shown to beimmune to mesh discontinuities. It also

requires simple treatment oftheboundary conditions.

• It can easily handle adaptive strategies, since refining orcoarsening thegridcanbe

achieved withoutconsidering thecontinuity restrictions typicalto theconforming

elements.

• It hasseveral useful mathematical properties with regards to stabilityandcon-

vergence.

• The method is compact, since each element is independent. This compactness

allows fora structured andsimplified coding forthemethod.

• Themethod allows forheterogeneity intheelements, thatis,theorder ofaccuracy,

shape, andeventhe choice of governing equations canvary fromelement to

element .....

Although the DGmethod is lesssusceptible to theproblems that arecommon in

finitedifference schemes, it is not withoutafewweaknesses [2].Themethod hasbeen

recognized asexpensive, in termsof bothcomputational operation countandstorage

requirements. Although theoretically themethod canbeapplied to anelement ofany

shape, the requirement of numerical quadrature for the integrals in the formulation

hasrestricted theapplications to hexahedral andquadrilateral elements. Therecently

developed [3] quadrature freeDG method triesto resolve theseproblems by using

polynomial basis functions, theproduct ofwhich canbeintegrated withoutnumerical

quadrature.

The restof this reportdealswith the quadrature freeDG methodappliedto

the scalaradvection equation, with periodicboundary conditions.Results for one-

dimensional linearandanonlinear _lvection equation arepresented, followed by solu-

tionsof thetwo-dimensional linearadvection equation ona square domain.

2 Formulation of the DG Method

Let the solution in an arbitrary domain be governed by a conservation equation of the form _, + v. # = 0 (1) Let the domain be divided into smaller elements _t(x, y, z) that span the domain.

Let the solution in each element be approximated using an expansion of a basis set given by N+I

_a = _ b._j (2)

j=l

B _ {bk,1 < k < _(p, d) + 1} (a)

The number of basis functions N + 1 depends on the order of expansion denoted by p, and the number of dimensions, d. Application of the traditional Galerkin method to equation (1) using (2) and (3) gives, (4)

£ bk(u, +v r)da: 0

for k= l...N + l.

Using (2) and integrating by parts, equation (4) can be recast as,

(5)

j=l for k= 1...N+I.

where d_' is the normal surface area element (or the line element in 2D), and /_R is the Riemann flux vector through the surface element, which will be approximated using the flux values of the element and the neighboring element.

It is convenient to carry out the integrations in the formulation (5) in transformed coordinates. In a conventional finite element method using quadrature, the integrals are calculated on a mapped element (denoted by A), but the final equation is evaluated and assembled in the real coordinates of f_. Since there is no assembling of the elements involved in the DG method, by mapping even the variables uj into the transformed coordinates, one can get a compact form of the equation (5). If the transformation from ft to A is linear, the storage requirements will be minimized. This will be explained in detail below.

Let the coordinates in the transformed space be (, r/ and C'. Let the unknown variables in the transformed element A(_, rh C) be denoted as vj. Then, N+I for k = 1 N + 1, where J -- o(x,_,z) and J8 is the Jacobian of the surface (or line) """ o(_,n,O ' coordinate transformation.

In particular, considering a linear advection equation in 2D, the flux is given by N+I N+I

(7)

j----I j=l and (6) can be written as,

(8)

M,j \ Ot,I - t(_j_ + f._ bkPR "IJ.la_= o

where

(9) M_j = fzx b_bklJIdA and iq Kij = JA J-1Vbkb_lJIA (10) for k = 1...N+ 1,j = 1...N+ 1.

Due to the presence of the term j-l, the matrix Kij is made up of a sum of several matrices. The Riemann flux is approximated as a Lax-Friedrich's flux of the form where b_ is the flux of the element to the left of the edge, and vecFs is that of the element to the right of the edge. _ is the unit normal to the edge from left to right, and a is a smooth positive function.

3 Computational Aspects of the DG Formula-

tion

The number of operations and the amount of storage required for the calculation of the integrals in the above formulation are determined by: .....

1. the order of the coordinate transformation and 2. the choice of basis functions.

If the transformation is linear, then J is a matrix made up of constants and IJ I is independent of { or 77 • This means that a single mass matrix Mij and component stiffness matrices Kij can be used for all the elements, thus minimizing the required storage. This is accomplished by using arbitrary triangles with straight edges for and an equilateral triangle for A, with the origin at its centroid. This is displayed in Fig. 1.

x-xo =Jr _ ,Jr= al a2

{ } {} [ ]

Y - Yo 77 bl b2 where al = xa - x2, a2 = 1/sqrt3(2xl - x2 - x8) bl=y3 - y2, b2=l/sqrt3(2yl - Y2 - Y3) where xl...x3 and Yl...Y3 are as shown in the figure 1, and x0 and !/0 are the coordinates of the centroid of the element _. The transformation matrix Jr will be

l

t m y ./x\ x Figure 1: Sketch of the transformation from physical to canonical coordinates different from, but related, to the transformation Jacobian d which appears in Eqn.

(8). ( JT = J' ).

This linear transformation also greatly reduces the computational time requirement.

Since a single mass matrix is involved for all elements, it can be inverted prior to the main computation. Hence time marching for each element reduces from solving a system of equations to multiplying and adding matrices and vectors.

3.1 Choice of the basis functions The choice of the expansion basis set greatly influences the computation time. Re- searchers [4] have worked on various polynomial basis function sets varying from simple monomials to orthogonal Legendre functions. A comparison of only the conventional Lagrange polynomials and simple monomials is presented here.

For a problem involving nonlinear flux , the integrals in Mij and Kij will involve products of more than two basis functions in their integrand. Choosing the conven- tional Lagrange polynomials for the basis set will necessitate excessive computation using numerical quadrature at each time step. Even if the integrals are computed an- alytically prior to the main computation, a huge storage will be required. The newly developed quadrature free approach aims at solving this problem by choosing simple monomial expansions for the basis set, thus making the evaluation of integrals simpler, and the matrices M and K sparser. The monomial expansions will be of the form _irt/. (For e.g. a second order basis set can be B = {1, _, r_, _2, _, r/2}. However, there are two drawbacks of the quadrature free approach. The first is that the condition number of the mass matrix is very high when monomials are used as basis functions.

The second is that the evaluation of boundary integrals in Eqn. (8) will involve edge transformations. Hence, evaluating the boundary integrals in the quadrature free ap- proach will require relatively complex coding and more computational effort. But this increase in computation effort will be offset by the decrease owing to the simplified calculation of the mass and stiffness matrix integrals.

On the other hand, if the flux is linear, using either Lagrange polynomials or the monomials, the integrals in M_j and Kij can be calculated analytically and stored before the main computation, without much usage of the memory. The mass matrix computed using Lagrange polynomials will be denser, but with a relatively lesser condi-

tionnumber compared tothemass matrixcomputed using monomials. Thecalculation

oftheboundary integral will beeasier andfaster using the Lagrange polynomials.

4 Results From the Implementation of the DG

Method

Several programs have been written to implement and test the quadrature free DG method in one and two dimensions. The results from these programs for the linear scalar advection equation in one and two dimensions using both the Lagrange polyno- mials and monomials are presented. A nonlinear advection equation in one dimension is solved using the quadrature free approach. A periodic boundary condition is imple- mented for both one and two-dimensional calculations. Time marching is accomplished using Crank Nicolson's 2 nd order scheme, a Fourth-order Runge-Kutta explicit scheme, and the Total Variation Bounded Runge-Kutta three (TVBRK3) stage method.

4.1 One dimensional calculations The test equations are cqu cqu 0--/ + a_ = 0 (12) and 0u a0(u2/2) 0--_ + Oxx = 0 (13) on the domain 0 < x < 10 with periodic boundary conditions.

The solution to the linear problem (12) is obtained with two initial conditions. The first is a half sine wave over a unit interval, and the second is a normal shock. While using the monomials, the initial condition was expanded as a Taylor series about the center of each element. The solution obtained using Lagrange polynomials is shown in figures 2a(1 st order approximation polynomials with 2 nd order time marching) and 2b (2 nd order approximation polynomials with 4 th order time marching).

The programs for the implementation of the DG method in one dimension are written in Matlab. The interval from x = (1, 10) is divided into 100 parts: that is, the value of Ax = 0.1. It is observed that the time step required for stability decreases as the order of spatial approximation increases.

The calculations for the linear advection equation with monomials are shown in Fig. 3a, b and c. It was observed that the use of monomials lead to oscillations in the approximate solution where the exact solution had sharp discontinuities, as can be seen in Fig. 3b. Finally, the solution for the nonlinear Eqn. 13 is shown in Fig 3d. The initial condition was a sine wave, which transformed into an 'N'-wave after sufficient time.

_2 ..... i .... ; ; i T ' , : : : : _ : -.- ¢_,W_", 1 s Js'_a !

..... _" - • * - '[..... F..... :- ..... :-- - - -_ ..... t h,zz'w_JI _a .....",..; ....._ ..... :. .....i,._ ....._ ....."...._ ......

i; ,! ! i !' i ! ! : o6 J _l 'i i i i i ! i i D4 / _" ti _ i :1 t:: ! i i / / D2 o2 t-----}----- i ..... i ..... :-..... i"-t ..... _..... "-_ .....

.... d

i i i i i [ i i i o_.::iii ii _ i, I 3 3 4 5 e., 7 el g lO a b Figure 2: One dimensional calculations using Lagrange polynomials t_ .... J-_._.: .... L_._; .... :......: .... [....L.... : ....

OJ ,,, .... -_--.: .... .... "_.-.-...;-.--...i ....

D.| :, ! : ,, ,, : : : : :i .: ' : : ; ; : : O_ ..... :,-- ._.} ........ } "--'_.... {" .... : .... 4 .... 0,4 ,: ; T : : : : : * T . " .... ;'"_"":'"', .... : ....

I o_ i i i i i i i i i .02 0 I 2 "1 4 6 I T _ a b :/ : -', ,,: ." : : : : os_. .... L...Y...; .... _... L....; .... ; .... _.... :...

<:/'i :_ ! i i i i OJ ........ : .... L....,....'.....' .... L... '.L...: ....

,k-..!..: ...._ ...._.:_...: ....i.... i.... i-

.o_ .... L,.._....; .... :.... L...,.' .... : .... --.'_'_.'...

o_ ....... : .... '. ....... '---._ .... "... '......: ....

i i i i i i !._ _

I i i i i i' i i i ,?'," , .... _i.... i .... i .... i...:-_::i:...i;...i .....

] : ; : : .* ; : : : 0 I 3 3 S I 7 _ _O • c d Figure 3: One dimensional calculations using monomials. (a) 3 rd order approximation mono- mials with 3rd order time marching, (b) Absolute error in the calculation, (c) Calculations with an initial condition of a normal shock wave, (d) Solution for the nonlinear advection equation.

=' ! i i ! i ! _ _ i "1.4 4t, ..41J 0 4 _ I II_ ,p.4 n_+ _lp a b Figure 4: The initial condition. (a) The sine wave pulse with unit maximum value at the center of the domain. (b) a contour map of the initial condition.

4.2 Two dimensional calculations 4.2.1 The advection equation The test equation is (]4) O"-U-U + a_x + b OU =0 Ot Oy u(0, x, y) = [sin(_rx) sin(Try)] 4 on a periodic square domain -1 < x, y < 1 . The grid is made up of 800 similar right angle triangles. The program is written in Fortran 90. The basis set is made up of Lagrange polynomials of up to second order and monomials of up to 4 th order. The initial condition is shown in Fig. 4.

Figure 5 shows the results for the case with the basis set made up of linear monomi- als. A similar result is obtained with all other choices of basis functions. The speeds of the wave in the x and y directions are taken to be the same. It is observed that the size of the time step required for a stable solution decreases as the order of approximation is increased.

4.3 Calculations for the linearized Euler equations in two dimensions Consider the propagation of an acoustic pulse in a constant mean flow. The lin- earized Euler equations are given by, c3U OE OF (15)

o--T + + N =o

with

U Mzu + P/Po , F = u ,E = M_u v Mzv Muv + P/Po

/

p Mxp + pou M_p + pov where p, u, v, p are the flow variables, and Mz and My are the values of the mean flow in the x and y directions respectively, non-dimensionalized by the speed of sound based on the mean thermodynamic values, and P0 is the mean flow density.

The finite element discretization of the above differential equation results in a sys- tem of equations that can be written as, M 0 0 OU 0 K22 0 K24 F2 (16) F3 0 0 0 0 /(33 K34[ U--

ii000

0 0 M 0 /(42 /(43 K44J F4 where M and Kij are the mass and stiffness matrices, and F1... F4 are the right hand side vectors of the corresponding individual equations.

4.3.1 Results The domain is divided with unstructured triangulation. Basis sets of degree 1 and 2 (order 2 and 3 respectively) are used. A square doniain is considered, with a uniform mean flow from the left to right, over a rigid bottom wall. An initial Gaussian distributed acoustic pulse is located near the lower boundary, as shown in Fig 6.

Reflective boundary conditions are imposed on the lower boundary, and non-reflecting boundary conditions are imposed on the three other boundaries. The radiating bound- ary conditions are implemented using characteristic variables [5]. Figures 7a and b show the computed propagation of the acoustic disturbance with the mean flow. Some reflections are observed near the right boundary, and these are due to the characteristic method of boundary condition used, as mentioned in [5].

5 Conclusions and Future Work

The linear advection equation with periodic boundary conditions has been solved using the DG method in one and two dimensions. A qualitative comparison of the Lagrange polynomials and simple monomials as choices for the basis set has been made.

A basis set comprised of simple monomials is more suited for h-p refinement of the solution. The choice of monomials for nonlinear convective equations is thought to be advantageous compared to the Lagrange polynomials under the formulation described, in which, the nonlinear flux is evaluated by multiplying or dividing the approximating expressions for the basic variables. However, the formulation involving the nonlinear flux could also be solved using iterative techniques. It remains to be investigated